MetaPhlAn丰度估计

学习资料

• Customizing the database
• MetaPhlAn 3.0 Tutorial（包含六个测试数据）
• metaphlan_databases
• 完整代码

markers2reads, n_metagenome_reads, avg_read_length = map2bbh(pars['inp'], pars['min_mapq_val'], pars['input_type'], pars['min_alignment_len'], pars['nreads'], pars['subsampling'], pars['subsampling_seed'])

使用isolate and metagenome-assembled genomes从头构建数据库

ref

为了生成将在 MetaPhlAn 4 中分析的 SGBs 数据库，分离基因组部分包括 NCBI提供的 236,620个细菌和原生动物基因组，并标注为 "reconstructed from isolate sequencing or single cells"。这些基因组与从人类（5 个不同的人体主要部位，164 个不同的人类队列）、动物宿主（包括 22 个非人灵长类物种）和非宿主相关环境（包括土壤、淡水和海洋；补充表 2 和 3）采集的样本中组装的 771,528 个 MAGs 集成在一起。

在剔除不符合质量控制标准的参考基因组和 MAGs（即基因组完整性高于 50%，污染低于 5%）后，该目录包括 729 195 个基因组（560 084 个 MAGs 和 169 111个参考基因组），并按 5%的序列相似性使用 Mash聚类为 SGB，最终形成 70,927 个 SGB 数据库，其中 47,190个在物种水平上分类不明（uSGBs）。该目录涵盖 95 个不同的门类，uSGBs 在这些门类中的富集程度相当一致。与最初的 SGBref目录相比，目前的收集所整合的来自高度多样化环境的 MAGs 多了 3.6 倍,所定义的 SGBs 多了 4.3 倍。虽然该资源库可用于基于基因组的更大规模研究，其规模远超迄今为止所描述的规模，但我们在此将重点放在元基因组分类群的鉴定和量化任务上。为此，并为了降低没有强力支持或极其罕见的 SGB 的潜在假阳性检测率，我们只保留了来自不同样本的至少含有五个 MAG 的 uSGB，以进行后续的元基因组剖析，最终得到了 29.4 千个质量受控的 SGB。

在 493,482 个新的 MAG 和参考基因组上应用 PhyloPhlAn 3（参考文献 81）的 "phylophlan_metagenomic "子程序，以确定它们最接近的 SGB、GGB 和 FGB 及其 MASH 距离。根据报告的距离，我们按照 Pasolli 等人定义的阈值（分别为 5%、15% 和 30% 的遗传距离） 将基因组分配到已有的 SGB、GGB 和 FGB 中。然后，我们使用 1.1.25 版的 python 软件包 "fastcluster"，对未分配到任何现有 SGB 的基因组的全对全 MASH 距离进行了平均关联的分层聚类。以 5%、15% 和 30% 的遗传距离为分界点，对得到的树枝图进行划分，分别定义了 54,596 个新的 SGB、37,546 个新的 GGB 和 18,211 个新的 FGB。

我们根据 NCBI 分类数据库（截至 2021 年 2 月）49 为所有 70,927 个 SGB 分配了分类标签。对于 kSGB，我们通过对每个 SGB 中包含的参考基因组的分类标签应用多数原则来分配分类标签。如果出现平局，则通过选择具有代表性的分类群，即字母排序在前的分类群来解决分类标签问题。对于 uSGB，我们采用了类似的多数原则，但对 GGB 级别所包含的参考基因组的分类学进行了处理，分配了直至属级的分类学标签。如果在 GGB 级没有参考基因组，我们就进一步在 FGB 级应用相同的程序。如果在 FGB 水平上没有找到参考基因组，我们就只分配到门一级的分类标签，方法是考虑在最接近的参考基因组的分类标签集中最经常出现的门，以及最多 100 个与 "phylophlan_metagenomic "确定的最接近的参考基因组的基因组距离在 5%以内的参考基因组。对于没有收到任何分类标签的分类级别，我们为所有内部分类节点分配了 SGB、GGB 和 FGB 标识符，以维持分类的所有级别，并对 uSGB 进行分类。

