基因组组装或者宏基因组binning获得的基因组草图,首先需要评估其质量,包括基因组完整度、污染度、序列分布等信息。
基因组评估最常用的软件是CheckM(https://ecogenomics.github.io/CheckM/)。CheckM提供了一系列工具用于评估从分离培养、单细胞、宏基因组获得的基因组质量,可以根据基因组在参考基因组发育树中的位置来推断其精确的单拷贝标记基因集(lineage-specificmarker set),同时也提供数据库可用的基于分类学的基因集(taxonomic-specificmarker set)。CheckM利用基因的单拷贝性来有效的估计基因组完整度和污染,同时能绘制基因组关键特征(例如GC含量、编码率)的图像来评估基因组的质量。
CheckM安装,所需要依赖的环境如下所示:
代码语言:javascript复制Python3
HMMER (>=3.1b1)
prodigal (2.60 or >=2.6.1)
pplacer (>=1.1,https://github.com/matsen/pplacer)
安装方法如下所示:
代码语言:javascript复制pip3 install numpy
pip3 install matplotlib
pip3 install pysam
#如果已安装可以忽略以上步骤
pip3 install checkm-genome
下载数据库并设置数据库路径:
代码语言:javascript复制wget -c https://data.ace.uq.edu.au/public/CheckM_databases/checkm_data_2015_01_16.tar.gz
tar -zxvf checkm_data_2015_01_16.tar.gz
checkm data setRoot $PATH/checkm_data
CheckM主要的工具命令如下所示:
代码语言:javascript复制Lineage-specific marker set:
tree 将bins放入参考基因组发育树
tree_qa 评估每个bin里的系统发育标记基因
lineage_set 推断每个bin的标记基因集
Taxonomic-specific marker set:
taxon_list 列出数据库可用的不同分类水平列表
taxon_set 指定一个分类水平制作基因集
以上为两种确定基因组标记基因集合(marker set)的方法,使用将bins放入系统发育树依据系统发育关系推断的标记集合为lineage-specificmarker sets,使用依据分类系统产生的为taxonomic-specificmarker set。
代码语言:javascript复制Apply marker set to genome bins:
analyze 识别bins中的标记基因
qa 评估bins完整度和污染度
下面两个命令为上述命令的集合流程:
代码语言:javascript复制 lineage_wf 运行tree、lineage_set、analyze、qa
taxonomy_wf 运行taxon_set、analyze、qa
一般情况下推荐使用基于系统发育的流程,其使用方法如下所示:
代码语言:javascript复制checkm lineage_wf <bin folder> <output folder>
其中bin folder为含有bins序列的路径,output folder为结果文件路径名称(程序会自动创建文件夹),如果所获得的draft基因组都是属于某个已知分类单元,那么使用基于分类学的方法更加便捷,使用方法如下所示:
代码语言:javascript复制checkm taxonomy_wf <rank> <taxon> <bin folder> <output folder>
其中rank为分类层级例如phylum,taxon为分类单元例如Cyanobacteria。下面使用lineage_wf流程进行分析,如下所示:
代码语言:javascript复制nohup checkm lineage_wf -t 20 -x fa --nt --tab_table -f bins_qa.txt metabat_bins bins_qa_result &
其中-x指定bins文件的拓展名,输入路径中其他拓展名的文件将被忽略;--nt输出每个bin中的基因序列(调用prodigal软件进行预测);-f将默认输出到标准输出的评估结果储存到指定结果文件;--tab_table结果文件中表格形式的结果以tab分隔。
运行结束后生成的bins_qa.txt结果文件中包含bin的谱系、基因组基因数目、marker基因数目、完整度、污染度等信息,如下所示:
在结果路径bins_qa_result/bins中为每个bin预测的基因序列,在bins_qa_result/storage中则为每个bin详细的评估信息,其中bin_stats.analyze.tsv为每个bin基础统计信息,bin_stats.tree.tsv为每个bin在发育树中的位置,bin_stats_ext.tsv为每个bin对应的marker基因集,marker_gene_stats.tsv为每个bin的序列上marker基因比对信息。
除了综合评估外,CheckM提供了一系列工具来计算基因组特征,具体如下所示:
代码语言:javascript复制unbinned 识别没有被分装(unbinned)的序列
coverage 计算序列的coverage
tetra 计算每条序列的四核苷酸频率
profile 计算map到每个bin的reads的百分率,可用比较bins丰度
join_tables 将tab分割的不同bin信息表文件整合
ssu_finder 识别序列中的核糖体小亚基RNA(SSU rRNAs),也即16S/18S
CheckM还提供了一系列作图工具,用于bins质量可视化,如下所示:
代码语言:javascript复制bin_qa_plot:绘制bin完整度、污染度和异质性条形图
gc_plot:绘制每个bin的不同序列GC含量分布直方图及误差图
coding_plot:绘制每个bin序列的编码密度(coding density,CD)直方图及误差图
tetra_plot:绘制bin每条序列与bin平均四核苷酸频率的距离(tetranucleotide distance,TD)直方图及误差图
dist_plot:将以上三个图形绘制在一起
其中dist_plot使用方法如下所示:
代码语言:javascript复制checkm dist_plot [Options] out_folder bin_folder plot_folder tetra_profile dist_value
out_folder CheckM评估bins的结果文件夹,也即前面生成的bins_qa_result
tetra_profile tetra命令计算的contigs序列四核苷酸频率
plot_folder 输出图像的文件夹,无需事先创建
dist_value 概率分布距离,也即展示contigs序列的置信区间,用于误差图
--image_type 输出图片格式,可选eps、pdf、png、ps、svg,默认为png
--dpi 输出图片的DPI,默认为600
--font_size 输出图片字体大小,默认为8
-x, --extension bins序列文件的拓展名,默认为fna,文件夹中其他后缀的文件将被忽略
--width 输出图片的宽度,默认为6.5
--height 输出图片的高度,默认为8
-a, --gc_window_size 计算GC含量时滑窗大小(window size),默认为5000
-b, --td_window_size 计算TD时滑窗大小,默认为5000
-c, --cd_window_size 计算CD时滑窗大小,默认为10000
-1, --gc_bin_width 图像中GC bars宽度,默认为0.01
-2, --td_bin_width 图像中TD bars宽度,默认为0.01
-3, --cd_bin_width 图像中CD bars宽度,默认为0.01
-q, --quiet 压缩输出结果
下面绘制bins质量评估图像,如下所示:
代码语言:javascript复制checkm dist_plot --image_type pdf -x fa bins_qa_result metabat_bins checkm_plots ../checkm_tetra.out 95
评估结果如下所示:
bin_qa_plot使用方法如下所示:
代码语言:javascript复制checkm bin_qa_plot --image_type pdf -x fa bins_qa_result metabat_bins checkm_qa_plots
部分结果如下所示:
不同的颜色分别代表单拷贝、丢失、杂合与污染的marker基因,每一个bar代表一个marker,多拷贝基因之间氨基酸匹配(amino acid identity,AAI)大于90%被认为是杂合的(同一个物种不同株的等位基因),而AAI小于90%被认为是其他物种污染。
为了进一步评估每个bin的拼接程度,可以绘制Nx图(x=0.5即为基因组评估的N50),如下所示:
代码语言:javascript复制checkm nx_plot --image_type pdf -x fa --font_size 12 metabat_bins checkm_Nx_plots
评估结果如下所示:
另外两个相似的作图命令:
代码语言:javascript复制len_plot:每个bin累积序列长度
len_hist:每个bin序列长度直方图
使用marker_plot命令可绘制marker基因在序列中的位置,如下所示:
代码语言:javascript复制checkm marker_plot --image_type pdf -x fa --font_size 10 bins_qa_result metabat_bins checkm_marker_plots
部分结果如下所示:
根据CheckM评估结果,可进行后续的基因组质量优化。
END