一：下载安装

wget http://eddylab.org/infernal/infernal-1.1.2-linux-intel-gcc.tar.gz
tar xf infernal-1.1.2-linux-intel-gcc.tar.gz
cd infernal-1.1.2
./configure  --prefix=.
make
make install
cd easel
make install
ls #可以用的软件
export PATH=${PATH}:/your_path/infernal-1.1.2/bin #添加环境变量

二：下载数据库

mkdir Rfam && cd Rfam    
wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/14.1/Rfam.cm.gz
gunzip Rfam.cm.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/14.1/Rfam14.1clanin
# 使用 Infernal的程序cmpress索引Rfam.cm # 大约需要一分钟
cmporess Rfam.cm


# 输出为
Working...    done.
Pressed and indexed 3016 CMs and p7 HMM filters (3016 names and 3016 accessions).
Covariance models and p7 filters pressed into binary file:  Rfam14.1.cm.i1m
SSI index for binary covariance model file:                 Rfam14.1.cm.i1i
Optimized p7 filter profiles (MSV part)  pressed into:      Rfam14.1.cm.i1f
Optimized p7 filter profiles (remainder) pressed into:      Rfam14.1.cm.i1p

三： 确定待查询的核苷酸序列或基因组的大小，作为后续命令的参数

esl-seqstat my-genome.fa

Z值

其输出结果中有一行，Total # of residues: 218502668是我们需要的。考虑到基因组为双链和下一步用到的参数的单位为Million。
计算公式为：Z = total * 2 * CMmumber/106
因此还要计算CM database中的模型的数量（Z=2636016）
在Rfam14.0版本中，使用

五：运行程序

# Rfam14.1.clanin 为下载的claninfo文件，需提供所在路径
# Rfam.cm 下载的cm文件
#nep_genome.fa 待查询序列
# nep.cmscan 输出结果
# nep.tblout 表格形式输出结果
# 对500M大小的输入序列，48线程，需要7个小时，最好放入后台
cmscan -Z 2636016 --cut_ga --rfam --nohmmonly --tblout mep_genome.tblout --fmt 2 --cpu 24 --clanin Rfam14.1.clanin Rfam.cm nep_genome.fa > nep.cmscan

建议加线程，别把服务器跑满了！！

六：输出结果

nep.tblout

#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    mir-7                RF00053   contig_1000          -         -          cm        1       88    51272    51359      +    no    1 0.38   0.0   58.2   1.9e-06  !   *       -      -      -      -      -      - mir-7 microRNA precursor
2    MIR1122              RF00906   contig_1000          -         -          cm        1      135   137126   137273      +    no    1 0.20  26.4   37.6       4.1  ?   *       -      -      -      -      -      - microRNA MIR1122
1    tRNA                 RF00005   contig_10006         -         CL00001    cm        1       71      616      707      +    no    1 0.53   0.0   49.7   0.00039  !   *       -      -      -      -      -      - tRNA
1    sno_ZL1              RF02723   contig_10010         -         -          cm       36      107   101457   101528      +    no    1 0.28   0.1   31.7       3.2  ?   *       -      -      -      -      -      - Small nucleolar RNA ZL1
1    Histone3             RF00032   contig_10016         -         -          cm        1       46   104440   104394      -    no    1 0.30   0.0   33.2     0.044  ?   *       -      -      -      -      -      - Histone 3' UTR stem-loop
1    LSU_rRNA_bacteria    RF02541   contig_10019         -         CL00112    cm        1     2925    51575    54339      +    no    1 0.44  53.8 2654.9         0  !   ^       -      -      -      -      -      - Bacterial large subunit ribosomal RNA
2    LSU_rRNA_archaea     RF02540   contig_10019         -         CL00112    cm        1     2990    51574    54338      +    no    1 0.44  54.3 1799.6         0  !   =       1  1.000  1.000      "      "      " Archaeal large subunit ribosomal RNA
## 每一行有27列，比较关键的是`target name`, `accession`, `query name`, `seq from`, `seq to`, `strand`, `E-value`, `score`。
`olp`列表示查询序列的重叠信息，`*`表示同一条链上，不存在与此查询序列重叠的序列也在Rfam数据库有匹配，这是需要保留的查询序列。`^`表示同一条链上，不存在比此查询序列与Rfam数据库匹配更好的序列，也需要保留。`=`表示同一条链上，存在比此查询序列与Rfam数据库匹配更好的序列，应忽略。
## 可通过下面的命令获取最终结果。
```bash
grep -v '=' nep.tblout >nep .deoverlapped.tblout

######################################################################

结果解析

#####################################################################

1.格式转换

得到tblout.final，转换下结果格式，提取必须得列和非重叠区域或重叠区域中得分高的区域。
=表示同一条链上，存在比此查询序列与Rfam数据库匹配更好的序列，应忽略。
去掉#开头的注释行 提取所需要的列。初始将分隔符设定为制表符。
如果是第一行，打印行名称 第3行开始，如果20列不是等号，且不以井号开始，打印2,3,4等列输出

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; }' nep.tblout >nep_genome.tblout.final.xls

2. 下载Rfam注释

首先下载Rfam家族的注释，点击Rfam网址，选择所有复选框，提交，把得到的表格拷贝下来，整理成TAB键分割的格式。将第三列的第二个分号信息提取出来，这就是各ncRNA class

cat rfamannotation.txt | awk 'BEGIN {FS=OFS="\t"}{split($3,x,";");class=x[2];print $1,$2,$3,class}' > rfam_anno_class.txt

3. 统计ncRNA数量

先读文件1，将列2和列4也就是名称和class按关系存入字典a。再读文件2，每一行，用自己的列1也就是名称去查询字典，得到class，存入type变量，并计数

awk 'BEGIN{OFS=FS="\t"}ARGIND==1{a[$2]=$4;}ARGIND==2{type=a[$1]; if(type=="") type="Others"; count[type]+=1;}END{for(type in count) print type, count[type];}' rfam_anno_class.txt nep_genome.tblout.final.xls

数量统计结果

4. 统计各类型RNA的总长度

awk 'BEGIN{OFS=FS="\t"}ARGIND==1{a[$2]=$4;}ARGIND==2{type=a[$1]; if(type=="") type="Others"; if($6=="-")len[type]+=$4-$5+1;if($6=="+")len[type]+=$5-$4+1}END{for(type in len) print type, len[type];}' rfam_anno_class.txt nep_genome.tblout.final.xls

长度统计结果