Rfam是用來鑒定non-coding RNAs的數(shù)據(jù)庫(kù)，常用于注釋新的核酸序列或者基因組序列厌秒。
infernal說明書:http://eddylab.org/infernal/
1 下載 infernal
wget http://eddylab.org/infernal/infernal-1.1.2-linux-intel-gcc.tar.gz
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ù)庫(kù)下載
wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/14.0/Rfam.cm.gz
cmpress Rfam.cm
wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/14.0/Rfam.clanin
4 準(zhǔn)備基因組.fa文件,確定其大小
esl-seqstat smallrna.fa
image.png
5 確定Z值
Z值的計(jì)算公式為 Z = total * 2 * CMnumber/1000000
CMnumber確定方法為
less Rfam.cm | grep "ACC" | wc -l
之前做的時(shí)候確定Z值都沒有乘CMnumber,偶然看到 jianshu.com/p/0bceb4c54474 的帖子,才發(fā)現(xiàn)原來是這樣算.
官方手冊(cè):For cmscan, Z is defined differently; it is the length of the current query sequence (again, multiplied by 2) in nucleotides multiplied by the number of models in the target CM database.
3 運(yùn)行cmcan
nohup cmscan --cpu 10 -Z (你得到的Z值) --cut_ga --rfam --nohmmonly --tblout smallrna.tblout --fmt 2 --clanin Rfam.clanin Rfam.cm smallrna.fa > smallrna.cmscan
4 處理得到的結(jié)果
上步得到cmscan和tblout輸出文件
tblout文件中
- 表示同一條鏈上限嫌，不存在與此查詢序列重疊的序列也在Rfam數(shù)據(jù)庫(kù)有匹配,保留。
^ 表示同一條鏈上册养，不存在比此查詢序列與Rfam數(shù)據(jù)庫(kù)匹配更好的序列,保留傲绣。
= 表示同一條鏈上掠哥，存在比此查詢序列與Rfam數(shù)據(jù)庫(kù)匹配更好的序列,刪除巩踏。
將分隔符設(shè)定為制表符
如果是第一行，打印行名稱
第3行開始续搀，如果20列不是等號(hào)塞琼，且不以井號(hào)開始，打印其他列
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; }' smallrna.tblout >smallrna.tblout.final.xls
5 下載注釋信息
https://rfam.xfam.org/
(1)search
(2)entry_type 全部勾選 submit
(3)
image.png
(4)全部復(fù)制
vi rfamannotation.txt
i
shift + insert
esc
shift + ;
wq
enter
6 處理注釋信息
less rfamannotation.txt | awk 'BEGIN {FS=OFS="\t"}{split($3,x,";");class=x[2];print $1,$2,$3,class}' > rfam_anno.txt
7 計(jì)算小RNA的種類
#寫一下思路
#將第四步得到的表格提取第一列計(jì)算每個(gè)id出現(xiàn)的次數(shù)
awk '{print $1}' smallrna.tblout.final.xls > id.xls
sort id.xls | uniq -c| sort -nr > id.frequency.txt
#將id.frequency.txt（a） 和 rfam_anno.txt（b） 讀入R語(yǔ)言目代，使用merge函數(shù)匹配id進(jìn)行注釋
#merge函數(shù)的用法
c=merge(a,b,by="a/b相同的列名")
write.csv(c,file="res.csv")
參考:
http://www.reibang.com/p/0bceb4c54474
https://blog.csdn.net/qazplm12_3/article/details/80160292
http://www.reibang.com/p/89d8b72d9bd5
