本地使用Rfam 12.0+
Rfam是用来鉴定non-coding RNAs的数据库，常用于注释新的核酸序列或者基因组序列。
最新版本的Rfam (12.2)，包含2588个RNA家族，其在线网站提供了便捷的查询使用功能，http://rfam.xfam.org/，尤其是对小批量数据。
对已经注释好的物种，建议在ENSEMBLE或UCSC直接下载官方的注释文件，直接下载GTF或使用bioMart或TableBrowser都可。具体可在本博客或微信公众号后台回复关键字获取使用方法。
最后确定要在本地使用了，如果是用之前的版本，请在线搜索，有几篇教程可用。如果有新版本，请参考此教程。
下载infernal软件
-
官网下载：infernal软件本来应该在http://eddylab.org/infernal/下载，但这个网站最近打不开了。
-
GitHub下载
# 如果不会或不能使用git，拷贝链接，在GitHub下载压缩文件，再按顺序解压即可 # https://github.com/EddyRivasLab git clone https://github.com/EddyRivasLab/infernal.git infernal cd infernal git clone https://github.com/EddyRivasLab/easel.git git clone https://github.com/EddyRivasLab/hmmer.git #如果当前目录有aclocal.m4，则不需要执行 #ln -s `pwd`/easel/aclocal.m4 `pwd`/ ln -s `pwd`/easel/aclocal.m4 hmmer # 如果没有autoconf，找管理员配置，或查看软件安装 autoconf (cd easel; autoconf) (cd hmmer; autoconf) ./configure --prefix=`pwd`/../infernal_bin make make install cd easel; make install cd ../../infernal_bin/bin ls # 就可以看到安装的程序了，加入环境变量即可，本站搜索环境变量查看具体配置方法 # 或微信公众号 生信宝典 后台回复 环境变量 查看 export PATH=${PATH}:`pwd`
下载数据库
wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/12.2/Rfam.cm.gz
gunzip Rfam.cm.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/12.2/Rfam12.2.claninfo
# 使用 Infernal的程序cmpress索引Rfam.cm
# 大约需要一分钟
cmporess Rfam.cm
# 输出为
Working... done.
Pressed and indexed 2588 CMs and p7 HMM filters (2588 names and 2588 accessions).
Covariance models and p7 filters pressed into binary file: Rfam.cm.i1m
SSI index for binary covariance model file: Rfam.cm.i1i
Optimized p7 filter profiles (MSV part) pressed into: Rfam.cm.i1f
Optimized p7 filter profiles (remainder) pressed into: Rfam.cm.i1p
确定待查询的核苷酸序列或基因组的大小，作为后续命令的参数
# 大约1分钟
esl-seqstat my-genome.fa
其输出结果中有一行，类似于Total # of residues: 3000000是我们需要的。考虑到基因组为双链和下一步用到的参数的单位为Million，我们使用公式3000000 * 2 / 1000000计算得出结果为6，作为下一步参数-Z的值.
使用cmscan程序注释基因组的ncRNAs
# Rfam12.2.claninfo 为下载的claninfo文件，需提供所在路径
# Rfam.cm 下载的cm文件
# my-genome.fa 待查询序列
# my-genome.cmscan 输出结果
# my-genome.tblout 有一个输出结果
# 对500M大小的输入序列，48线程，需要7个小时，最好放入后台
cmscan -Z `esl-seqstat my-genome.fa | awk '{if($0~/^Total/) print int($4/2000000);}''` --cut_ga --rfam --nohmmonly --tblout my-genome.tblout --fmt 2 --clanin Rfam12.2.claninfo Rfam.cm my-genome.fa > my-genome.cmscan
参数解释:
-Z: 查询序列的大小，以M为单位。由esl-seqstat算出或自己写程序计算，记得乘以2，除以10^6.
--cut_ga: 输出不小于Rfam GA阈值的结果。这是Rfam认证RNA家族的阈值，不低于这个阈值的序列得分被认为是真同源序列。The bit score gathering threshold (GA cutoff), set by Rfam curators when building the family. All sequences that score at or above this threshold will be included in the full alignment and are believed to be true homologs to the model.
--rfam: run in “fast” mode, the same mode used for Rfam annotation and determination of GA thresholds.
--nohmmonly: all models, even those with zero basepairs, are run in CM mode (not HMM mode). This ensures all GA cutoffs, which were determined in CM mode for each model, are valid.
--tblout: table输出。
--fmt 2: table输出的一种格式。
--clanin: 下载的clan信息。This file lists which models belong to the same clan. Rfam clans are groups of models that are homologous and therefore it is expected that some hits to these models will overlap. For example, the LSU rRNA archaea and LSU rRNA bacteria models are both in the same clan.
结果解释
my-genome.cmscan
cmscan命令的标准输出，使用>重定向。
- 结果的第一部分是运行的命令的参数记录
-
第二部分是查询的FASTA文件中每个序列的top hits, 根据
E-value排序。Query: scaffold12 [L=2429009] Hit scores: rank E-value score bias modelname start end mdl trunc gc description ---- --------- ------ ----- --------- ------- ------- --- ----- ---- ----------- (1) ! 5.4e-20 86.6 0.0 mir-166 961624 961752 + cm no 0.43 - (2) ! 9.6e-14 70.9 0.0 tRNA 2369877 2369805 - cm no 0.59 - (3) ! 3.7e-13 68.8 0.0 tRNA 1344161 1344232 + cm no 0.53 - (4) ! 2.3e-12 65.9 0.0 tRNA 973203 973274 + cm no 0.54 - (5) ! 5.8e-12 64.5 0.0 tRNA 1241524 1241595 + cm no 0.50 -E-value: 统计显著性，依赖于查询数据库的大小。score: Log-odds得分，独立于查询序列数据库的大小。在使用了--cut_ga后所有输出的结果都是高于Rfam GA得分的。modelname: Rfam家族或模型的名字。start,stop: 查询序列匹配的区域。后面跟着的是链的信息，对于+，起始位置小于终止位置；对于-，其实位置大于终止位置。
互有重叠的查询区域可能会匹配到不同的Rfam家族或模型。推荐保留具有最低的E-value，最高的bit scaore的部分。
- 结果的下一部分是比对结果。具体可查看文后的参考网址链接的内容。
my-genome.tblout
表格输出包含了cmscan标准输出的大部分内容，并且便于对结果的进一步处理。
#idx target name accession query name accession clan name mdl mdl from mdl to seq from seq to strand trunc pass gc bias score E-value inc olp anyidx afrct1 afrct2 winidx wfrct1 wfrct2 description of target
#--- -------------------- --------- -------------------- --------- --------- --- -------- -------- -------- -------- ------ ----- ---- ---- ----- ------ --------- --- --- ------ ------ ------ ------ ------ ------ ---------------------
1 tRNA RF00005 scaffold20 - CL00001 cm 1 71 1399503 1399576 + no 1 0.57 0.0 58.8 2.4e-10 ! * - - - - - - -
2 tRNA RF00005 scaffold20 - CL00001 cm 1 71 186338 186267 - no 1 0.54 0.0 55.4 2e-09 ! * - - - - - - -
1 mir-166 RF00075 scaffold12 - - cm 1 126 961624 961752 + no 1 0.43 0.0 86.6 5.4e-20 ! * - - - - - - -
2 tRNA RF00005 scaffold12 - CL00001 cm 1 71 2369877 2369805 - no 1 0.59 0.0 70.9 9.6e-14 ! * - - - - - - -
3 tRNA RF00005 scaffold12 - CL00001 cm 1 71 1344161 1344232 + no 1 0.53 0.0 68.8 3.7e-13 ! * - - - - - - -
4 tRNA RF00005 scaffold12 - CL00001 cm 1 71 973203 973274 + no 1 0.54 0.0 65.9 2.3e-12 ! * - - - - - - -
5 tRNA RF00005 scaffold12 - CL00001 cm 1 71 1241524 1241595 + no 1 0.50 0.0 64.5 5.8e-12 ! * - - - - - - -
1 Plant_SRP RF01855 scaffold15 - CL00003 cm 1 305 1511627 1511325 - no 1 0.62 0.7 249.5 1.1e-70 ! ^ - - - - - - -
每一行有27列，比较关键的是`target name`, `accession`, `query name`, `seq from`, `seq to`, `strand`, `E-value`, `score`。
`olp`列表示查询序列的重叠信息，`*`表示同一条链上，不存在与此查询序列重叠的序列也在Rfam数据库有匹配，这是需要保留的查询序列。`^`表示同一条链上，不存在比此查询序列与Rfam数据库匹配更好的序列，也需要保留。`=`表示同一条链上，存在比此查询序列与Rfam数据库匹配更好的序列，应忽略。
可通过下面的命令获取最终结果。
```bash
grep -v '=' my-genome.tblout >my-genome.deoverlapped.tblout
结果解析
首先转换下结果格式，提取必须得列和非重叠区域或重叠区域中得分高的区域。
awk 'BEGIN{OFS="\t";}{if(FNR==1) print "target_name\taccession\tquery_name\tquery_start\tquery_end\tstrand\tscore\tEvalue"; if(FNR>2 && $20!="=" && $0!~/^#/) print $2,$3,$4,$10,$11,$12,$17,$18; }' my-genome.tblout >my-genome.tblout.final.xls
统计预测出的ncRNA的类型。
首先下载Rfam家族的注释，点击http://rfam.xfam.org/search#tabview=tab5，选择所有复选框，提交，把得到的表格拷贝下来，整理成TAB键分割的格式。并吧第三列拆开，取出类型, 存储为 Rfam_anno.txt。
Accession ID Type Description class
RF00001 5S_rRNA Gene; rRNA 5S ribosomal RNA rRNA
RF00002 5_8S_rRNA Gene; rRNA 5.8S ribosomal RNA rRNA
RF00003 U1 Gene; snRNA; splicing U1 spliceosomal RNA snRNA
awk 'BEGIN{OFS=FS="\t"}ARGIND==1{a[$2]=$5;}ARGIND==2{type=a[$1]; if(type=="") type="Others"; count[type]+=1;}END{for(type in count) print type, count[type];}' Rfam_anno.txt my-genome.tblout.final.xls
最终输出
rRNA 61
snRNA 397
sRNA 1
Others 55
antisense 1
tRNA 427
miRNA 95
ribozyme 1
riboswitch 1
