一：下載安裝

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 #添加環(huán)境變量

二：下載數(shù)據(jù)庫

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

三：確定待查詢的核苷酸序列或基因組的大小矫户，作為后續(xù)命令的參數(shù)

esl-seqstat my-genome.fa

Z值

其輸出結(jié)果中有一行，Total # of residues: 218502668是我們需要的银舱。考慮到基因組為雙鏈和下一步用到的參數(shù)的單位為Million跛梗。
計算公式為：Z = total * 2 * CMmumber/106
因此還要計算CM database中的模型的數(shù)量（Z=2636016）
在Rfam14.0版本中纵朋，使用

五：運(yùn)行程序

# Rfam14.1.clanin 為下載的claninfo文件，需提供所在路徑
# Rfam.cm 下載的cm文件
#nep_genome.fa 待查詢序列
# nep.cmscan 輸出結(jié)果
# nep.tblout 表格形式輸出結(jié)果
# 對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

建議加線程宪祥，別把服務(wù)器跑滿了Ｄ粜健！

六：輸出結(jié)果

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列蝗羊，比較關(guān)鍵的是`target name`, `accession`, `query name`, `seq from`, `seq to`, `strand`, `E-value`, `score`藏澳。
`olp`列表示查詢序列的重疊信息，`*`表示同一條鏈上耀找，不存在與此查詢序列重疊的序列也在Rfam數(shù)據(jù)庫有匹配翔悠，這是需要保留的查詢序列。`^`表示同一條鏈上野芒，不存在比此查詢序列與Rfam數(shù)據(jù)庫匹配更好的序列蓄愁，也需要保留。`=`表示同一條鏈上狞悲，存在比此查詢序列與Rfam數(shù)據(jù)庫匹配更好的序列撮抓，應(yīng)忽略。
## 可通過下面的命令獲取最終結(jié)果摇锋。
```bash
grep -v '=' nep.tblout >nep .deoverlapped.tblout

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

結(jié)果解析

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

1.格式轉(zhuǎn)換

得到tblout.final，轉(zhuǎn)換下結(jié)果格式，提取必須得列和非重疊區(qū)域或重疊區(qū)域中得分高的區(qū)域琅攘。
=表示同一條鏈上浪腐，存在比此查詢序列與Rfam數(shù)據(jù)庫匹配更好的序列，應(yīng)忽略。
去掉#開頭的注釋行提取所需要的列。初始將分隔符設(shè)定為制表符。
如果是第一行媳纬，打印行名稱第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家族的注釋，點(diǎn)擊Rfam網(wǎng)址七芭，選擇所有復(fù)選框素挽，提交，把得到的表格拷貝下來狸驳，整理成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. 統(tǒng)計ncRNA數(shù)量

先讀文件1耙箍，將列2和列4也就是名稱和class按關(guān)系存入字典a撰糠。再讀文件2，每一行辩昆，用自己的列1也就是名稱去查詢字典阅酪，得到class，存入type變量汁针，并計數(shù)

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

數(shù)量統(tǒng)計結(jié)果

4. 統(tǒng)計各類型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

長度統(tǒng)計結(jié)果

infernal預(yù)測ncRNA

infernal預(yù)測ncRNA

一：下載安裝

二：下載數(shù)據(jù)庫

三：確定待查詢的核苷酸序列或基因組的大小矫户，作為后續(xù)命令的參數(shù)

五：運(yùn)行程序

六：輸出結(jié)果

結(jié)果解析

1.格式轉(zhuǎn)換

2. 下載Rfam注釋

3. 統(tǒng)計ncRNA數(shù)量

4. 統(tǒng)計各類型RNA的總長度

infernal預(yù)測ncRNA

一：下載安裝

二：下載數(shù)據(jù)庫

三： 確定待查詢的核苷酸序列或基因組的大小矫户，作為后續(xù)命令的參數(shù)

五：運(yùn)行程序

六：輸出結(jié)果

結(jié)果解析

1.格式轉(zhuǎn)換

2. 下載Rfam注釋

3. 統(tǒng)計ncRNA數(shù)量

4. 統(tǒng)計各類型RNA的總長度

三：確定待查詢的核苷酸序列或基因組的大小矫户，作為后續(xù)命令的參數(shù)