用cmscan挖掘ncRNA信息

0.准备工作：

获取Rfam种子
获取Rfam的claninfo
软件安装
待处理物种的基因组文件

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

获取种子和chanin文件

下载Rfam种子：

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

下载clanin文件：

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

软件安装

用conda安装infernal软件

conda install infernal

源码安装

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的值.

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 &

因为比较耗时所以建议使用nohup命令来跑
-Z：查询序列的大小，以M为单位。由esl-seqstat算出或自己写程序计算，记得乘以2，除以10^6
--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

参考：

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