然后使用 UniRef90 数据库对基因组目录进行注释，并使用 UniClust90（参考文献 59）标准（同源性大于 90%，聚类中心点覆盖率大于 80%）对每个 SGB 中无法归入 UniRef90 基因家族的基因进行重新聚类。利用由此产生的 UniRef 和 UniClust90 注释，我们为每个质量受控的 SGB 定义了一组核心基因（存在于组成 SGB 的几乎所有基因组中的基因），在将所有核心基因与整个基因组目录进行映射后，我们为总共 21,978 个 kSGB 和 4,992 个 uSGB 定义了一组 510 万个 SGB 特异性标记基因（不存在于任何其他 SGB 中的核心基因）,平均每个 SGB 有 189 ± 34 个独特的标记基因。

对筛选出的 729 195 个 MAGs 和参考基因组进行了注释工作流程：(1)用 Prokka（1.14 版）86 处理 FASTA 文件，以检测和注释编码序列（CDS）；(2)随后用基于 DIAMOND 的管道（可在 https://github.com/biobakery/uniref_annotator 上查阅）将 CDS 分配到 UniRef90 集群58。基于 DIAMOND 的管道根据 UniRef90 数据库（2019_06 版）对蛋白质序列进行序列搜索（DIAMOND 0.9.24 版）87 ，然后对映射结果应用 UniRef90 的纳入标准来注释输入序列（>90% 的同一性和>80% 的聚类中心点覆盖率）。在每个 SGB 中，未被归入任何 UniRef90 聚类的蛋白质序列按照 Uniclust90 标准（"-c 0.80-min-seq-id 0.9 "参数）59 使用 MMseqs2（参考文献 88）进行聚类。

对于每个 SGB，根据 UniRef90 和 UniClust90 注释，收集该 SGB 至少一个基因组中存在的所有 UniRef/UniClust90 簇，生成一个泛基因组。对于每个聚类，在所有基因组中随机选取代表性序列，并根据该 SGB 的 2 k 个最高质量基因组中聚类的普遍程度计算出核心度值。在接下来的步骤中，包含少于五个 MAG 的 uSGB 将被丢弃。我们之所以实施这一限制，是因为我们发现有证据表明，一些小的 uSGB 可能包含了组装假象或嵌合基因组，而且它们也更有可能因未能省略后来被证明是模棱两可的潜在标记而产生假阳性。在这一步骤中，70,927 个 SGB 中的 41,498 个 uSGB 被丢弃，而所有 kSGB 则被保留，因为它们代表了理论上更可靠的序列。

鉴定每个SGB内的核心基因以及核心基因的SGB特异性筛选。

为了识别核心基因，该程序首先根据 SGB 泛基因组定义了核心度百分比阈值（即一个基因在 SGB 中的百分率）。具体来说，我们选择的最大核心度阈值允许检索到至少 800 个核心基因（长度在 450 到 4500 个核苷酸之间）。对于少于 100 个基因组的 SGB，最小核心度阈值为 60%，对于其他 SGB，最小核心度阈值为 50%。使用推断出的核心度阈值为每个 SGB 生成核心基因集。我们平均为每个 SGB 检索到 2,985 个核心基因（中位数，2,687；标准差，1,861）。核心基因少于 200 个的 SGB 将被舍弃，不再考虑（9 个 SGB）。
为了检测 SGB 特异性标记基因，每组核心基因都使用 Bowtie 2（版本 2.3.5.1；-敏感参数）83 与其他 SGB 的基因组进行了比对。出于计算原因，每个 SGB 都选择了质量最高的 100 个基因组子集进行比对。每个核心基因被分割成 150-nt 长的片段，以模拟元基因组读数，然后与 SGB 的代表性基因组子集进行比对。片段的比对命中被视为相应核心基因的命中。没有命中任何其他 SGB 基因组的核心基因（完全唯一标记）或命中率低于 1%（准标记），且命中的 SGB 基因组数量高于或等于其核心度阈值的核心基因被选为标记基因。最重要的是，与以前的 MetaPhlAn 版本相比，这一唯一性程序要严格得多，因为与原始物种分类分配相比，SGB 的一致性得到了提高。

