用cmscan挖掘ncRNA信息(2019-11-4更新)

更新日志:
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/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

2019年11月4日更新

结果解读

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

U开头的系列

剪接体
U1 RF00003 U1 spliceosomal RNA
U2 RF00004 U2 spliceosomal RNA
U3 RF00012 U3 spliceosomal RNA
U4 RF00015 U4 spliceosomal RNA
U5 RF00020 U5 spliceosomal RNA
U6 RF00026 U6 spliceosomal RNA

tmRNA?

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

SSU_rRNA系列

SSU_rRNA_microsporidia RF02542 微孢子虫小亚基核糖体RNA
SSU_rRNA_eukarya RF01960 真核……
SSU_rRNA_bacteria RF00177 细菌……
SSU_rRNA_archaea RF01959 古菌……

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

S开头的一系列玩意儿

Name ID Description Translation
sRNA162 RF02792 Archaeal Small RNA 162 古细菌小RNA
SorX RF02784 Singlet oxygen resistance RNA X 单线态抗氧性RNA X
snR77 RF01181 Small nucleolar RNA snR77 小核仁RNA
snoZ267 RF00344 Small nucleolar RNA Z267 也是一个小核仁RNA
SmY RF01844 SmY spliceosomal RNA 也是一种剪接体
SL1 RF00198 SL1 RNA 这玩意儿就叫做SL1 RNA,家族里就只有俩,一个SL1,一个SL2
rli47 RF01478 Listeria sRNA rli47 李斯特菌的小RNA

参考:

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

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 206,839评论 6 482
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 88,543评论 2 382
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 153,116评论 0 344
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 55,371评论 1 279
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 64,384评论 5 374
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,111评论 1 285
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,416评论 3 400
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,053评论 0 259
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 43,558评论 1 300
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,007评论 2 325
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,117评论 1 334
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,756评论 4 324
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,324评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,315评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,539评论 1 262
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,578评论 2 355
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,877评论 2 345

推荐阅读更多精彩内容

  • 今天因为一大早要去为人大代表换届选举服务,所以取消了早上L5的课程。(在寒冷的天气下,户外工作三个小时,也没能扫出...
    尾八女士阅读 585评论 0 51
  • 前言:被病魔折磨近一个月的舅妈今天终于去了。弥留状态的舅妈,直直地盯着我妈__她的大姑姐与她的丈夫__我...
    厌南瓜症患者阅读 420评论 0 2
  • (1)A:记得上次心理课时,考试布置了一个作业,要求我们做一个最后的宝典设计的提纲。我们没有及时交,导致平时成绩扣...
    灿生_c9b2阅读 275评论 0 1
  • 威腾智库:个人品牌营销学 梁红兵威腾,2018年7月 一、个人品牌的核心系统 1 品类: 最细分的族群类别、阶层类...
    鼎优威腾智库阅读 543评论 0 0