【6.3.1】16SrRNA的预测与鉴定
16S rRNA作为微生物的特征片段，在鉴定微生物物种的过程中起着重要的作用。我们在做宏基因组分析的时候可以找一下scaffolds中的16S rRNA，然后跟数据库比对，看里面都有什么物种；或者我们直接对环境样品进行16S rRNA 高通量测序，然后来获知该环境样品的群落组成。总之16srRNA的分析很重要。本博文将围绕16S rRNA的预测与鉴定两部分展开。
一、预测
Rnammer、Infernal、Usearch 可以用来预测rRNA。
1.1. RNAmmer
网页版偶尔不好使，可以本地化该软件
RNAmmer本地化
1.1.1 安装要求
Unix系统，perl，hmmsearch(需要安装hmmer-2.2g版本，使用最新版本会出错)
1.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.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.1.4. RNAmmer的结果
rRNA.fasta，rRNA.gff2，rRNA.hmmreport，rRNA.xml。
1.2. infernal
官网http://infernal.janelia.org/
1.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,是这样吗？ 我怎么感觉仅仅是针对属来说的呢？在这里留下一个悬念。
参考资料
