更新日志：
2019年11月4日: 了解了一下预测的结果都是些啥
2018年12月27日: 对部分内容进行了修正，添加了新版本Rfam的下载链接。对Z参数的计算进行了修改。
2019年2月8日: 增加了一个使用示例

0.准备工作：

1. 获取Rfam种子
2. 获取Rfam的claninfo
3. 软件安装
4. 待处理物种的基因组文件

新建一个专门用于处理RNA的文件夹mkdir Cmscan

获取种子和chanin文件

下载Rfam种子：

 axel -q ftp://ftp.ebi.ac.uk/pub/databases/Rfam/13.0/Rfam.cm.gz

2018年12月27日 update：Rfam已更新至14.0.请用下面的链接：
wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/14.0/Rfam.cm.gz

下载clanin文件：

wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/13.0/Rfam.clanin

2018年12月27日 update：Rfam已更新至14.0.请用下面的链接：
wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/14.0/Rfam.clanin

软件安装

用conda安装infernal软件

conda create -n miRNA
conda activate miRNA
conda install -c bioconda infernal=1.1.2
conda install -y hmmer=3.2.1

源码安装

http://eddylab.org/infernal/

官网提供可下载的源码和User's Guide，相当人性化了。并且mac系统也可以使用brew安装infernal。

准备基因组文件

将待处理的基因组文件软链接到Cmscan文件夹下：

ln -s /path/to/file.fa

1.建库

cmpress Rfam.cm

2.确定基因组大小

esl-seqstat my-genome.fa

其输出结果中有一行，类似于Total # of residues: 3000000是我们需要的。考虑到基因组为双链和下一步用到的参数的单位为Million，我们使用公式3000000 * 2 / 1000000计算得出结果为6，作为下一步参数-Z的值.

== 2018年12月27日 update： ==
根据本篇评论中@代码0019旁友的提醒，回去查了一下tutorial，确实有提到cmsearch和cmscan的Z值在计算时候是有所不同的。原文如下:

For cmscan, Z is the length of the current query sequence multiplied by 2 (because both strands of the sequence are searched) and multiplied again by the number of CMs in the target CM database

正确的计算公式为：Z = total * 2 * CMmumber/10⁶
因此还要计算CM database中的模型的数量
在Rfam14.0版本中，使用

less Rfam.cm | grep 'NAME'|sort|wc -l

得到结果为5582

tips:esl-seqstat命令是hmmer的一个插件，如果没法全局调用则建议直接locate esl-seqstat查找绝对路径。在我的服务器上它在的位置是/media/newdisk/interproscan-5.28-67.0/bin/hmmer/hmmer3/3.1b1/easel/miniapps/esl-seqstat

当然，有可能是因为我的interproscan没装好导致没法直接使用。。

3.运行程序（举个栗子）

nohup cmscan -Z 208 --cut_ga --rfam --nohmmonly --tblout kfl.tblout --fmt 2 --clanin /media/newdisk/Cmscan/Rfam.clanin  Rfam.cm /media/newdisk/lunzao/KFL/120824_klebsormidium_Scaffolds_v1.0.fna > kfl.cmscan &

nohup cmscan -Z 3503 --cut_ga --rfam --nohmmonly --tblout chara.tblout --fmt 2 --clanin /media/newdisk/Cmscan/Rfam.clanin  Rfam.cm /media/newdisk/lunzao/Chara/dailydata/chara_genome.fasta > chara.cmscan &

2019年2月8日更新：很多的参数都需要自己摸索一下，不能随意添加上去。添加了--cpu的参数，这样就不会跑满整个服务器不用担心被kill掉了~

cmscan -Z 3825680 --rfam --tblout papaya.tblout --fmt 2 --cpu 16  --clanin ./Rfam.clanin  Rfam.cm ../data/papaya.genome.fa

• 因为比较耗时所以建议使用nohup命令来跑

• -Z：查询序列的大小，以M为单位。由esl-seqstat算出或自己写程序计算，记得乘以2，乘以CM model 的数量，除以10⁶

• --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.

4.结果处理

在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; }' my-genome.tblout >my-genome.tblout.final.xls

tip：若不设置--cpu的话会默认使用全部线程。

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 my-genome.tblout.final.xls

2019年11月4日更新

结果解读

一直想知道预测的结果都包含了哪些东西，那些不常见的都是些啥玩意儿~今天来理一下
类似于rRNA和tRNA之类的就不管了，太常见。整理一下不常见的。

U开头的系列

tmRNA？

tmRNA RF00023 transfer-messenger RNA 传递信使RNA? 神奇

SSU_rRNA系列

还有一个LSU_rRNA系列，large subunit ribosomal RNA，大亚基核糖体RNA。哈哈原来如此。

S开头的一系列玩意儿

参考：

https://www.jianshu.com/p/89d8b72d9bd5

http://rfam.xfam.org/search/type

https://blog.csdn.net/qazplm12_3/article/details/73380016

https://mp.weixin.qq.com/s/5OIRHA22ZLr5Z8bEhDiBqg

http://blog.genesino.com/2017/06/Rfam/

http://rfam.readthedocs.io/en/latest/genome-annotation.html