用cmscan挖掘ncRNA信息
更新日志：
2018年12月27日: 对部分内容进行了修正，添加了新版本Rfam的下载链接。对Z参数的计算进行了修改。
2019年2月8日: 增加了一个使用示例
0.准备工作：
- 获取Rfam种子
- 获取Rfam的claninfo
- 软件安装
- 待处理物种的基因组文件
新建一个专门用于处理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
源码安装
官网提供可下载的源码和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/106
因此还要计算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 的数量，除以106 -
--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
参考：
