Rfam是用来鉴定non-coding RNAs的数据库，常用于注释新的核酸序列或者基因组序列。
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` #添加环境变量
3 数据库下载
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 准备基因组.fa文件,确定其大小
esl-seqstat smallrna.fa
image.png
5 确定Z值
Z值的计算公式为 Z = total * 2 * CMnumber/1000000
CMnumber确定方法为
less Rfam.cm | grep "ACC" | wc -l
之前做的时候确定Z值都没有乘CMnumber,偶然看到 jianshu.com/p/0bceb4c54474 的帖子,才发现原来是这样算.
官方手册: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 运行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 处理得到的结果
上步得到cmscan和tblout输出文件
tblout文件中
- 表示同一条链上，不存在与此查询序列重叠的序列也在Rfam数据库有匹配,保留。
^ 表示同一条链上，不存在比此查询序列与Rfam数据库匹配更好的序列,保留。
= 表示同一条链上，存在比此查询序列与Rfam数据库匹配更好的序列,删除。
将分隔符设定为制表符
如果是第一行，打印行名称
第3行开始，如果20列不是等号，且不以井号开始，打印其他列
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)全部复制
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 计算小RNA的种类
#写一下思路
#将第四步得到的表格提取第一列计算每个id出现的次数
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语言，使用merge函数匹配id进行注释
#merge函数的用法
c=merge(a,b,by="a/b相同的列名")
write.csv(c,file="res.csv")
参考:
https://www.jianshu.com/p/0bceb4c54474
https://blog.csdn.net/qazplm12_3/article/details/80160292
https://www.jianshu.com/p/89d8b72d9bd5
