一：下载安装

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

长度统计结果

infernal预测ncRNA

infernal预测ncRNA

一：下载安装

二：下载数据库

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

五：运行程序

六：输出结果

结果解析

1.格式转换

2. 下载Rfam注释

3. 统计ncRNA数量

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

infernal预测ncRNA

一：下载安装

二：下载数据库

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

五：运行程序

六：输出结果

结果解析

1.格式转换

2. 下载Rfam注释

3. 统计ncRNA数量

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

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