drep：微生物基因组快速去冗余-文章解读+帮助文档+实战教程

在微生物分离培养、分箱中获得的大量的基因组、宏基因组拼接的基因组(MAG)，如何确定到底有多少种非冗余的细菌基因组呢?

来自加州大学伯克利分校Jillian F Banfield组开发的dRep可以帮助我们实现此分析。

摘要

每年测序的微生物基因组的数量正在迅速增加，部分原因是通过宏基因组学研究可以获得数百个质量合格的基因组。已经开发了快速算法来全面比较大型基因组集，但是对于草图质量的基因组来说它们并不准确。在这里，我们开发了dRep程序，通过依次应用快速、不准确的基因组距离估算和较慢准确的平均核苷酸同一性测量，可以减少成对基因组比较的计算时间。与先前开发的算法相比，dRep的速度提高了28倍，具有完美的查全率(recall)和精确度。我们演示了使用dRep从时间序列数据集中恢复基因组。每个宏基因组分别组装，并使用dRep识别基本相同的基因组，并从每个重复组中选择最佳基因组。与使用共装配法回收的基因组相比，这导致回收了更多，更高质量的基因组。

文章主要结果

图1. 与共组装相比，使用dRep进行组装和去重复可产生更多和更高质量的基因组箱。

Assembly and de-replication with dRep results in more and higher-quality genome bins as compared to co-assembly

（a）完整的大肠杆菌基因组以10％（10％，20％，30％等）的增量分为子集10次。使用ANIm，MASH和gANI这三种算法，以成对的方式将子集彼此进行比较（总共进行100次比较）。对于每对子集，由MUMmer确定的两个基因组之间的比对覆盖率显示在x轴上（比对长度/平均基因组长度），而每种算法报告的ANI则显示在y轴上。ANIm和gANI在基因组不完整时是准确的，但是MASH仅在基因组基本完整时才是准确的。

（b）使用先前报告的算法运行时，我们估算了去重复各种大小的基因组所需的时间。gANI表现出急剧的指数攀升，限制了其在更大的基因组集上的使用；MASH和dRep仍然线性增长。

（c）从单样本组装和多样本混合组装中（dRep组装方法）比单独多样本混合组装导获得更多的分箱（完整 > 75％，污染 < 5％）。

（d和e）由dRep生成的基因组相关性数字示例。红色虚线是簇中每个基因组的自对自我比对所产生的最低ANI值。

软件使用

软件主页：https://github.com/MrOlm/drep

安装

conda安装，版本和安装代码详见 http://bioconda.github.io/recipes/drep/README.html

需要在Python3环境下，此次安装版本为2.6.2

conda install drep

如果安装存在问题，可以新建环境。

conda create -n drep
conda activate drep
conda install drep -c bioconda

快速开始

https://drep.readthedocs.io/en/latest/

默认参数计算，指定去冗余功能(dereplicate)，然后是输出目录(outout_directory)和通配符(*)指定多个fasta输入文件

