CheckM:基因组质量评估

2022-05-05 13:44:15 浏览数 (1)

基因组组装或者宏基因组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

0 人点赞