Rfam是用来鉴定非编码RNA的数据库，常用于注释新的核酸序列或基因组序列。
最新版本的Rfam14.4，已经包含了3941个家族。
Rfam官方网址：https://rfam.xfam.org/#searchTypeBlock
Rfam用户手册：http://eddylab.org/infernal/Userguide.pdf
Rfam的安装和使用流程：
1 下载infernal
wget http://eddylab.org/infernal/infernal-1.1.2-linux-intel-gcc.tar.gz
infernal官方网址：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` #添加环境变量
3 数据库下载
wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/14.4/Rfam.cm.gz
cmpress Rfam.cm（建立索引）
我自己的应用位置：
./infernal_bin/bin/cmpress ./Rfam.cm
索引生成的文件：
wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/14.0/Rfam.clanin
4 准备基因组的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
输出的结果如下：
其中输出的一列结果中，有一列Total # residues :500736760是我们所要的。有的博客中写道“考虑到基因组为双链和下一步用到的参数的单位为million，所以使用公式
Z = total * 2 * CMmumber/106
在这个过程中还要计算CM database中的模型的数量
在Rfam14.4中使用：
less Rfam.cm | grep 'NAME'|sort|wc -l
我得到的结果是7874
在这一步我的是单链的，所以我没有乘2；
在这个过程中，会出现一个问题就是你使用esl-seqstat命令的时候会出现如下情况：
esl-seqstat命令是hmmer的一个插件，如果没法全局调用则建议直接locate esl-seqstat查找绝对路径。
我查询到的路径是在：
/home/student/grh/Software/hmmer-3.2.1/easel/miniapps/esl-seqstat
当然，我参看博客有说，是因为我的interproscan没装好导致没法直接使用。
5 运行程序
nohup cmscan --cpu 10 -Z (你得到的Z值) --cut_ga --rfam --nohmmonly --tblout smallrna.tblout --fmt 2 --clanin Rfam.clanin Rfam.cm smallrna.fa > smallrna.cmscan
我自己运行的例子：
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 &
因为比较耗时，所以就使用nohup命令放在后台运行着；
各个参数说明：
-Z：查询序列的大小，以M为单位。由esl-seqstat算出或自己写程序计算，记得乘以2，乘以CM model 的数量，除以1000000
--cut_ga: 输出不小于Rfam GA阈值的结果。这是Rfam认证RNA家族的阈值，不低于这个阈值的序列得分被认为是真同源序列。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 结果处理
上一步得到cmscan和tblout输出文件
在filename.tblout文件中，有一栏是olp，表示查询序列的重叠信息:
. 表示同一条链上，不存在与此查询序列重叠的序列也在Rfam数据库有匹配，这是需要保留的查询序列。
^表示同一条链上，不存在比此查询序列与Rfam数据库匹配更好的序列，也需要保留。
= 表示同一条链上，存在比此类查询序列与Rfam数据库匹配更好的序列，因此应将搜索到=的行给去除掉。
grep -v '=' my-genome.tblout >my-genome.deoverlapped.tblout
将文件转成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数据库进行下载；
（1）进入网页后选择search；
（2）entry_type 全部勾选后选择Submit；
进入页面示例：
在文件的底部有一句话描述了进行复制，粘贴的方法；选择unformatted text，复制粘贴下来。
显示如下示例：
复制下来后，保存为一个txt文件；
8 处理注释信息
把文件中的第三列拆开，取出其类型，存储为HJ_rfam_anno.txt
运行命令：
less rfamannotation.txt | awk 'BEGIN {FS=OFS="\t"}{split($3,x,";");class=x[2];print $1,$2,$3,class}' > HJ_rfam_anno.txt
9 计算小RNA的种类
运行命令如下:
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
这个时候就能统计出小RNA的类型：
统计各个小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];}' /home/Lab/huangjin/16_miRNA_software/rfam/HJ_rfam_anno.txt ./CK-1.tblout.final.xls
参考博客：
https://www.jianshu.com/p/844066501c32
https://www.jianshu.com/p/0bceb4c54474
https://blog.csdn.net/qazplm12_3/article/details/80160292
https://www.jianshu.com/p/89d8b72d9bd5
对于该软件的安装和使用，参考以上博客进行汇总并自己进行实践，流程如上！
Hjin 2021.1.8