产生少于 100 个标记基因的小部分 SGB（810 个 SGB）,
采用了以下工作流程：
对一个 SGB（目标）提取 marker genes；如果提取出来的 marker genes 少于 100 → 说明可能有问题，比如太像别的菌；于是你检查有哪些 external SGBs “抢占”了你的 marker gene 空间；如果某个 external SGB 匹配到了 目标 SGB 的 200 个以上核心基因，而且占比还不高（<10%）→ 就认为这个 external SGB 是个“噪声”或污染 → 把它排除；然后你再对目标 SGB 重新筛 marker gene，可能就够 100 个了。

如果目标SGB的200多个核心基因与外部SGB（属于同一物种或USGB的KSGB或USGB）匹配，并且外部SGB在目标SGB中的基因组的10％少于10％，则丢弃了外部SGB（这发生在392 KSGB和150 USGB中）。 每次删除外部SGB之前，都会重复此步骤，直到目标SGB产生100个标记基因，或者不再可以评估外部SGB。 在后一种情况下，将外部SGB的去除倒退。

Reconstruct the marker sequences (in fasta format) from the MetaPhlAn BowTie2 database

bowtie2-inspect /data/databases/metaphlan_databases/mpa_vOct22_CHOCOPhlAnSGB_202403  >  mpa_vOct22_CHOCOPhlAnSGB_202403_markers.fasta

根据名称从fasta文件提取序列

samtools faidx  mpa_vOct22_CHOCOPhlAnSGB_202403_markers.fasta
grep ">" mpa_vOct22_CHOCOPhlAnSGB_202403_markers.fasta | head
samtools faidx  mpa_vOct22_CHOCOPhlAnSGB_202403_markers.fasta  "VDB|0046-020C-0-0001|M623-c1000-c0-c0"

Statistical approach for converting marker abundances into clade(分支) abundances

 wnreads = sorted([(float(n)/(np.absolute(r-self.avg_read_length)+1),(np.absolute(r - self.avg_read_length)+1) ,n) for r,n in rat_nreads], key=lambda x:x[0])

metaphlan自带的数据库

  --bowtie2db METAPHLAN_BOWTIE2_DB
                        Folder containing the MetaPhlAn database. You can specify the location by exporting the DEFAULT_DB_FOLDER variable in the shell.[default /opt/miniforge3/envs/gmwi2_env/lib/python3.8/site-packages/metaphlan/metaphlan_databases]

metaphlan_databases
tar -tf mpa_vJan21_CHOCOPhlAnSGB_202103.tar

mpa_vJan21_CHOCOPhlAnSGB_202103.pkl
mpa_vJan21_CHOCOPhlAnSGB_202103_SGB.fna.bz2
mpa_vJan21_CHOCOPhlAnSGB_202103_VINFO.csv
mpa_vJan21_CHOCOPhlAnSGB_202103_VSG.fna.bz2

import pickle
import bz2

db = pickle.load(bz2.open('mpa_vJan21_CHOCOPhlAnSGB_202103.pkl', 'r'))

>>> type(db)
<class 'dict'>
>>> db.keys()
dict_keys(['taxonomy', 'markers', 'merged_taxon'])

>>> list(db['taxonomy'].items())[0]
(
    'k__Bacteria|p__Fusobacteria|c__CFGB15|o__OFGB15|f__FGB15|g__GGB18|s__GGB18_SGB19|t__SGB19', 
    ('2|32066||||||', 1057586)
)

>>> list(db['markers'].items())[0]
(
    'SGB231__BALAOFEP_02157', 
    {
        'clade': 't__SGB231', 
        'ext': [], 
        'len': 3879, 
        'taxon': 'k__Archaea|p__Euryarchaeota|c__Halobacteria|o__Halobacteriales|f__Haloarculaceae|g__Halomicrobium|s__Halomicrobium_katesii|t__SGB231'
    }
)

