StrainPhlAn
学习资料
介绍
StrainPhlAn 使用的是目标菌种中的特异性标记基因库 (MetaPhlAn)，StrainPhlAn 则 侧重利用 SNP 信息计算样品间的进化关系，并不反馈各样品内的具体菌株丰度。StrainPhlAn 有一个固定假设：一个样品中一个菌种只有一个 代表菌株，因此该代表菌株的丰度等于该菌种丰度。(谭 et al., 2020, p. 2614) (pdf)
一个例子
StrainPhlAn是一种用于在大量样本中跟踪个体菌株的计算工具。StrainPhlAn的输入是一组宏基因组样本，对于每个物种，输出是从样本中直接重建的所有物种菌株的多序列比对（multiple sequence alignment, MSA）文件。根据这个MSA，StrainPhlAn调用PhyloPhlAn来构建显示样本菌株进化的系统发育树。
对于每个样本，StrainPhlAn能够通过合并和连接映射到物种特异（species-specific）的MetaPhlAn标记的所有读取来识别该物种的菌株
详细来说，让我们从一个玩具例子开始，包括6个HMP肠道宏基因组样本（SRS055982-subjectID_638754422，SRS022137-subjectID_638754422，SRS019161-subjectID_763496533，SRS013951-subjectID_763496533，SRS014613-subjectID_763840445，SRS064276-subjectID_763840445）来自3个受试者（每个受试者在两个时间点采样），以及一个粪便拟杆菌(Bacteroides caccae)基因组G000273725。我们想要：
- 从这些样本中提取粪便拟杆菌菌株（G000273725），并将它们与参考基因组进行比较，构建一个系统发育树。
- 了解这些菌株与参考基因组之间的SNP数目。
通过在这些样本上运行StrainPhlAn，我们将获得粪便拟杆菌菌株（SGB1877）的系统发育树以及其多序列比对，如下图所示（使用ete2和Jalview生成）：
我们可以看到来自同一受试者的菌株被分组在一起。该树还突显了每个受试者的菌株在两个采样时间点之间没有发生进化。从树中我们还可以得知，主体“ SRS019161-subjectID_763496533,”的菌株与参考基因组更为接近，而其他主体的菌株与参考基因组的距离较远。
此外，下表显示了基于StrainPhlAn返回的菌株比对结果，样本菌株与参考基因组之间的SNP数目。
在接下来的部分，我们将逐步说明如何在这个玩具例子中运行StrainPhlAn，以重现上述图形。
- 从这里下载6个HMP肠道宏基因组样本、metadata.txt文件和一个参考基因组，将这些文件夹放置在"strainer_tutorial"文件夹下的"fastqs"和"reference_genomes"文件夹中。
步骤2. 通过将这些样本与MetaPhlAn数据库进行比对，获得它们的SAM文件
这一步将运行MetaPhlAn：对宏基因组样本进行与MetaPhlAn标记数据库的比对，并生成SAM文件（*.sam.bz2）。每个SAM文件（以SAM格式）包含针对MetaPhlAn标记数据库的映射读取。运行的命令如下：
mkdir -p sams
mkdir -p bowtie2
mkdir -p profiles
for f in fastq/SRS*
do
echo "Running MetaPhlAn on ${f}"
bn=$(basename ${f})
metaphlan ${f} --input_type fastq -s sams/${bn}.sam.bz2 --bowtie2out bowtie2/${bn}.bowtie2.bz2 -o profiles/${bn}_profiled.tsv
done
完成这一步骤后，您将在"sams"文件夹中找到包含SAM文件（*.sam.bz2）和其他MetaPhlAn输出文件的文件夹，这些文件夹包括"bowtie2"和"profiles"文件夹。
#mpa_vJan21_CHOCOPhlAnSGB_202103
#/ssd1/wy/.conda/envs/metaphlan/bin/metaphlan /ssd1/wy/workspace/metagenomics/strainphlan/raw_data/fastqs/SRS013951.fastq --input_type fastq --nproc 20 -s sams/SRS013951.fastq.sam.bz2 --bowtie2out bowtie2/SRS013951.fastq.bowtie2.bz2 -o profiles/SRS013951.fastq_profiled.tsv
#12847 reads processed
#SampleID Metaphlan_Analysis
#clade_name NCBI_tax_id relative_abundance additional_species
k__Bacteria 2 100.0
k__Bacteria|p__Bacteroidetes 2|976 100.0
k__Bacteria|p__Bacteroidetes|c__Bacteroidia 2|976|200643 100.0
k__Bacteria|p__Bacteroidetes|c__Bacteroidia|o__Bacteroidales 2|976|200643|171549 100.0
k__Bacteria|p__Bacteroidetes|c__Bacteroidia|o__Bacteroidales|f__Bacteroidaceae 2|976|200643|171549|815 100.0
k__Bacteria|p__Bacteroidetes|c__Bacteroidia|o__Bacteroidales|f__Bacteroidaceae|g__Bacteroides 2|976|200643|171549|815|816 100.0
k__Bacteria|p__Bacteroidetes|c__Bacteroidia|o__Bacteroidales|f__Bacteroidaceae|g__Bacteroides|s__Bacteroides_caccae 2|976|200643|171549|815|816|47678 100.0
k__Bacteria|p__Bacteroidetes|c__Bacteroidia|o__Bacteroidales|f__Bacteroidaceae|g__Bacteroides|s__Bacteroides_caccae|t__SGB1877 2|976|200643|171549|815|816|47678| 100.0
步骤3. 生成consensus-marker，这些文件是StrainPhlAn的输入
对于每个样本，这一步将重建其中发现的所有物种菌株，并将它们存储在一个pickle文件（*.pkl）中。这些菌株被称为样本重建的菌株。运行的命令如下：
mkdir -p consensus_markers
sample2markers.py -i sams/*.sam.bz2 -o consensus_markers -n 8
完成这一步后，您将拥有一个名为consensus_markers的文件夹，其中包含所有的样本标记文件（*.pkl）。
步骤4. 从MetaPhlAn数据库中提取Bacteroides_caccae（SGB1877）的标记（以便稍后添加其参考基因组）
这一步将从数据库中提取Bacteroides_caccae（SGB1877）的标记，然后StrainPhlAn将使用blast在参考基因组中识别最接近这些标记的序列（在下一步中）。这些序列将被连接起来，并称为参考基因组重建的菌株。运行的命令如下：
mkdir -p db_markers
extract_markers.py -c t__SGB1877 -o db_markers/
完成这一步后，您应该在"db_markers"文件夹中拥有一个文件："t__SGB1877.fna"，其中包含Bacteroides caccae（SGB1877）的标记。
步骤5. 构建多序列比对和系统发育树
这一步将根据在步骤3生成的标记文件中的样本重建菌株（sample-reconstructed strains）以及参考基因组中的存在情况（如果指定了参考基因组），对选定的类群标记进行筛选。同时，将根据选定的类群标记的存在情况，对样本重建菌株和参考基因组进行筛选。从筛选后的标记和样本中，StrainPhlAn将调用PhyloPhlAn来生成多序列比对（MSA），然后构建系统发育树。
mkdir -p output
strainphlan -s consensus_markers/*.pkl -m db_markers/t__SGB1877.fna -r reference_genomes/G000273725.fna.bz2 -o output -n 8 -c t__SGB1877 --mutation_rates
完成这一步后，您将在"output/RAxML_bestTree.t__SGB1877.StrainPhlAn4.tre"中找到生成的树。
为了添加metadata，我们提供了一个名为add_metadata_tree.py的脚本，可以按照以下方式使用：
add_metadata_tree.py -t output/RAxML_bestTree.t__SGB1877.StrainPhlAn4.tre -f fastq/metadata.txt -m subjectID --string_to_remove .fastq.bz2
脚本"add_metadata_tree.py"可以接受多个元数据文件（用空格分隔，也可以使用通配符），以及多个树文件。元数据文件是一个以制表符分隔的文件，其中第一行是元数据的标题，接下来的行包含每个样本的元数据。多个元数据文件在以下情况下使用：如果您的样本来自多个数据集，并且不想合并元数据文件。