Rfam是用來鑒定非編碼RNA的數(shù)據(jù)庫，常用于注釋新的核酸序列或基因組序列。
最新版本的Rfam14.4，已經(jīng)包含了3941個(gè)家族。
Rfam官方網(wǎng)址：https://rfam.xfam.org/#searchTypeBlock
Rfam用戶手冊(cè)：http://eddylab.org/infernal/Userguide.pdf
Rfam的安裝和使用流程：
1 下載infernal
wget http://eddylab.org/infernal/infernal-1.1.2-linux-intel-gcc.tar.gz
infernal官方網(wǎng)址：http://eddylab.org/infernal/
2 安裝infernal
tar xf infernal-1.1.2-linux-intel-gcc.tar.gz
cd infernal-1.1.2-linux-intel-gcc.tar.gz
./configure? --prefix=`pwd`/../infernal_bin
make
make install
cd easel
make install
cd ../../infernal_bin/bin
ls #查看可以用的軟件
export PATH=${PATH}:`pwd` #添加環(huán)境變量
3 數(shù)據(jù)庫下載
wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/14.4/Rfam.cm.gz
cmpress Rfam.cm（建立索引）
我自己的應(yīng)用位置：
./infernal_bin/bin/cmpress ./Rfam.cm
索引生成的文件：
wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/14.0/Rfam.clanin
4 準(zhǔn)備基因組的fasta文件
可以將待處理的fasta文件軟鏈接到cmscan文件夾下：
例如：ln -s /path/to/file.fa
以我的為例：
ln -s /home/Lab/huangjin/16_miRNA_software/miRDP_2/CSS_miReads/CK-1.1.fa
確定查詢fasta文件的大小
esl-seqstat my-genome.fa
以我的為例：
esl-seqstat ./CK-1.1.fa
輸出的結(jié)果如下：
其中輸出的一列結(jié)果中，有一列Total # residues :500736760是我們所要的。有的博客中寫道“考慮到基因組為雙鏈和下一步用到的參數(shù)的單位為million，所以使用公式
Z = total * 2 * CMmumber/106
在這個(gè)過程中還要計(jì)算CM database中的模型的數(shù)量
在Rfam14.4中使用：
less Rfam.cm | grep 'NAME'|sort|wc -l
我得到的結(jié)果是7874
在這一步我的是單鏈的，所以我沒有乘2；
在這個(gè)過程中，會(huì)出現(xiàn)一個(gè)問題就是你使用esl-seqstat命令的時(shí)候會(huì)出現(xiàn)如下情況：
esl-seqstat命令是hmmer的一個(gè)插件，如果沒法全局調(diào)用則建議直接locate esl-seqstat查找絕對(duì)路徑。
我查詢到的路徑是在：
/home/student/grh/Software/hmmer-3.2.1/easel/miniapps/esl-seqstat
當(dāng)然，我參看博客有說，是因?yàn)槲业膇nterproscan沒裝好導(dǎo)致沒法直接使用。
5 運(yùn)行程序
nohup cmscan --cpu 10 -Z (你得到的Z值) --cut_ga --rfam --nohmmonly --tblout smallrna.tblout --fmt 2 --clanin Rfam.clanin Rfam.cm smallrna.fa > smallrna.cmscan
我自己運(yùn)行的例子：
nohup cmscan -Z 3942801 --rfam --tblout CK_1.tblout --fmt 2 --cpu 20 --clanin /home/Lab/huangjin/16_miRNA_software/rfam/Rfam.clanin /home/Lab/huangjin/16_miRNA_software/rfam/Rfam.cm /home/Lab/huangjin/16_miRNA_software/miRDP_2/CSS_miReads/CK-1.1.fa? > CK-1.cmscan &
因?yàn)楸容^耗時(shí)，所以就使用nohup命令放在后臺(tái)運(yùn)行著；
各個(gè)參數(shù)說明：
-Z：查詢序列的大小，以M為單位。由esl-seqstat算出或自己寫程序計(jì)算，記得乘以2，乘以CM model 的數(shù)量，除以1000000
--cut_ga: 輸出不小于Rfam GA閾值的結(jié)果。這是Rfam認(rèn)證RNA家族的閾值，不低于這個(gè)閾值的序列得分被認(rèn)為是真同源序列。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.
6 結(jié)果處理
上一步得到cmscan和tblout輸出文件
在filename.tblout文件中，有一欄是olp，表示查詢序列的重疊信息:
. 表示同一條鏈上，不存在與此查詢序列重疊的序列也在Rfam數(shù)據(jù)庫有匹配，這是需要保留的查詢序列。
^表示同一條鏈上，不存在比此查詢序列與Rfam數(shù)據(jù)庫匹配更好的序列，也需要保留。
= 表示同一條鏈上，存在比此類查詢序列與Rfam數(shù)據(jù)庫匹配更好的序列，因此應(yīng)將搜索到=的行給去除掉。
grep -v '=' my-genome.tblout >my-genome.deoverlapped.tblout
將文件轉(zhuǎn)成excel形式，只保留所需要的信息。
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; }' CK_1.tblout >CK-1.tblout.final.xls
7 下載注釋信息
注釋信息主要通過Rfam數(shù)據(jù)庫進(jìn)行下載；
（1）進(jìn)入網(wǎng)頁后選擇search；
（2）entry_type 全部勾選后選擇Submit；
進(jìn)入頁面示例：
在文件的底部有一句話描述了進(jìn)行復(fù)制，粘貼的方法；選擇unformatted text，復(fù)制粘貼下來。
顯示如下示例：
復(fù)制下來后，保存為一個(gè)txt文件；
8 處理注釋信息
把文件中的第三列拆開，取出其類型，存儲(chǔ)為HJ_rfam_anno.txt
運(yùn)行命令：
less rfamannotation.txt | awk 'BEGIN {FS=OFS="\t"}{split($3,x,";");class=x[2];print $1,$2,$3,class}' > HJ_rfam_anno.txt
9 計(jì)算小RNA的種類
運(yùn)行命令如下:
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];}' /home/Lab/huangjin/16_miRNA_software/rfam/HJ_rfam_anno.txt ./CK-1.tblout.final.xls
這個(gè)時(shí)候就能統(tǒng)計(jì)出小RNA的類型：
統(tǒng)計(jì)各個(gè)小RNA 的總長(zhǎng)度
運(yùn)行命令如下：
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];}' /home/Lab/huangjin/16_miRNA_software/rfam/HJ_rfam_anno.txt ./CK-1.tblout.final.xls
參考博客：
http://m.itdecent.cn/p/844066501c32
http://m.itdecent.cn/p/0bceb4c54474
https://blog.csdn.net/qazplm12_3/article/details/80160292
http://m.itdecent.cn/p/89d8b72d9bd5
對(duì)于該軟件的安裝和使用，參考以上博客進(jìn)行匯總并自己進(jìn)行實(shí)踐，流程如上！
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? Hjin 2021.1.8