>>> list(db['merged_taxon'].items())[0]
(
    ('k__Archaea|p__Euryarchaeota|c__Methanomicrobia|o__Methanomicrobiales|f__Methanomicrobiaceae|g__Methanoculleus|s__Methanoculleus_marisnigri|t__SGB177_group', '2157|28890|224756|2191|2194|45989|2198|'), 
    [
        ('k__Archaea|p__Euryarchaeota|c__Methanomicrobia|o__Methanomicrobiales|f__Methanomicrobiaceae|g__Methanoculleus|s__Methanoculleus_sp_MH98A', '2157|28890|224756|2191|2194|45989|1495314', 1), 
        ('k__Archaea|p__Euryarchaeota|c__Methanomicrobia|o__Methanomicrobiales|f__Methanomicrobiaceae|g__Methanoculleus|s__Methanoculleus_chikugoensis', '2157|28890|224756|2191|2194|45989|118126', 1), 
        ('k__Archaea|p__Euryarchaeota|c__Methanomicrobia|o__Methanomicrobiales|f__Methanomicrobiaceae|g__Methanoculleus|s__Methanoculleus_horonobensis', '2157|28890|224756|2191|2194|45989|528314', 1), 
        ('k__Archaea|p__Euryarchaeota|c__Methanomicrobia|o__Methanomicrobiales|f__Methanomicrobiaceae|g__Methanoculleus|s__Methanoculleus_sediminis', '2157|28890|224756|2191|2194|45989|1550566', 1)
    ]
)

创建marker和taxonomy

db['markers']["test_new_marker"] = {
                                   'clade': "t__SGB231",
                                   'ext': {},
                                   'len': 3879,
                                   'taxon': 'k__Archaea|p__Euryarchaeota|c__Halobacteria|o__Halobacteriales|f__Haloarculaceae|g__Halomicrobium|s__Halomicrobium_katesii|t__SGB231'
                                }
db['markers']["test_new_marker"]

db['taxonomy']['k__Archaea|p__Euryarchaeota|c__Halobacteria|o__Halobacteriales|f__Haloarculaceae|g__Halomicrobium|s__Halomicrobium_katesii|t__SGB231'] = ('2|32066||||||', 1057586)

研究一下mpa_vJan21_CHOCOPhlAnSGB_202103.pkl

import bz2
import pickle

with bz2.BZ2File( '/ssd1/wy/workspace/MetaPhlAn/test_data/databases/mpa_vJan21_TOY_CHOCOPhlAnSGB_202103.pkl', 'r' ) as a:
    mpa_pkl = pickle.load( a )

mpa = mpa_pkl
df = pd.DataFrame(columns=['clade','value'])
for clade, value in mpa['taxonomy'].items():
    # print(clade, value)
    df.loc[len(df), df.columns] = clade, value
df.to_csv("taxonomy.tsv", sep='\t',index=False)
df

db['taxonomy']['7-levels taxonomy with clade names of genome1'] = ('7-levels NCBI taxonomy id of genome1', length of genome1)
关于NCBI taxonomy id https://www.ncbi.nlm.nih.gov/taxonomy/?term=60133

df = pd.DataFrame(columns=['marker','clade','ext','len','taxon'])
print(df.columns)
for k,p in mpa['markers'].items():
    # print(p)
    # clade,ext,len,taxon = tuple(p.values())
    df.loc[len(df), df.columns] = k,p['clade'],p['ext'],p['len'],p['taxon']
df.to_csv("marker.tsv", sep='\t',index=False)
df

db['markers'][new_marker_name] = {
'clade': the clade that the marker belongs to,
'ext': {the GCA of the first external genome where the marker appears,
the GCA of the second external genome where the marker appears,
},
'len': length of the marker,
'taxon': the taxon of the marker
}

参数

tree = TaxTree( mpa_pkl, ignore_markers )
tree.set_min_cu_len( pars['min_cu_len'] )
tree.set_stat( pars['stat'], pars['stat_q'], pars['perc_nonzero'], avg_read_length, pars['avoid_disqm'])

• --min_cu_len: minimum total nucleotide length for the markers in a clade for estimating the abundance without considering sub-clade abundances(最小总核苷酸长度的标记在一个进化枝估计丰度而不考虑子进化枝丰度)
• --stat': Statistical approach for converting marker abundances into clade abundances(default tavg_g)(标记丰度转换为分支丰度的统计方法)
  • avg_g : clade global (i.e. normalizing all markers together) average
  • avg_l : average of length-normalized marker counts
  • tavg_g : truncated clade global average at --stat_q quantile
  • tavg_l : truncated average of length-normalized marker counts (at --stat_q)
  • wavg_g : winsorized clade global average (at --stat_q)
  • wavg_l : winsorized average of length-normalized marker counts (at --stat_q)
  • med : median of length-normalized marker counts
