16S  rRNA作为微生物的特征片段,在鉴定微生物物种的过程中起着重要的作用。我们在做宏基因组分析的时候可以找一下scaffolds中的16S  rRNA,然后跟数据库比对,看里面都有什么物种;或者我们直接对环境样品进行16S  rRNA 高通量测序,然后来获知该环境样品的群落组成。总之16srRNA的分析很重要。本博文将围绕16S  rRNA的预测与鉴定两部分展开。

16s

一、预测

Rnammer、Infernal、Usearch 可以用来预测rRNA。

1.RNAmmer

网页版偶尔不好使,可以本地化该软件

RNAmmer本地化

1.1安装要求:
Unix系统,perl,hmmsearch(需要安装hmmer-2.2g版本,使用最新版本会出错)

1.2安装

安装hmmer-2.2g

$ wget ftp://selab.janelia.org/pub/software/hmmer/2.2g/hmmer-2.2g.tar.gz
$ tar zxf hmmer-2.2g.tar.gz
$ cd hmmer-2.2g
$ ./configure
$ make;
$make install(这一步应该是改环境变量吧,这一步我没有弄。有一个哥们提到:到make这一步都完全正常,只是make install 时出现了错误,但是最后binaries文件夹中出现了一系列程序,直接将程序cp到/usr/local/bin/目录种就可以用了。未验证他的这种说法,保留意见。)

安装RNAmmer
该软件需要使用edu邮箱去申请。下载地址:http://www.cbs.dtu.dk/cgi-bin/nph-sw_request?rnammer (需要学校的邮箱来认证)
下载,解压缩,然后进入rnammer修改部分地址, 一个是你安装Rnammer路径,一个是hmmsearch的路径,安装在哪里就改成哪里就可以了

(两处修改)
my $INSTALL_PATH = "/usr/cbs/bio/src/rnammer-1.2";
修改为my $INSTALL_PATH = "/sam/hmmer/binaries/rnammer";
HMMSEARCH_BINARY = "/usr/cbs/bio/bin/linux64/hmmsearch";
修改为HMMSEARCH_BINARY = "/sam/hmm/binaries /hmmsearch";
(一哥们安装经验,为验证,放在这供参考: 申请到的RNAmmer软件包是现有的程序,所以最佳的方式就是直接将它解压到/usr/local/bin/目录中就好了。修改一下rnammer程序。首先应该修改rnammer的权限,为了方便我直接用 chmod 777 rnammer 命令改成了最高权限了。然后用vi nammer 对其进行修改,如博主所述,修改的两个地方。一个是你安装nammer路径,一个是hmmsearch的路径,安装在哪里就改成哪里就可以了。ps: /usr/local/bin/ 是默认的环境变量路径,所以安装在这个目录中的软件你在任何一个目录下都可以使用,就不需要每次都加上程序的绝对路径了。)

检测是否正常运行

cd /sam/rnammer/example;
../rnammer -S bac -m lsu,ssu,tsu -xml ecoli.xml -gff ecoli.gff -h ecoli.hmmreport < ecoli.fsa
(建议复制这个example文件夹,这样在这个文件夹里面运行结束后,可以比较两个的结果)

1.3. RNAmmer的运行与参数

$ rnammer [options] sequence.fasta
-S 指定输入序列的物种所属的界: arc bac 或 euk
-gff 生成的gff文件结果
-m 所需要预测的moleculers: 'tsu' for 5/8s rRNA, 'ssu' for 16/18s rRNA, 'lsu' for 23/28s rRNA。如果全部进行预测,则该参数后为'tsu,ssu,lsu'。
-multi 并行运算,预测正反两条链上所有的moleculers。最多并行运行6个计算。使用该参数,则不需要上一个参数。
-f 生成的rRNA的fasta结果文件
-h 生成的hmm报告结果
-gff 生成的rRNA的gff2文件
-x 生成的xml结果文件

对真核生物基因组进行rRNA预测的一个示例命令:

$ /opt/biosoft/rnammer-1.2/rnammer -S euk -multi -f rRNA.fasta\
-h rRNA.hmmreport -xml rRNA.xml -gff rRNA.gff2 genome.fasta

细菌的一条染色体,运行一下命令即可

rnammer -S bac -multi -f rRNA.fasta -h rRNA.hmmreport -xml rRNA.xml -gff rRNA.gff2 < genome.fasta

1.4. RNAmmer的结果
rRNA.fasta,rRNA.gff2,rRNA.hmmreport,rRNA.xml。

2. infernal

官网http://infernal.janelia.org/

3.usearch寻找16S rRNA

当然,也可参照我之前的博客http://blog.sina.com.cn/s/blog_670445240101njwi.html

二.鉴定

预测出来的16S rRNA用RDP即可知道该rRNA是属于哪个物种的,由此我们推测该cluster是什么物种。

Rdp (ribosomal database project)用于对16S rRNA或者真菌LSU序列的分类地位上的鉴定。RDP naïve Bayesian rRNA Classifier的一种分层分类方法。
官网服务器:http://rdp.cme.msu.edu/classifier/classifier.jsp
在官网上将你的序列上传上去,如果是16S rRNA,就选择模式:16S rRNA training set 9,然后submit,一会的功夫就会出来结果。

两种模式:
16S rRNA training set 9
Fungal LSU

计算方式:
这个方法不需要比对。8个碱基一组(称为word),这样就有6万多种排列组合。对已知的880个属(现在数目应该大于这个)按照这个这种方式统计6万多种排列的比例情况。 上传的序列的组合数也求出来,然后随即选择word,跟参考库的来计算,(具体怎么计算,这个就是贝叶斯算法的核心罗,先不管),进行100次,100次计算都在同一个分类地位的数目,则为bootstrap cutoff。分类结果的后面都会跟上一个bootstrap cutoff。设定一个阈值,root的置信阈值一般选择的是80%。这个bootstrap cutoff又叫confidence estimate

上传序列的格式:
Fasta, GenBank and EMBL格式
长度至少为50个碱基,序列数要小于100000条。序列数超过100000条,联系
rdpstaff@msu.edu

置信阈值
通过100次抽样,分配到该分类地位的次数,则为boot cutoff ,极为最后结果分类地位后面跟着的那个百分数.
对于序列长度小于250bp,长于50bp.bootstrap cutoff 为50%则认为可以反映序列属的水平

结果:
序列后面跟着的“-”意味是根据序列反向来得出的结果。
根据置信阈值,如果小于50%,则在前面的分类地位前加上unclassfied,写了一个脚本。
allrank_rep_otu_nosingle.fasta_classified.txt为RDP反馈的结果,有需求脚本rdpformat.pl留下邮箱。

sed '1,7d' allrank_rep_otu_nosingle.fasta_classified.txt|sed 's/%//g' |sed 's/"//g' >abc;
perl /sam/usefulscript/rdp_format.pl abc abd;

PS:
这里面还有一点我不是太确定,所有的分类地位后面都有一个bootstrap cutoff,大于50%则是一类,小于50%则是unclassfied,是这样吗? 我怎么感觉仅仅是针对属来说的呢?在这里留下一个悬念。


参考资料:
http://rdp.cme.msu.edu/classifier/class_help.jsp#model
http://www.chenlianfu.com/?p=1979
http://www.cbs.dtu.dk/services/doc/rnammer-1.2.readme

发表评论

电子邮件地址不会被公开。 必填项已用*标注