dRep dereplicate outout_directory -g path/to/genomes/*.fasta

结果描述

输出结果的 figures/ 目录中有图片

Clustering_scatterplots.pdf
Cluster_scoring.pdf
Primary_clustering_dendrogram.pdf
Secondary_clustering_dendrograms.pdf
Winning_genomes.pdf
……

初级聚类图 Primary_clustering_dendrogram

基因组间成对mash距离。虚线是初级ANI值，它决定二级簇的起点，越大可显著减少两两比较的数量计算时间，但准确度也容易降低。图中显示形成2个次级簇，分别为蓝色和绿色文字标签。

注：在相同的初级簇的基因组，将会使用更敏感的算法(gANI或ANIm)进行两两比对，以便获得次级簇。不同初级簇的基因组将不再比较。

次级聚类图 Secondary_clustering_dendrograms

每个初级聚类存在多个簇，将在“次级聚类树状图”文件中具有单独的图。在此示例中，只有一个具有> 1个成员的主簇(如果基因组间非冗余，将没有结果)。

该树状图总结了由次级聚类算法（gANI / ANIm）确定的每个初级簇中所有生物基因组之间的成对距离。在最顶部显示了初级簇编号(Primary cluster 1)，在树状图上方显示了有关次级聚类算法参数的信息(gANI,average, 0.1)：次级聚类是使用gANI执行的，最小比对覆盖率为10％，层次聚类方法是平均值。

黑色虚线显示次级聚类的ANI阈值（默认为99％）。该值确定哪些基因组最终位于同初、次级簇中，因此被认为是“相同”的。在上图中，形成了两个次级聚类。每个次级簇的“最佳”基因组都标记为星号。

红线显示的是用于基因组列表中所有基因组“自我”比较的最低ANI。也就是说，当将该主要簇中的每个基因组与其自身进行比较时，红线是您获得的最低ANI。这代表了某种“检测极限”。进行自我与自我比较时，gANI始终会产生100％的ANI，但ANIm不会（如下图所示）。

簇得分 Cluster_scoring

每个次级簇在“簇评分”图中将有其自己的页面。包括完整度(completeness)、污染率(contamination)、株水平杂合度(stran heterogenetity)、总长(length)、N50和得分(score)。在此示例中，有三个次级簇-其中2个来自初级簇1，其中1个是初级簇2的唯一成员。

这些数字显示了每个基因组的得分，以及可以用来确定分数的所有指标。这有助于用户可视化所有基因组的质量，并确保它们与被选为“最佳”的基因组一致。“最佳”基因组带有星号，并且始终是得分最高的基因组。

从每个次级簇中选择一个基因组，将其包括在去重复的基因组中。因此，在以上示例中，我们将在去重复的基因组集合中拥有3个基因组。这是因为该算法确定簇1_1中的所有基因组都是“相同”的，并选择GCA_900083945作为“最佳”。

其它图 Other figures

Clustering scatterplots (Clustering_scatterplots.pdf)：

基因组比对的统计
Winning genomes(Winning_genomes.pdf): 每个重复集的“最佳”基因组，以及几个快速的整体统计数据。

实战数据

我们这里以Metawrap分批次样本/组分箱的结果为例(同一个分箱的结果基本没有重复)，本质上输入文件为多个基因组或基因组草图(bins)。测试数据和结果可以在中下载：https://github.com/YongxinLiu/Note/tree/master/Meta/dRep

# 我们使用了4个样本分别分箱获得的6个草图为例，位于bin目录下的*.fq.gz，先解压数据
gunzip bin/*.fa.gz# 默认参数，输出目录为当前，输入为bin中的fa文件
dRep dereplicate ./ -g bin/*.fa

用时10分，结果计算过程显示如下：中文注释见标题下

***************************************************..:: dRep dereplicate Step 1. Filter ::..
***************************************************
第一步过滤
Will filter the genome list
Calculating genome info of genomes
按长度过滤<50k的基因组，全部通过
100.00% of genomes passed length filtering
Running prodigal
Running checkM
按完整度>75，污染<25，仅有6*1/3=2个通过
33.33% of genomes passed checkM filtering
***************************************************..:: dRep dereplicate Step 2. Cluster ::..
***************************************************
第二步聚类
Clustering Step 1. Parse Arguments
Clustering Step 2. Perform MASH (primary) clustering
2a. Run pair-wise MASH clustering
2b. Cluster pair-wise MASH clustering
1 primary clusters made
Step 3. Perform secondary clustering
Running 4 ANImf comparisons- should take ~ 0.3 min
Step 4. Return output
***************************************************..:: dRep dereplicate Step 3. Choose ::..
***************************************************
第三步：选择最优基因组
Loading work directory
........:: dRep dereplicate finished ::..
完成后输出结果如下：
非冗余结果 Dereplicated genomes................. ./dereplicated_genomes/
Dereplicated genomes information..... ./data_tables/Widb.csv
图片是重点 Figures.............................. ./figures/
Warnings............................. ./log/warnings.txt

统计输出结果，去冗余后仅剩2个，4个被质量过滤，1个冗余

ls temp/dRep/dereplicated_genomes/bin*.fa|wc -l

图. 初级聚类结果，两个菌间ANI为97%

我们需要根据实际情况调整一些参数。不让基因组被过滤掉。

常用参数

比对参考NBT 2020最新人类基因集的参数。按95% ANI聚类(默认0.99)，最小重叠为30%(默认0.1)，并采用多线程加速

MAG在完整度较低，污染率较高，默认为75和25，可修改中等质量阈值>50%完整度，且<10%污染率

dRep dereplicate ./ -g bin/*.fa -sa 0.95 -nc 0.30 -p 24 -comp 50 -con 10

统计输出非冗余结果

ls ./dereplicated_genomes/bin*.fa|wc -l

查看日志

less -S log/logger.log

去冗余后3个(3/6=50%)，有1半菌ANI<95%被去重复。

图Primary_clustering_dendrogram.pdf. 6个分箱的初级聚类结果，虚拟显示生成了3个初级簇

图Secondary_clustering_dendrograms.pdf. 初级簇1中3个成员，小于95%ANI，挑选了最下方基因组为最佳

图Winning_genomes.pdf. 按完整度、污染率对簇进行评估，并显示选择最优的依据。

dRep dereplicate参数详解

dRep dereplicate -h

帮助如下：重点参数有中文翻译和注释

usage: dRep dereplicate [-p 线程数] [-d] [-h] [-l LENGTH][-comp 完整度] [-con 污染率][--ignoreGenomeQuality] [-ms MASH_SKETCH][--S_algorithm {gANI,ANIn,ANImf,fastANI,goANI}][--n_PRESET {normal,tight}] [-pa P_ANI] [-sa S_ANI][--SkipMash] [--SkipSecondary] [-nc COV_THRESH][-cm {total,larger}] [--clusterAlg CLUSTERALG][-comW COMPLETENESS_WEIGHT][-conW CONTAMINATION_WEIGHT][-strW STRAIN_HETEROGENEITY_WEIGHT] [-N50W N50_WEIGHT][-sizeW SIZE_WEIGHT] [--run_tax][--tax_method {percent,max}] [-per PERCENT][--cent_index CENT_INDEX] [--warn_dist WARN_DIST][--warn_sim WARN_SIM] [--warn_aln WARN_ALN][-g [GENOMES [GENOMES ...]]][--checkM_method {lineage_wf,taxonomy_wf}][--genomeInfo GENOMEINFO]work_directorypositional arguments:work_directory        输出目录Directory where data and output*** USE THE SAME WORK DIRECTORY FOR ALL DREP OPERATIONS ***SYSTEM PARAMETERS:-p 线程程数，默认6, --processors PROCESSORSthreads (default: 6)-d, --debug           make extra debugging output (default: False)-h, --help            显示帮助 show this help message and exitFILTERING OPTIONS:-l LENGTH, --length LENGTH最小基因组长度，默认50k，Minimum genome length (default: 50000)-comp COMPLETENESS, --completeness COMPLETENESS最小的基因组完整度，默认75 Minumum genome completeness (default: 75)-con CONTAMINATION, --contamination CONTAMINATION最大的基因组污染率，默认25 Maximum genome contamination (default: 25)--ignoreGenomeQuality不运行checkM进行质量过滤 Don't run checkM or do any quality filtering. NOTRECOMMENDED! This is useful for use withbacteriophages or eukaryotes or things where checkMscoring does not work. Will only choose genomes basedon length and N50 (default: False)GENOME COMPARISON PARAMETERS:-ms MASH_SKETCH, --MASH_sketch MASH_SKETCHMASH sketch size (default: 1000)--S_algorithm {gANI,ANIn,ANImf,fastANI,goANI}聚类比较算法选择，速度和精度不同 Algorithm for secondary clustering comaprisons:fastANI = 超过的ANI计算方法 Kmer-based approach; very fastANImf   = 默认方法，全基因组比对，过滤比对区域，比较比对区域 (DEFAULT) Align whole genomes with nucmer; filter alignment; compare aligned regionsANIn    = 不过滤比对区域 Align whole genomes with nucmer; compare aligned regionsgANI    = 鉴定ORFs并比对 Identify and align ORFs; compare aligned ORFSgoANI   = Open source version of gANI; requires nsmimscan(default: ANImf)--n_PRESET {normal,tight}Presets to pass to nucmertight   = 只比对高度保守区 only align highly conserved regionsnormal  = 默认 default ANIn parameters (default: normal)CLUSTERING PARAMETERS:-pa P_ANI, --P_ani P_ANI默认聚类参数90% ANI threshold to form primary (MASH) clusters(default: 0.9)-sa S_ANI, --S_ani S_ANI二级聚类为99% ANI threshold to form secondary clusters (default:0.99)--SkipMash            Skip MASH clustering, just do secondary clustering onall genomes (default: False)--SkipSecondary       Skip secondary clustering, just perform MASHclustering (default: False)-nc COV_THRESH, --cov_thresh COV_THRESH最小的重叠是10% Minmum level of overlap between genomes when doingsecondary comparisons (default: 0.1)-cm {total,larger}, --coverage_method {total,larger}Method to calculate coverage of an alignment(for ANIn/ANImf only; gANI and fastANI can only do larger method)total   = 2*(aligned length) / (sum of total genome lengths)larger  = max((aligned length / genome 1), (aligned_length / genome2))(default: larger)--clusterAlg CLUSTERALGAlgorithm used to cluster genomes (passed toscipy.cluster.hierarchy.linkage (default: average)SCORING CRITERIA
Based off of the formula:
A*Completeness - B*Contamination + C*(Contamination * (strain_heterogeneity/100)) + D*log(N50) + E*log(size)A = completeness_weight; B = contamination_weight; C = strain_heterogeneity_weight; D = N50_weight; E = size_weight:-comW COMPLETENESS_WEIGHT, --completeness_weight COMPLETENESS_WEIGHTcompleteness weight (default: 1)-conW CONTAMINATION_WEIGHT, --contamination_weight CONTAMINATION_WEIGHTcontamination weight (default: 5)-strW STRAIN_HETEROGENEITY_WEIGHT, --strain_heterogeneity_weight STRAIN_HETEROGENEITY_WEIGHTstrain heterogeneity weight (default: 1)-N50W N50_WEIGHT, --N50_weight N50_WEIGHTweight of log(genome N50) (default: 0.5)-sizeW SIZE_WEIGHT, --size_weight SIZE_WEIGHTweight of log(genome size) (default: 0)TAXONOMY:--run_tax             generate taxonomy information (Tdb) (default: False)--tax_method {percent,max}Method of determining taxonomypercent = The most descriptive taxonimic level with at least (per) hitsmax     = The centrifuge taxonomic level with the most overall hits (default: percent)-per PERCENT, --percent PERCENTminimum percent for percent method (default: 50)--cent_index CENT_INDEXpath to centrifuge index (for example,/home/mattolm/download/centrifuge/indices/b+h+v(default: None)WARNINGS:--warn_dist WARN_DISTHow far from the threshold to throw cluster warnings(default: 0.25)--warn_sim WARN_SIM   Similarity threshold for warnings between dereplicatedgenomes (default: 0.98)--warn_aln WARN_ALN   Minimum aligned fraction for warnings betweendereplicated genomes (ANIn) (default: 0.25)I/O PARAMETERS:-g [GENOMES [GENOMES ...]], --genomes [GENOMES [GENOMES ...]]genomes to cluster in .fasta format; can pass a numberof .fasta files or a single text file listing thelocations of all .fasta files (default: None)--checkM_method {lineage_wf,taxonomy_wf}Either lineage_wf (more accurate) or taxonomy_wf(faster) (default: lineage_wf)--genomeInfo GENOMEINFOlocation of .csv file containing quality informationon the genomes. Must contain: ["genome"(basename of.fasta file of that genome), "completeness"(0-100value for completeness of the genome),"contamination"(0-100 value of the contamination ofthe genome)] (default: None)Example: dRep dereplicate output_dir/ -g /path/to/genomes/*.fast

参考文献

Olm, M. R., Brown, C. T., Brooks, B. & Banfield, J. F. dRep: a tool for fast and accurate genomic comparisons that enables improved genome recovery from metagenomes through de-replication. The ISME Journal 11, 2864-2868, doi:10.1038/ismej.2017.126 (2017).

猜你喜欢

10000+：菌群分析宝宝与猫狗梅毒狂想曲提DNA发Nature Cell专刊肠道指挥大脑

系列教程：微生物组入门 Biostar 微生物组宏基因组

专业技能：学术图表高分文章生信宝典不可或缺的人

一文读懂：宏基因组寄生虫益处进化树

必备技能：提问搜索 Endnote

文献阅读热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读分析流程统计绘图

16S功能预测 PICRUSt FAPROTAX Bugbase Tax4Fun

在线工具：16S预测培养基生信绘图

科研经验：云笔记云协作公众号

编程模板: Shell R Perl

生物科普: 肠道细菌人体上的生命生命大跃进细胞暗战人体奥秘

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍未解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读