• --stat_q: Quantile value for the robust average
• --perc_nonzero: Percentage of markers with a non zero relative abundance for misidentify a species
• --avoid_disqm: Deactivate the procedure of disambiguating the quasi-markers based on the marker abundance pattern found in the sample. It is generally recommended to keep the disambiguation procedure in order to minimize false positives

from metaphlan.metaphlan import TaxTree,TaxClade
tree = TaxTree( mpa_pkl, {} )

for key in tree.all_clades:
    print(key,tree.all_clades[key].nreads)

for marker,reads in sorted(markers2reads.items(), key=lambda pars: pars[0]):
    if marker not in tree.markers2lens:
        continue
    print(marker,reads)
    print(marker,len(reads))

提取核心代码

import bz2
import pickle
from metaphlan.metaphlan import map2bbh
import numpy as np
import pandas as pd

with bz2.BZ2File( '/ssd1/wy/workspace/MetaPhlAn/test_data/databases/mpa_vJan21_TOY_CHOCOPhlAnSGB_202103.pkl', 'r' ) as a:
    mpa_pkl = pickle.load( a )
mpa = mpa_pkl
len(mpa['taxonomy'])
from metaphlan.metaphlan import TaxTree,map2bbh
tree = TaxTree( mpa_pkl, {} )

tree.set_min_cu_len( 2000 )



markers2reads, n_metagenome_reads, avg_read_length = map2bbh("/ssd1/wy/workspace/MetaPhlAn/test_data/01.fasta.bowtie2out.txt", 
                                                           5, "bowtie2out", None, None, None, '1992')
tree.set_stat( 'tavg_g',0.2, 0.33, avg_read_length,False)

map_out = []
for marker,reads in sorted(markers2reads.items(), key=lambda pars: pars[0]):
    if marker not in tree.markers2lens:
        continue
    tax_seq, ids_seq = tree.add_reads( marker, len(reads),
                                add_viruses = False,
                                ignore_eukaryotes = False,
                                ignore_bacteria = False,
                                ignore_archaea = False,
                                ignore_ksgbs = False,
                                ignore_usgbs = False
                                )
    if tax_seq:
        map_out +=["\t".join([r,tax_seq, ids_seq]) for r in sorted(reads)]

tax_lev =None
SGB_ANALYSIS = True
tax_units = "kpcofgst"

clade2abundance_n = dict([(tax_label, clade) for tax_label, clade in tree.all_clades.items()
                    if tax_label.startswith("k__") and not clade.uncl])

clade2abundance, clade2est_nreads, tot_ab, tot_reads = {}, {}, 0.0, 0



def compute_abundance( clade ):
    if clade.abundance is not None: return clade.abundance

    sum_ab = sum([compute_abundance(c) for c in clade.children.values()])


    rat_nreads, removed = [], []
    for marker, n_reads in sorted(clade.markers2nreads.items(),key=lambda x:x[0]):
        misidentified = False

        # if not clade.avoid_disqm:
        #     for ext in clade.markers2exts[marker]:
        #         ext_clade = clade.taxa2clades[ext]
        #         m2nr = ext_clade.markers2nreads
                
        #         tocladetmp = ext_clade
        #         while len(tocladetmp.children) == 1:
        #             tocladetmp = list(tocladetmp.children.values())[0]
        #             m2nr = tocladetmp.markers2nreads

        #         nonzeros = sum([v>0 for v in m2nr.values()])
        #         if len(m2nr):
        #             if float(nonzeros) / len(m2nr) > clade.perc_nonzero:
        #                 misidentified = True
        #                 removed.append( (clade.markers2lens[marker],n_reads) )
        #                 break
        if not misidentified:
            rat_nreads.append( (clade.markers2lens[marker],n_reads) )
        # if not clade.avoid_disqm and len(removed):
        #     n_rat_nreads = float(len(rat_nreads))
        #     n_removed = float(len(removed))
        #     n_tot = n_rat_nreads + n_removed
        #     n_ripr = 10

        #     if len(clade.get_terminals()) < 2:
        #         n_ripr = 0

        #     if "k__Viruses" in clade.get_full_name():
        #         n_ripr = 0

        #     if n_rat_nreads < n_ripr and n_tot > n_rat_nreads:
        #         rat_nreads += removed[:n_ripr-int(n_rat_nreads)]       
    rat_nreads = sorted(rat_nreads, key = lambda x: x[1])
    # rat_v： marker 基因的长度  nreads_v： 每一个marker基因对应的reads个数
    rat_v,nreads_v = zip(*rat_nreads) if rat_nreads else ([],[])
    # rat： marker 基因的长度总和 nrawreads：reads的个数
    rat, nrawreads, loc_ab = float(sum(rat_v)) or -1.0, sum(nreads_v), 0.0
    quant = int(clade.quantile*len(rat_nreads))
    ql,qr,qn = (quant,-quant,quant) if quant else (None,None,0)


    if rat < 0.0:
        pass
    elif clade.stat == 'tavg_g':
        # n: 每一个marker基因对应的reads个数; r: 每一个marker基因的长度; avg_read_length: 平均reads的长度
        wnreads = sorted([(float(n)/(np.absolute(r-clade.avg_read_length)+1),(np.absolute(r - clade.avg_read_length)+1) ,n) for r,n in rat_nreads], key=lambda x:x[0])
        den,num = zip(*[v[1:] for v in wnreads[ql:qr]])
        if sum(num) >0:
            print("111")
        loc_ab = float(sum(num))/float(sum(den)) if any(den) else 0.0

    clade.abundance = loc_ab
    if rat < clade.min_cu_len and clade.children:
        clade.abundance = sum_ab
    elif loc_ab < sum_ab:
        clade.abundance = sum_ab

    if clade.abundance > sum_ab and clade.children: # *1.1??
        clade.uncl_abundance = clade.abundance - sum_ab
    clade.subcl_uncl = not clade.children and clade.name[0] not in tax_units[-2:]

    return clade.abundance


