机构类网站有哪些,html5新手做的网页,免费私人网站,logo标志文章目录 Install依赖关系 常用命令常见问题pplacer线程超过30报错当比较基因组很多#xff08;4096#xff09;有了Bdv.csv文件后无需输入基因组list 超多基因组为什么需要界定种#xff1f;dRep重要概念次级ANI的选择Minimum alignment coverage3. 选择有代表性的基因… 文章目录 Install依赖关系 常用命令常见问题pplacer线程超过30报错当比较基因组很多4096有了Bdv.csv文件后无需输入基因组list 超多基因组为什么需要界定种dRep重要概念次级ANI的选择Minimum alignment coverage3. 选择有代表性的基因组4. 使用贪婪的算法5. 基因组完整性的重要性7. 基因组比较算法概述8. 对非细菌基因组进行比较和去重复 dRep结果处理1.获得OTU members2.比较两种生境基因组ANI距离分布密度图 关于代表基因组选择使用测试去重复 参考 Install
(base) [yutio02 01_Morchella]$ mamba create -n drep -c bioconda drep依赖关系
dRep 需要其他程序才能运行。并非所有操作都需要所有依赖程序要检查系统中安装了哪些依赖项dRep 可以访问这些依赖项请运行
(drep) [yutio02 01_Morchella]$ dRep check_dependencies
mash.................................... all good (location /datanode02/yut/Software/Miniconda3/envs/drep/bin/mash)
nucmer.................................. all good (location /datanode02/yut/Software/Miniconda3/envs/drep/bin/nucmer)
checkm.................................. !!! ERROR !!! (location None)
ANIcalculator........................... !!! ERROR !!! (location None)
prodigal................................ all good (location /datanode02/yut/Software/Miniconda3/envs/drep/bin/prodigal)
centrifuge.............................. !!! ERROR !!! (location None)
nsimscan................................ !!! ERROR !!! (location None)
fastANI................................. !!! ERROR !!! (location /datanode02/yut/Software/Miniconda3/envs/drep/bin/fastANI)安装其他依赖
# install
$ mamba install -c bioconda checkm-genome# download database
$ wget https://data.ace.uq.edu.au/public/CheckM_databases/checkm_data_2015_01_16.tar.gz
$ [/datanode02/yut/Database/checkm_data_2015_01_16] tar -zxvf checkm_data_2015_01_16.tar.gz
(drep) [yutnode02 checkm_data_2015_01_16]$ checkm data setRoot /datanode02/yut/Database/checkm_data_2015_01_16
[2023-11-10 20:06:35] INFO: CheckM v1.2.2
[2023-11-10 20:06:35] INFO: checkm data setRoot /datanode02/yut/Database/checkm_data_2015_01_16
[2023-11-10 20:06:35] INFO: CheckM data: /datanode02/yut/.checkm
[2023-11-10 20:06:35] INFO: [CheckM - data] Check for database updates. [setRoot]Path [/datanode02/yut/Database/checkm_data_2015_01_16] exists and you have permission to write to this folder.
(re) creating manifest file (please be patient).Path [/datanode02/yut/Database/checkm_data_2015_01_16] exists and you have permission to write to this folder.
(re) creating manifest file (please be patient).常用命令
输入文件可以是.fna或者fna.gz或者混合的
(gtdbtk) [yutaoglobin test]$ ls *fna*
fna.fa GCA_000685155.1_ANME2D_V10_genomic.fna.gz test.fna(gtdbtk) [yutaoglobin test]$ time dRep dereplicate drep_out --ignoreGenomeQuality -p 20 -g *fna*(gtdbtk) [yutaoglobin test]$ cat drep_out/data_tables/genomeInformation.csv
location,length,N50,genome,centrality
/home/yutao/test/fna.fa,3203386,545726,fna.fa,0.0
/home/yutao/test/GCA_000685155.1_ANME2D_V10_genomic.fna.gz,3203386,545726,GCA_000685155.1_ANME2D_V10_genomic.fna.gz,0.0
/home/yutao/test/test.fna,199347,46542,test.fna,0.0(drep) [uh3241TGG_High_Quality_Genomes]$ nohup time dRep dereplicate 3241Genomes_dRep_0.99 -comp 0 -con 100 -p 28 -sa 0.99 -g 3241TGG_genomes/*.fa drep_0.99.log(drep) [uhGEM_env_sets]$ nohup dRep dereplicate GEM_ENV_dRep0.95 --ignoreGenomeQuality -p 28 -sa 0.95 -g *.fa gem0.95.log常见问题
pplacer线程超过30报错
dRep v3.2.2 超过30线程pplacer进程一直在但是没有计算结果。
当比较基因组很多4096
当比较基因组很多时某些系统会限制参数打开个数需要将基因组的路径写到一个文件跟在-g
(gtdbtk) [yutaomyosin release202]$ find $PWD/tmp -iname *fna* -type f 48862genome_path.tsv
(gtdbtk) [yutaomyosin release202]$ nohup time dRep dereplicate 968OTU_vs_RS202_47894OTUs_dRep0.95 --ignoreGenomeQuality -p 80 -sa 0.95 -g 48862genome_path.tsv drep.log注意一定要将路径写对有一个不对则会导致checkM错误
有了Bdv.csv文件后无需输入基因组list
如果有了Bdv.csv文件则w无需输入基因组list直接去掉-g选项及参数
超多基因组
使用--ignoreGenomeQuality -p 150 --S_algorithm fastANI --multiround_primary_clustering --greedy_secondary_clustering快速比较5000基因组
(gtdbtk) [yutaoglobin GCS_OTUs_vs_OtherDB]$ nohup time dRep dereplicate GCSvsGEM_FastANI_dRep95 --ignoreGenomeQuality -p 100 --S_algorithm fastANI --multiround_primary_clustering --greedy_secondary_clustering -sa 0.95 -g 47516_GEM.path 47516_GEM_fastANI.drep.log
# 近48k基因组耗时69h# 可通过以下命令查看进度
$ grep -c running cluster GCSvsGEM_FastANI_dRep95/log/logger.log
39141 #已经聚类的数量
$ grep primary clusters made GCSvsGEM_FastANI_dRep95/log/logger.log
10-25 05:54 INFO 39141 primary clusters made # 总的二级聚类数量nohup time dRep dereplicate GCSvsGEM_dRep95 --ignoreGenomeQuality --multiround_primary_clustering -p 80 -sa 0.95 -g 47516_GEM.path 47516_GEM_drep95.log
# 不做checkM超过30核
# 耗时超过6天未出结果256threads 500G因为secondary cluster是成对的并且不是使用FastaANI非常慢线程大于30的时候--ignoreGenomeQuality要放到最前面否则报错
为什么需要界定种
在许多情况下确定微生物之间的关系是研究问题的中心。 居住在建筑物表面的微生物是否与居住在其租户中的微生物相同 医院病房中的微生物是否与新生婴儿中的微生物相同 生活在木制物表面的大肠杆菌与生活在塑料的大肠杆菌一样吗常常通过平均核酸相似性Average Nucleotide Identity, ANI)来衡量。 基本思想是比对两个基因组并计算比对中错配的数量。 例如ANI为99的基因组每100个碱基之间有1个错配而ANI为95的基因组每100个碱基之间有5个错配依此类推。 有许多计算平均核苷酸同一性的方法主要区别在于用于比对基因组的算法差异。2–4 通过计算许多系统中基因组之间的ANI已记录了一些松散的和一般的ANI断点 96% ANI Same 16S cluster (using standard 97% clustering)596% ANI Same bacterial species498% ANI Same E. coli clade698.8% ANI Same Prochlorococcus clade799.9% ANI Same K. pneumoniae outbreak strain8
在哪个ANI阈值处将基因组称为“相同”更为合理这取决于研究问题。 如果问题是弗拉格斯塔夫办公室中的微生物与圣地亚哥办公室中的微生物是否相同那么属于同一个species微生物即可因此ANI为95或16S测序 将足以解决这个问题确实如此 Chase20169。 如果问题是两个不同身体部位的微生物是否来自同一来源那么95ANI太低了。 仅仅因为大肠杆菌位于两个身体部位并不意味着它们来自同一地点 一种菌株可能来自土壤另一种菌株来自隔壁的邻居。 绝对需要95以上的ANI才能显示两种菌株都来自同一来源但是到底需要多高的ANI这其实是另一个问题在最近的出版物中使用99.9的ANI来解决这个问题 Olm201610。
为特定问题选择ANI阈值时可视化基因组之间的关系通常会很有帮助。 dRep是最近在bioRxiv11上发布的python程序旨在实现这一目的。 例如 上图显示了匹兹堡同一重症监护室中居住在不同婴儿中的链霉菌Streptomyces 菌株之间的ANI。 从图中可以看到conN3_174_037G1_concoct_13和conN1_023_029G1_concoct_18之间的ANI约为99.25conN3_174_023G1_concoct_19和conN3_174_021G1_concoct_4之间的ANI约为100依此类推。 该图还清楚地表明不同的ANI阈值将得出关于哪些婴儿具有“相同”菌株的不同结论。 例如如果其ANI 98.5如黑色虚线所示则称基因组为“相同”same将得出结论即所有婴儿都只有一个链霉菌菌株。 如果其ANI 99.5如红色虚线所示则称基因组为“相同”将得出结论有4个不同的链霉菌菌株上面2个可视为同一种中间三个可视为同一种最下面2个是2个种其中两个conN3_174_037G1_concoct_13和conN1_023_029G1_concoct_18同在一个婴儿中。 在此示例中将ANI阈值更改单个百分点会完全更改从数据得出的结论从而突出显示了谨慎选择阈值的重要性。
dRep的大致流程 1.Assemble each sample separately using your favorite assembler. You can also perform a co-assembly to catch low-abundance microbes2.Bin each assembly (and co-assembly) separately using your favorite binner3.Pull the bins from all assemblies together and run dRep on them4.Perform downstream analysis on the de-replicated genome list
dRep重要概念
在使用dRep进行基因组去重复和基因组比较的时候需要理解以下几个重要概念
1.选择合适的次级ANI阈值这个阈值规定了基因组需要有多相似才能被认为是相同的。本节将介绍如何选择适当的S_ANI阈值2.最小比对覆盖率。在计算基因组之间的ANI时我们还确定了基因组中进行比较的覆盖度fraction。该值可用于建立基因组之间所需的最小重叠水平(cov_thresh)但在使用此参数时需要注意很多问题3.选择具有代表性的基因组。在去重复过程中第一步是识别相似基因组的类别第二步是为每个类别选择一个代表性基因组(RG)4.使用贪婪算法。贪心算法是走捷径来获得解的算法。dRep可以在几个阶段使用贪婪算法来加快运行速度;本节将介绍如何使用它们以及需要注意什么5.基因组完整性的重要性。由于它对Mash基因组的依赖dRep必须在一定程度上完成对它们的处理。本节解释原因6.奇怪的分层聚类。 在比较所有基因组之后dRep使用分层聚类将成对的ANI值转换为基因组簇。本节介绍了这件事以及它可能出错的地方7.基因组比较算法概述。dRep有很多不同的基因组比较算法可供选择。本节简要介绍它们的工作原理和它们之间的区别8.比较和去重复非细菌基因组。在比较噬菌体、质粒和真核生物时一些关于使用适当参数的想法
次级ANI的选择
1.对于两个 相同 的基因组之间共享的平均核苷酸身份ANI没有一个定义。这是一个用户必须根据自己的具体应用情况自行决定的。ANI是由二级聚类算法决定的最小二级ANI是被认为是 相同 的基因组之间的最小ANI最小对齐部分是相信报告的ANI值的最小基因组重叠量。见7. 基因组比较算法概述中关于二级聚类算法的描述和2. 关于调整最小对齐分数的概念见7.最小对齐覆盖率。基因组去重复的一个用例是把一大群基因组拿出来挑出一套物种级的代表基因组SRG。最近Almeida等人使用dRep完成了这项工作其目的是为生活在人类肠道中的微生物产生一套全面的高质量参考基因组。对于划分不同物种的细菌的确切阈值存在争议但大多数研究都认为95%的ANI是物种水平去重的适当阈值。关于如何选择这个阈值的背景以及关于细菌物种在自然界中是作为一个连续体还是作为离散单位存在的一些争论见Olm等人2020。去重复的另一个用途是产生一组足够不同的基因组以避免错误的映射。如果我们将元基因组读数通常是~150bp映射到两个99.9%相同的基因组上大多数读数将同样映射到两个基因组上而你将无法分辨哪个读数真正 属于 一个基因组或另一个。如果去重复的目标是在映射短读数时产生一组不同的基因组那么98%的ANI是一个合适的阈值。关于98%的数值是如何得出的请看下面的页面。 dRep中的默认值是99%这是对二级ANI阈值的精确程度的限制。作为一个经验法则阈值高达99.9%可能是安全的但高于这个值就超出了检测的极限。将基因组与自身进行比较通常会产生~99.99%的ANI值因为算法并不完美会被基因组重复区域和类似的东西甩开。参见InStrain程序的出版物了解有关测试dRep检测极限的细节以及如果需要的话如何进行更详细的染色剂水平比较。
Minimum alignment coverage
最小对齐分数是基因组必须重叠的最小数量以使报告的ANI值 “计数”。这个值是作为ANIm/gANI算法的一部分报告的。 想象一下这样的情景两个远缘的基因组共享一个相同的转座子。当基因组比较算法运行时转座子可能是基因组中唯一对齐的部分并且对齐后的ANI值为100%。这将导致报告的ANI为100%而报告的对齐率为~0.1%。最小对齐分数是为了处理上述情况–任何少于最小对齐分数的基因组对齐都会使ANI变为0。 dRep处理最小对齐分数的方式可以说明上述情况但总体上是非常笨拙的。当一个比较低于最小对齐分数时dRep实际上只是将ANI从算法报告的任何内容改为 “0”。这可能会给接下来的分层聚类带来问题见7.基因组比较算法的概述。我试着思考更好的处理方法但这是一个难以解释的问题。你如何处理一组在20%的基因组上共享100% ANI的基因组然而在使用dRep多年后我得出这样的结论在实际情况下这很少是一个问题。在大多数情况下对齐的部分要么非常高50%要么非常低10%所以10%的默认值在大多数情况下是有效的。参见Olm等人2020年的图1描述了典型的物种内基因组共享95% ANI的比对覆盖率值注意95% ANI的基因组几乎不会有低比对覆盖率。 有人建议物种级别的分类定义应该采用最低60%的对齐比例Varghese 2015。然而当使用不完整的基因组时这可能太严格了正如基因组去重复时经常出现的情况。
3. 选择有代表性的基因组
dRep使用一个基于分数的系统来挑选代表基因组。聚类中的每个基因组都被赋予一个分数分数最高的基因组被选为代表。这个分数是基于以下公式。
A∗Completeness−B∗ContaminationC∗(Contamination∗(strainheterogeneity/100))D∗log(N50)E∗log(size)F∗(centrality−Sani)其中A-F是命令行参数默认值分别为1、5、1、0.5、0和1。调整A-F可以让你决定在选择代表性基因组时对特定特征的权重。例如如果你真的关心有低污染和高N50你可以增加B和D。
完整性、污染性和菌株异质性由用户提供或用checkM计算。N50是衡量构成基因组的片段有多大。大小是基因组的总长度。中心度是衡量一个基因组与它的群组中所有其他基因组的相似程度。这个指标有助于挑选与所有其他基因组相似的基因组并避免挑选相对离群的基因组。
一些出版物在挑选有代表性的基因组时在其评分中加入了其他指标如基因组是否来自分离物。例子见《人类肠道微生物组204,938个参考基因组的统一目录》和《细菌和古细菌的完整领域-物种分类法》。
4. 使用贪婪的算法
通俗地说贪婪算法是那些走捷径的算法运行速度更快得出的解决方案可能不是最优的但却 “足够接近”。贪婪算法越好最优解和贪婪解之间的差异就越小。由于成对比较迅速扩展到需要数十年计算的水平dRep使用一些贪婪算法来加速。
dRep使用的一种贪婪算法是初级聚类。执行这一步骤可以极大地减少必须进行的基因组比较的数量减少运行时间。这样做的代价是如果基因组最终出现在不同的主聚类中它们将永远不会被比较因此也永远不会出现在相同的最终聚类中。这就是为什么下面的章节基因组完整性的重要性很重要。 在2021年dRep v3引入了几个额外的贪婪算法描述如下。这些都是相对较新的功能所以如果你发现问题或有建议请不要犹豫地联系我们。
-multiround_primary_clustering在一系列的组中进行初级聚类然后在最后用单链路聚类进行合并。这极大地减少了RAM的使用提高了速度但代价是精度略有损失而且不能绘制primary_clustering_dendrograms。在对5000多个基因组进行聚类时尤其有帮助。
-greedy_secondary_clustering在进行二次聚类时使用启发式方法来避免成对比较。其工作方式是随机选择一个基因组来代表一个聚类。然后将下一个基因组与该基因组进行比较。如果它与该基因组的ANI阈值低它将被放入该群组。如果不是它将被放入一个新的集群并成为新集群的代表基因组。然后第3个基因组将与所有集群的代表进行比较以此类推。这基本上导致了单链路聚类而不需要进行成对的比较。不幸的是由于FastANI需要不断地重新绘制基因组图谱这并没有像你期望的那样提高速度。这个选项目前只适用于FastANI的S_algorithm。 -run_tertiary_clustering不是一种贪婪的算法而是一种处理贪婪算法所带来的潜在不一致的方法。一旦聚类完成并选择了有代表性的基因组这个选项就会对最终的基因组集再运行一轮聚类。这在进行贪婪聚类和/或处理类似的基因组最终出现在不同的主聚类中的情况时特别有用。这基本上是一个检查以确保所有的基因组都是基于给定的参数而彼此不同的。
5. 基因组完整性的重要性
这个决定要比前者复杂得多。从本质上讲在计算效率和最小基因组完整性之间存在着一种权衡。 图A五个基因组被细分为10%-100%的部分并对同一基因组的部分进行比较。X轴是允许的最小基因组完整度。这个值越宽松对齐的分数范围就越大。 如上图A所示基因组完整性的极限越低两个基因组可能对齐的部分就越低。这是有道理的–如果你随机抽取一个基因组的20%然后再做同样的事情当你比较这两个随机的20%的子集时你不会期望它们中有多少是对齐的。当你考虑到它对Mash的影响时这个 对齐的部分 就真正成为一个问题。 图B: 一个相同的大肠杆菌基因组被子集到10%-100%的分量并对分量进行比较。当较低数量的基因组对齐时由于不完整Mash ANI会受到严重的影响
如上图B所示对于相同的基因组对齐的部分越低报告的Mash ANI越低。
请记住–基因组首先用Mash分为一级簇然后每个一级簇被分为 相同 基因组的二级簇。因此符合 相同 定义的基因组必须最终出现在同一个主簇中否则程序永远不会意识到它们是相同的。由于更多的不完整的基因组具有较低的Mash值即使这些基因组是真正相同的见图B你允许进入你的基因组列表的不完整的基因组越多你必须减少主簇的阈值。
有一个较低的主集群阈值这将导致更大的主集群这将导致更多需要的二级比较。这将导致更长的运行时间。
还听我的吗
例如假设我将最低基因组完整性设置为50%。如果我取一个大肠杆菌基因组将其50%的基因组子集2次然后将这2个子集基因组放在一起比较Mash将报告96%的ANI。因此主聚类阈值必须至少为96%否则这两个基因组可能最终处于不同的主聚类中因此永远不会有二级算法在它们之间运行因此也就不会被去重。
然而你不想把主集群的阈值设置得过低因为这将导致更多的基因组被包括在每个主集群中从而导致更多的二级比较这很慢因此运行时间也会更高。
把这些放在一起我们可以得到一个数字即相同的基因组被子集到不同部分的最低报告ANI。这个数字只考虑到了5个不同的基因组但给出了一个大致的限制概念。 最后要考虑的是当真正运行dRep时用户实际上并不知道他们的基因组有多不完整。他们必须依靠单拷贝基因清单这样的指标来告诉他们。这就是dRep目前不支持噬菌体和质粒的原因–没有办法知道它们有多完整因此也没有办法过滤掉那些太不完整的分类。不过总的来说checkM在访问基因组的完整性方面是非常好的。 挑选基因组完整度阈值的一些一般准则。
不建议低于50%的完整度。所得到的基因组无论如何都是非常糟糕的甚至二级算法在这一点上也会崩溃。 降低二级ANI应该导致MASH ANI的全面降低。这是因为你想让Mash将非相似和不完整的基因组分组。
7. 基因组比较算法概述
一级聚类总是用Mash进行这是一种极快但有些不准确的算法。
有几种支持的二级聚类算法。这些算法计算出基因组之间准确的平均核苷酸身份ANI用于将基因组聚成二级聚类。目前从第3版开始支持以下算法。
ANInRichter 2009。这是用nucmer对整个基因组进行比对并对比对的区域进行比较。 ANImf (DEFAULT)。这与ANIn相同但对排列进行过滤使基因组1的每个区域只与基因组2的一个区域对齐。这需要稍多的时间但在有重复区域的基因组上更准确。 gANIVarghese 2015。这是对Prodigal调用的基因ORFs进行对齐而不是对整个基因组进行对齐。这个算法比基于ANIm的算法快一点但只对齐编码区。 goANI。这是我自己对gANI的开源实现gANI是不开源的被问及时作者也不会分享源代码。我写了这个算法以便我可以为这项研究计算对齐的基因之间的dN/dS你也可以使用dnds_from_drep.py。需要程序NSimScan。 FastANIJain 2018。一个真正快速的基于Mash的算法也可以处理不完整的基因组。似乎与基于对齐的算法一样准确。当你关心运行时间的时候可能应该是默认的算法。 这些算法都不是完美的特别是在容易重复的基因组中。基因组中非同源的区域可以相互对齐人为地降低ANI。事实上当一个基因组与它自己进行比较时由于这个原因这些算法经常报告的数值100%。
8. 对非细菌基因组进行比较和去重复
dRep是针对细菌解重复的使用情况而开发的在非细菌实体上运行时有一些事情需要注意。
需要注意的一个主要问题是初级聚类。正如在5.基因组完整性的重要性中所描述的基因组需要50%的完整性才能使初级聚类起作用。如果你正在比较那些你无法评估完整性的实体或者你想比较那些只共享有限数量基因的基因组如噬菌体或质粒这就是一个问题。最简单的处理方法是用参数-SkipMash完全避免一级聚类或者用-pa降低一级聚类的阈值。
还要考虑对准覆盖率2.最小对准覆盖率对分层聚类的影响6.分层聚类的怪异性。如果你的工作对象是那些特别具有镶嵌性的实体如噬菌体这可能是一个比细菌更大的问题。
基因组的过滤和评分也是一个主要因素。如果你的基因组不能被checkM评估你可以用标志-ignoreGenomeQuality关闭质量过滤以及在挑选基因组时使用完整性和污染性。
在为我自己的研究Olm 2019和Olm 2020考虑这些选项时我采用了以下dRep命令来对噬菌体和质粒进行聚类。请把这看作是一个人为他们的具体使用情况处理这些参数的尝试不要害怕在你认为合适的时候做额外的调整。
dRep dereplicate --S_algorithm ANImf -nc .5 -l 10000 -N50W 0 -sizeW 1 --ignoreGenomeQuality --clusterAlg singledRep结果处理
1.获得OTU members awk -F, NRFNR{a[$1]$1}NRFNR{if($1 in a ){print a[$1]*,$0}else{print $1,$0} } Wdb.csv Cdb.csv OTUs_members_Relationship.csv第一列为所有基因组名称代表基因组用星号表示第三列为OUT id相同id表示同一个OTU
2.比较两种生境基因组ANI距离分布密度图
Mdb.csv的第4列为两个基因组之间的相似性随后基于基因组metadata信息仅保留两个生境之间基因组的比较结果用于绘制密度图
# Mdb.csv示例
(base) [yutaomyosin data_tables]$ head Mdb.csv
genome1,genome2,dist,similarity
RS_4.bin.50.fa,RS_4.bin.50.fa,0.0,1.0
8_GM_sbin_oily_55.fa,RS_4.bin.50.fa,1.0,0.0
SRR13892588_me2_bin.130.fa,RS_4.bin.50.fa,1.0,0.0
7_SB_sbin_7_41.fa,RS_4.bin.50.fa,1.0,0.0
HTR7_me2_bin.71.fa,RS_4.bin.50.fa,1.0,0.0
RS_1.bin.2.fa,RS_4.bin.50.fa,1.0,0.0# 包含两种生境的基因组metadata示例
$ sed -i -e s/_genomic.fa//g;s/.fa//g Mdb_processed.csv# 存在AB及BA的情况
(base) [yutaomyosin data_tables]$ awk -F, $1RS_4.bin.50 $22_HM_cbin_9 Mdb_processed.csv
RS_4.bin.50,2_HM_cbin_9,1.0,0.0
(base) [yutaomyosin data_tables]$ awk -F, $2RS_4.bin.50 $12_HM_cbin_9 Mdb_processed.csv
2_HM_cbin_9,RS_4.bin.50,1.0,0.0# 保留其中一个即可
#echo -e a,b\nb,a\na,c|awk -F, {if(!a[$1,$2]!a[$2,$1])a[$1,$2]1}END{for(i in a){if(a[i])print i}}
$ awk -F, NR1{print}{if(!a[$1,$2,$3,$4]!a[$2,$1,$3,$4])a[$1,$2,$3,$4]1}END{for(i in a){if(a[i])print i}} Mdb_processed.csv Mdb_dereplicate.csv
(gtdbtk) [yutaomyosin data_tables]$ wc -l Mdb_processed.csv Mdb_dereplicate.csv9180901 Mdb_processed.csv4591967 Mdb_dereplicate.csv# 添加分组信息选择不在同一个组且相似性值不为0的结果
(base) [yutaomyosin data_tables]$ awk -F\t BEGIN{print MAG1,group1,MAG2,group2,dist,similarity}NRFNR{a[$1]$2}NRFNR{FS,;OFS,;if($1!~/genome/)print $1,a[$1],$2,a[$2],$3,$4} MAGs_metadata.tsv Mdb_dereplicate.csv|awk -F, $2!$4 $NF!0|head
MAG1,group1,MAG2,group2,dist,similarity
MAG_1907NJYTS1_metabat2_bin.13,Glacier,RS_4.bin.50,Cold seep,0.295981,0.704019
SRR5275918_metabat2_bin.57,Glacier,RS_4.bin.50,Cold seep,0.295981,0.704019
ISO_B518-1,Glacier,RS_4.bin.50,Cold seep,0.295981,0.704019
ISO_GM2-5-1,Glacier,RS_4.bin.50,Cold seep,0.295981,0.704019
MAG_KQGRI2.20.08_metabat2_bin.34,Glacier,8_GM_sbin_oily_55,Cold seep,0.295981,0.704019
MAG_18PLSC_vamb_S1C8874,Glacier,8_GM_sbin_oily_55,Cold seep,0.295981,0.704019
GCA_015751965.1_ASM1575196v1,Glacier,SRR13892588_me2_bin.130,Cold seep,0.295981,0.704019
ISO_B1382-2,Glacier,SRR13892588_me2_bin.130,Cold seep,0.295981,0.704019
ISO_N1997-1-1,Glacier,SRR13892588_me2_bin.130,Cold seep,0.295981,0.704019# 或可在最后去重
(gtdbtk) [yutaomyosin data_tables]$ awk -F, {if(!a[$1,$3]!a[$3,$1])a[$1,$3]1}END{for(i in a){if(a[i])print $0}} GGG_and_CSMD_ANI_similarity.csv
根据drep ANI 0.99的聚类结果文件dRep_0.99/data_tables/Widb.csv的第11列cluster_members绘制Pie图
Cdb.csv包含了每个OTU下的isolates原始数据dRep_0.99/data_tables/Widb.csv
[uhdata_tables]$ head Widb.csv
genome,score,completeness,contamination,strain_heterogeneity,size,N50,cluster,taxonomy,tax_confidence,cluster_members,closest_cluster_member,furthest_cluster_member,completeness_metric,contamination_metric
MAG_1611DDSC2_vamb_S1C744.fa,1.7918826841424995,NA,NA,NA,NA,NA,1_1,NA,NA,2,100.00,100.00,NA,NA
MAG_1908TGLC3_vamb_S1C1195.fa,1.9215229052672849,NA,NA,NA,NA,NA,1_2,NA,NA,11,100.00,99.85,NA,NA
MAG_1807DKMD2_vamb_S1C96.fa,1.9449308606290945,NA,NA,NA,NA,NA,2_1,NA,NA,5,99.44,99.39,NA,NA查看同一个OTU下的基因组来源
[yut.yut-PC] ➤ awk {print $2} tab.txt |sort -u|while read id; do type_num$(awk -v id$id $2id{split($1,a,_);print a[1]} tab.txt|sort -u|wc -l); otu_num$(awk -v id$id $2id tab.txt|wc -l);if [ $type_num -eq 2 ];then printf $id\t$type_num\t$otu_num\n;fi ;done
OTU212 2 21
OTU229 2 4
OTU259 2 7
OTU292 2 2
OTU349 2 16
OTU358 2 3
OTU359 2 7
OTU395 2 3work directory file-tree
workDirectory
./data
...../checkM/
...../Clustering_files/
...../gANI_files/
...../MASH_files/
...../ANIn_files/
...../prodigal/
./data_tables
...../Bdb.csv # Sequence locations and filenames
...../Cdb.csv # Genomes and cluster designations
...../Chdb.csv # CheckM results for Bdb
...../Mdb.csv # Raw results of MASH comparisons
...../Ndb.csv # Raw results of ANIn comparisons
...../Sdb.csv # Scoring information
...../Wdb.csv # Winning genomes
...../Widb.csv # Winning genomes checkM information
./dereplicated_genomes
./figures
./log
...../cluster_arguments.json
...../logger.log
...../warnings.txt其他输出结果
gtdbtk) [yutaoglobin F_Public_Database]$ ls 6023_MIMAG_MAGs_dRep95
data # 中间结果
data_tables # 数据统计结果
dereplicated_genomes # 去冗余后的基因组
figures # pdf图
log # (gtdbtk) [yutaoglobin 6023_MIMAG_MAGs_dRep95]$ ls data
ANImf_files # NUCMER 运行的每个基因组和相近基因组的delta.filtered结果
Clustering_files # 聚类的pickle文件secondary_linkage_cluster_1000.pickle
MASH_files # mash 结果
prodigal # prodigal预测结果包含fna,faa但不是所有基因组的关于代表基因组选择
dRep在使用checkM确定基因组质量以后会选择基因组完整度/污染度最优的基因组(Qality scorecomp-5cont)作为代表序列在cluster_score.pdf可以看到每个聚类及质量情况。 *代表OTU
使用
基因组去重复是识别基因组集合中“相同”的基因组并从每个冗余集合中除去“最佳”基因组之外的所有基因组的过程。
测试去重复
两个基因组相似性0.98
--S_algorithm fastANI -p PROCESSORS, --processors PROCESSORSthreads (default: 6)
--ignoreGenomeQualityDont 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)--S_algorithm {goANI,fastANI,ANImf,ANIn,gANI}Algorithm for secondary clustering comaprisons:fastANI Kmer-based approach; very fastANImf (DEFAULT) Align whole genomes with nucmer; filter alignment; compare aligned regionsANIn Align whole genomes with nucmer; compare aligned regionsgANI Identify and align ORFs; compare aligned ORFSgoANI Open source version of gANI; requires nsmimscan(default: ANImf)
-pa P_ANI, --P_ani P_ANIANI threshold to form primary (MASH) clusters(default: 0.9)-sa S_ANI, --S_ani S_ANIANI threshold to form secondary clusters (default:0.99)
(drep) [uhtest]$ tree -L 2 drep_out/
drep_out/
├── data
│ ├── ANImf_files
│ ├── Clustering_files
│ └── MASH_files
├── data_tables
│ ├── Bdb.csv
│ ├── Cdb.csv
│ ├── genomeInformation.csv 所有基因组信息N50,Size以及CheckM评估结果若执行checkm
│ ├── Mdb.csv
│ ├── Ndb.csv
│ ├── Sdb.csv
│ ├── Wdb.csv
│ └── Widb.csv
├── dereplicated_genomes 去重后的基因组
│ └── Ecoli_k12.fasta
├── figures 相关图片
│ ├── Clustering_scatterplots.pdf
│ ├── Cluster_scoring.pdf
│ ├── Primary_clustering_dendrogram.pdf
│ ├── Secondary_clustering_dendrograms.pdf
│ ├── Secondary_clustering_MDS.pdf
│ └── Winning_genomes.pdf
└── log├── cluster_arguments.json├── logger.log└── warnings.txt
de-replication
dRep dereplicate -p 20 outout_directory -g path/to/genomes/*.fasta
# 设置完整度/污染度
dRep dereplicate -p 50 -comp 0 -con 100000 dRep_out -g metabat2_bins/*fa positional arguments: work_directory Directory where data and output *** USE THE SAME WORK DIRECTORY FOR ALL DREP OPERATIONS *** -p PROCESSORS, --processors PROCESSORS threads (default: 6) -l LENGTH, --length LENGTH Minimum genome length (default: 50000) -comp COMPLETENESS, --completeness COMPLETENESS Minumum genome completeness (default: 75) -con CONTAMINATION, --contamination CONTAMINATION Maximum genome contamination (default: 25) CLUSTERING PARAMETERS: -pa P_ANI, --P_ani P_ANI ANI threshold to form primary (MASH) clusters (default: 0.9) -sa S_ANI, --S_ani S_ANI ANI threshold to form secondary clusters (default: 0.99) 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) 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) 参考
dRep doc Are these microbes the same? Konstantinidis, K. T., Ramette, A. Tiedje, J. M. The bacterial species definition in the genomic era. Philos. Trans. R. Soc. B Biol. Sci. 361, 1929–1940 (2006). Goris, J. et al. DNA-DNA hybridization values and their relationship to whole-genome sequence similarities. Int. J. Syst. Evol. Microbiol. 57, 81–91 (2007). Richter, M. Rosselló-Móra, R. Shifting the genomic gold standard for the prokaryotic species definition. Proc. Natl. Acad. Sci. 106, 19126–19131 (2009). Varghese, N. J. et al. Microbial species delineation using whole genome sequences. Nucleic Acids Res. 43, 6761–6771 (2015). Kim, M., Oh, H.-S., Park, S.-C. Chun, J. Towards a taxonomic coherence between average nucleotide identity and 16S rRNA gene sequence similarity for species demarcation of prokaryotes. Int. J. Syst. Evol. Microbiol. 64, 346–351 (2014). Luo, C. et al. Genome sequencing of environmental Escherichia coli expands understanding of the ecology and speciation of the model bacterial species. Proc. Natl. Acad. Sci. 108, 7200–7205 (2011). Kashtan, N. et al. Single-cell genomics reveals hundreds of coexisting subpopulations in wild Prochlorococcus. Science 344, 416–420 (2014). Snitkin, E. S. et al. Tracking a hospital outbreak of carbapenem-resistant Klebsiella pneumoniae with whole-genome sequencing. Sci. Transl. Med. 4, 148ra116–148ra116 (2012). Chase, J. et al. Geography and Location Are the Primary Drivers of Office Microbiome Composition. mSystems 1, (2016). Olm, M. R. et al. Identical bacterial populations colonize premature infant gut, skin, and oral microbiomes and exhibit different in situ growth rates. Genome Res. gr-213256 (2017). Olm, M. R., Brown, C. T., Brooks, B. Banfield, J. F. dRep: A tool for fast and accurate genome de-replication that enables tracking of microbial genotypes and improved genome recovery from metagenomes. bioRxiv (2017). doi:10.1101/108142 IF: NA NA NA
