在前两期的推送稿中,我们为大家介绍了宏基因组组装的基本原理和操作。基于组装序列,我们可以实现基因预测、物种注释、功能注释等相关分析,从而研究微生物菌群结构、菌属功能及作用机制。因此,本期我们将从基因预测的原理和操作两个部分出发,为大家介绍基于组装序列的基因预测。
1
基因预测原理
宏基因组基因预测一般包括同源预测和从头预测。同源预测是通过与基因的同源序列比对,从而获得与已知基因序列最大匹配,其预测依赖于已知的基因信息,且不能注释出在数据库中缺少功能相似性序列的基因和新基因,计算资源消耗过大,时间花费过长。而从头预测是根据给定的序列特征来预测,主要依赖于在编码区和非编码区拥有不同特征的信息,并在统计学上进行描述,构建概率模型,用以区别编码与非编码区。从头预测能够预测出已知的和未知的基因,且计算资源消耗小,时间花费少,常用软件包括:GeneMark,MetaGeneMark,MetaGene等。本期我们主要通过基于从头预测原理的MetaGeneMark来预测基因。
预测过程包括
2
基因预测实现
(1)软件
MetaGeneMark(预测的范围是细菌和古菌),下载地址:
http://exon.gatech.edu/license_download.cgi。
CD-HIT去除冗余序列,下载地址:http://www.bioinformatics.org/cd-hit。
(2)输入文件
SOAPdenovo-63mer对单个样本进行组装后,筛选出长度不小于500bp的scaftigs,得到结果文件sample1.cut500.scafSeq,sample2.cut500.scafSeq,sample3.cut500.scafSeq;SOAPdenovo-63mer对所有样本进行混合组装,筛选出长度不小于500bp的scaftigs,得到结果文件mix.cut500.scafSeq。文件格式如图所示,包括scaffold编号,长度及序列信息。
(3)基因预测
基于单个样本的基因预测
gmhmmp -a -d -f G -m MetaGeneMark_v1.modsample1.cut500.scafSeq -A sample1_protein.fasta -D sample1_nucleotide.fasta
gmhmmp -a -d -f G -m MetaGeneMark_v1.modsample2.cut500.scafSeq -A sample2_protein.fasta -D sample2_nucleotide.fasta
gmhmmp -a -d -f G -m MetaGeneMark_v1.modsample3.cut500.scafSeq -A sample3_protein.fasta -D sample3_nucleotide.fasta
基于混合组装的基因预测
gmhmmp -a -d -f G -m MetaGeneMark_v1.modmix.cut500.scafSeq -A mix_protein.fasta -D mix_nucleotide.fasta
参数说明:
-a 显示预测得到的基因的蛋白序列
-A 蛋白序列输出文件
-d 显示预测得到的基因的核酸序列
-D 核酸序列输出文件
-f显示输出格式,L=LST,G=GFF
-m 用于基因预测的模型文件,MetaGeneMark提供的MetaGeneMark_v1.mod适用于宏基因组预测
(4)基因去冗余
A. 将上一步得到的所有核酸序列(sample1_nucleotide.fasta,sample2_nucleotide.fasta,sample3_nucleotide.fasta,mix_nucleotide.fasta)合并成一个核酸序列total.gene.nucl.fasta; 将所有蛋白序列合并成一个total.gene.prot.fasta。
B. 蛋白质序列去冗余:
cd-hit-c 0.9 -n 5 -M 1600 -d 0 -T 8 -i total.gene.prot.fasta -o NonRundant.gene.prot.fasta
参数说明:
-c 相似性阈值,默认值为0.9
-n word长度值,基于word filter方法使用不同word长度值产生的去冗余水平不同,例如长度为2的word得到相似性在50%以上的序列,长度为3的word得到相似性在66.7%以上的序列
-M 分配的内存
-d 聚类信息文件中序列名长度,0代表显示完整序列名。
-T 线程数
-i 输入文件,fasta格式
-o 输出文件前缀,有两个输出文件,分别为fasta格式和以.clstr结尾的聚类信息文件。
输出文件:去冗余后的fasta格式文件NonRundant.gene.prot.fasta;以.clstr结尾的聚类信息文件NonRundant.gene.prot.fasta.clstr(见下图)。
注:每一个聚类组以“>”区分,且每个聚类组有不同的聚类序列,每个聚类序列中百分比代表该序列与代表序列的相似度,“*”代表该序列即为代表序列。
C. 核酸序列去冗余:
由于不同的核酸序列翻译后可能产生相同的蛋白质序列,因此对核酸序列的去冗余需要基于蛋白质序列。具体做法是:基于蛋白质序列文件与核酸序列文件中序列号的对应关系,根据去冗余后的蛋白质序列文件的序列号筛选核酸序列。
(5)基因丰度计算
使用soap将clean data比对到非冗余的核酸序列NonRundant.gene.nucl.fasta。根据soap比对结果,获得比对到某基因的read数。Soap比对工具的使用在之前的微信推送稿中介绍过,这里就不做说明。然后,基于read数获得基因的丰度值,计算方法如下:
S1 | S2 | S3 | |
G1 | . | . | . |
G2 | . | . | . |
G3 | . | . | . |
样本i中基因j表示为SiGj, 样本i中比对到基因j的read数表示为SiGj_read,基因j长度表示为Gj_len。
SiGj丰度值=(SiGj_read/Gj_len)/(SiG1_read/G1_len+SiG2_read/G2_len+SiG3_read/G3_len ...)
宏基因组基于组装的基因预测就介绍到这里,大家可以动起手来试一试,在下面一期小编将为大家带来宏基因组物种注释和功能注释部分,敬请期待!
供稿人:微生物事业部 王先月
联系客服