for tax_label, clade in clade2abundance_n.items():
    tot_ab += compute_abundance(clade)

def get_all_abundances(clade ):

    ret = [(clade.name, clade.tax_id, clade.abundance)]
    if clade.uncl_abundance > 0.0:
        lchild = list(clade.children.values())[0].name[:3]
        ret += [(lchild+clade.name[3:]+"_unclassified", "", clade.uncl_abundance)]
    if clade.subcl_uncl and clade.name[0] != tax_units[-2]:
        cind = tax_units.index( clade.name[0] )
        ret += [(   tax_units[cind+1]+clade.name[1:]+"_unclassified","",
                    clade.abundance)]
    for c in clade.children.values():
        ret += get_all_abundances(c)
    return ret


for tax_label, clade in clade2abundance_n.items():
    # print(clade.get_all_abundances())

    # print(get_all_abundances(clade))

    for clade_label, tax_id, abundance in sorted(clade.get_all_abundances(), key=lambda pars:pars[0]):
       clade2abundance[(clade_label, tax_id)] = abundance

ret_d = dict([( tax, float(abundance) / tot_ab if tot_ab else 0.0) for tax, abundance in clade2abundance.items()])
outpred = [(taxstr, taxid,round(relab*100.0,5)) for (taxstr, taxid), relab in ret_d.items() if relab > 0.0]
for clade, taxid, relab in sorted(  outpred, reverse=True,
                                        key=lambda x:x[2]+(100.0*(8-(x[0].count("|"))))):
    add_repr = ''
    print( "\t".join( [clade, 
                                taxid, 
                                str(relab*1.0), 
                                add_repr
                            ] ) + "\n" )

for marker,reads in sorted(markers2reads.items(), key=lambda pars: pars[0]):
    if marker not in tree.markers2lens:
        continue
    print(marker,reads)
    print(marker,len(reads))

for key in tree.all_clades:
    print(key,tree.all_clades[key].markers2nreads)

for k,p in mpa['markers'].items():
    print(k,p)

Marker level analysis

metaphlan -t marker_ab_table /ssd1/wy/workspace/metagenomics/strainphlan/output/bowtie2/SRS013951.fastq.bowtie2 --input_type bowtie2out -o marker_abundance_table.txt

mark count

metaphlan  --bowtie2db /ssd1/wy/workspace/MetaPhlAn/test_data/databases --index mpa_vJan21_TOY_CHOCOPhlAnSGB_202103 /ssd1/wy/workspace/MetaPhlAn/test_data/SRS014476-Supragingival_plaque.fasta --input_type fasta -t marker_counts