生信分析、meta分析、数据挖掘
TCGA、GEO、SEER、Oncomine

宏基因组分析——基因预测篇

在前两期的推送稿中,我们为大家介绍了宏基因组组装的基本原理和操作。基于组装序列,我们可以实现基因预测、物种注释、功能注释等相关分析,从而研究微生物菌群结构、菌属功能及作用机制。因此,本期我们将从基因预测的原理和操作两个部分出发,为大家介绍基于组装序列的基因预测。

 基因预测原理 

宏基因组基因预测一般包括同源预测和从头预测。同源预测是通过与基因的同源序列比对,从而获得与已知基因序列最大匹配,其预测依赖于已知的基因信息,且不能注释出在数据库中缺少功能相似性序列的基因和新基因,计算资源消耗过大,时间花费过长。而从头预测是根据给定的序列特征来预测,主要依赖于在编码区和非编码区拥有不同特征的信息,并在统计学上进行描述,构建概率模型,用以区别编码与非编码区。从头预测能够预测出已知的和未知的基因,且计算资源消耗小,时间花费少,常用软件包括:GeneMarkMetaGeneMarkMetaGene等。本期我们主要通过基于从头预测原理的MetaGeneMark来预测基因。

预测过程包括

 


2  基因预测实现


(1)软件

MetaGeneMark(预测的范围是细菌和古菌),下载地址:

http://exon.gatech.edu/license_download.cgi。

CD-HIT去除冗余序列,下载地址:http://www.bioinformatics.org/cd-hit。

2)输入文件

SOAPdenovo-63mer对单个样本进行组装后,筛选出长度不小于500bpscaftigs,得到结果文件sample1.cut500.scafSeqsample2.cut500.scafSeqsample3.cut500.scafSeqSOAPdenovo-63mer对所有样本进行混合组装,筛选出长度不小于500bpscaftigs,得到结果文件mix.cut500.scafSeq。文件格式如图所示,包括scaffold编号,长度及序列信息。

(3)基因预测

基于单个样本的基因预测

gmhmmp -a -d -f G -m MetaGeneMark_v1.mod  sample1.cut500.scafSeq -A sample1_protein.fasta -D sample1_nucleotide.fasta

gmhmmp -a -d -f G -m MetaGeneMark_v1.mod  sample2.cut500.scafSeq -A sample2_protein.fasta -D sample2_nucleotide.fasta

gmhmmp -a -d -f G -m MetaGeneMark_v1.mod  sample3.cut500.scafSeq -A sample3_protein.fasta -D sample3_nucleotide.fasta

基于混合组装的基因预测

gmhmmp -a -d -f G -m MetaGeneMark_v1.mod  mix.cut500.scafSeq -A mix_protein.fasta -D mix_nucleotide.fasta

参数说明:

-a 显示预测得到的基因的蛋白序列

-A 蛋白序列输出文件

-d 显示预测得到的基因的核酸序列

-D 核酸序列输出文件

-f  显示输出格式,L=LSTG=GFF

-m 用于基因预测的模型文件,MetaGeneMark提供的MetaGeneMark_v1.mod适用于宏基因组预测

(4)基因去冗余

A. 将上一步得到的所有核酸序列(sample1_nucleotide.fastasample2_nucleotide.fastasample3_nucleotide.fastamix_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)基因丰度计算

使用soapclean data比对到非冗余的核酸序列NonRundant.gene.nucl.fasta。根据soap比对结果,获得比对到某基因的read数。Soap比对工具的使用在之前的微信推送稿中介绍过,这里就不做说明。然后,基于read数获得基因的丰度值,计算方法如下:


S1

S2

S3

G1

.

.

.

G2

.

.

.

G3

.

.

.

样本i中基因j表示为SiGj,  样本i中比对到基因jread数表示为SiGj_read,基因j长度表示为Gj_len

SiGj丰度值=SiGj_read/Gj_len/SiG1_read/G1_len+SiG2_read/G2_len+SiG3_read/G3_len …

宏基因组基于组装的基因预测就介绍到这里,大家可以动起手来试一试,在下面一期小编将为大家带来宏基因组物种注释和功能注释部分,敬请期待!


供稿人:微生物事业部 王先月


推荐阅读
  1. 微生物分析——R语言绘图篇1

  2. 微生物分析——R语言绘图篇2

  3. 微生物分析——R语言绘图篇3

  4. 六步教会你基因组组装!

  5. 手把手教您学会宏基因组组组装

生信圈致力于每天推送生物信息干货,让大家了解生信行业。旨在通过更多的交流促进行业的发展。我们一直在寻找志同道合的伙伴!投稿邮箱:bioinfor_club@163.com

生信圈

微信ID:bioinfor-club

1.点击历史信息,查看更多内容

2.长按右侧二维码,关注更多生物信息干货

长按二维码关注

赞(0) 打赏
未经允许不得转载:医学SCI科研之家 » 宏基因组分析——基因预测篇
分享到: 更多 (0)

评论 抢沙发

  • 昵称 (必填)
  • 邮箱 (必填)
  • 网址

meta分析、生信分析

meta、生信交流群综合科研交流群