在微生物分离培养、分箱中获得的大量的基因组、宏基因组拼接的基因组(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/ 目录中有图片

1. Clustering_scatterplots.pdf
2. Cluster_scoring.pdf
3. Primary_clustering_dendrogram.pdf
4. Secondary_clustering_dendrograms.pdf
5. Winning_genomes.pdf
6. ……

初级聚类图 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_directory

positional arguments:
  work_directory        输出目录Directory where data and output    
                        *** USE THE SAME WORK DIRECTORY FOR ALL DREP OPERATIONS ***

SYSTEM PARAMETERS:
  -p 线程程数，默认6, --processors PROCESSORS
                        threads (default: 6)
  -d, --debug           make extra debugging output (default: False)
  -h, --help            显示帮助 show this help message and exit

FILTERING 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. NOT
                        RECOMMENDED! This is useful for use with
                        bacteriophages or eukaryotes or things where checkM
                        scoring does not work. Will only choose genomes based
                        on length and N50 (default: False)

GENOME COMPARISON PARAMETERS:
  -ms MASH_SKETCH, --MASH_sketch MASH_SKETCH
                        MASH sketch size (default: 1000)
  --S_algorithm {gANI,ANIn,ANImf,fastANI,goANI}
                        聚类比较算法选择，速度和精度不同 Algorithm for secondary clustering comaprisons:
                        fastANI = 超过的ANI计算方法 Kmer-based approach; very fast
                        ANImf   = 默认方法，全基因组比对，过滤比对区域，比较比对区域 (DEFAULT) Align whole genomes with nucmer; filter alignment; compare aligned regions
                        ANIn    = 不过滤比对区域 Align whole genomes with nucmer; compare aligned regions
                        gANI    = 鉴定ORFs并比对 Identify and align ORFs; compare aligned ORFS
                        goANI   = Open source version of gANI; requires nsmimscan
                         (default: ANImf)
  --n_PRESET {normal,tight}
                        Presets to pass to nucmer
                        tight   = 只比对高度保守区 only align highly conserved regions
                        normal  = 默认 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 on
                        all genomes (default: False)
  --SkipSecondary       Skip secondary clustering, just perform MASH
                        clustering (default: False)
  -nc COV_THRESH, --cov_thresh COV_THRESH
                        最小的重叠是10% Minmum level of overlap between genomes when doing
                        secondary 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 CLUSTERALG
                        Algorithm used to cluster genomes (passed to
                        scipy.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_WEIGHT
                        completeness weight (default: 1)
  -conW CONTAMINATION_WEIGHT, --contamination_weight CONTAMINATION_WEIGHT
                        contamination weight (default: 5)
  -strW STRAIN_HETEROGENEITY_WEIGHT, --strain_heterogeneity_weight STRAIN_HETEROGENEITY_WEIGHT
                        strain heterogeneity weight (default: 1)
  -N50W N50_WEIGHT, --N50_weight N50_WEIGHT
                        weight of log(genome N50) (default: 0.5)
  -sizeW SIZE_WEIGHT, --size_weight SIZE_WEIGHT
                        weight of log(genome size) (default: 0)

TAXONOMY:
  --run_tax             generate taxonomy information (Tdb) (default: False)
  --tax_method {percent,max}
                        Method of determining taxonomy
                        percent = The most descriptive taxonimic level with at least (per) hits
                        max     = The centrifuge taxonomic level with the most overall hits (default: percent)
  -per PERCENT, --percent PERCENT
                        minimum percent for percent method (default: 50)
  --cent_index CENT_INDEX
                        path to centrifuge index (for example,
                        /home/mattolm/download/centrifuge/indices/b+h+v
                        (default: None)

WARNINGS:
  --warn_dist WARN_DIST
                        How far from the threshold to throw cluster warnings
                        (default: 0.25)
  --warn_sim WARN_SIM   Similarity threshold for warnings between dereplicated
                        genomes (default: 0.98)
  --warn_aln WARN_ALN   Minimum aligned fraction for warnings between
                        dereplicated genomes (ANIn) (default: 0.25)

I/O PARAMETERS:
  -g [GENOMES [GENOMES ...]], --genomes [GENOMES [GENOMES ...]]
                        genomes to cluster in .fasta format; can pass a number
                        of .fasta files or a single text file listing the
                        locations 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 GENOMEINFO
                        location of .csv file containing quality information
                        on the genomes. Must contain: ["genome"(basename of
                        .fasta file of that genome), "completeness"(0-100
                        value for completeness of the genome),
                        "contamination"(0-100 value of the contamination of
                        the genome)] (default: None)

Example: dRep dereplicate output_dir/ -g /path/to/genomes/*.fast

参考文献

1. 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).