本地使用Rfam 12.0+

2018-04-29 本文已影响207人生信宝典

欢迎关注天下博客：http://blog.genesino.com/2017/06/Rfam/
Jump to...

Rfam是用来鉴定non-coding RNAs的数据库，常用于注释新的核酸序列或者基因组序列。

最新版本的Rfam (12.2)，包含2588个RNA家族，其在线网站提供了便捷的查询使用功能，http://rfam.xfam.org/，尤其是对小批量数据。

对已经注释好的物种，建议在ENSEMBLE或UCSC直接下载官方的注释文件，直接下载GTF或使用bioMart或TableBrowser都可。具体可在本博客或微信公众号后台回复关键字获取使用方法。

最后确定要在本地使用了，如果是用之前的版本，请在线搜索，有几篇教程可用。如果有新版本，请参考此教程。

下载infernal软件

官网下载：infernal软件本来应该在http://eddylab.org/infernal/下载，但这个网站最近打不开了。

GitHub下载

# 如果不会或不能使用git，拷贝链接，在GitHub下载压缩文件，再按顺序解压即可
# https://github.com/EddyRivasLab 
git clone https://github.com/EddyRivasLab/infernal.git infernal
cd infernal
git clone https://github.com/EddyRivasLab/easel.git
git clone https://github.com/EddyRivasLab/hmmer.git

#如果当前目录有aclocal.m4，则不需要执行
#ln -s `pwd`/easel/aclocal.m4 `pwd`/

ln -s `pwd`/easel/aclocal.m4 hmmer

# 如果没有autoconf，找管理员配置，或查看软件安装
autoconf
(cd easel; autoconf)
(cd hmmer; autoconf)

./configure  --prefix=`pwd`/../infernal_bin
make 
make install
cd easel; make install

cd ../../infernal_bin/bin
ls

# 就可以看到安装的程序了，加入环境变量即可，本站搜索环境变量查看具体配置方法
# 或微信公众号 生信宝典 后台回复 环境变量 查看 
export PATH=${PATH}:`pwd`

下载数据库

wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/12.2/Rfam.cm.gz
gunzip Rfam.cm.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/12.2/Rfam12.2.claninfo
# 使用 Infernal的程序cmpress索引Rfam.cm 
# 大约需要一分钟
cmporess Rfam.cm

# 输出为
Working...    done.
Pressed and indexed 2588 CMs and p7 HMM filters (2588 names and 2588 accessions).
Covariance models and p7 filters pressed into binary file:  Rfam.cm.i1m
SSI index for binary covariance model file:                 Rfam.cm.i1i
Optimized p7 filter profiles (MSV part)  pressed into:      Rfam.cm.i1f
Optimized p7 filter profiles (remainder) pressed into:      Rfam.cm.i1p

确定待查询的核苷酸序列或基因组的大小，作为后续命令的参数

# 大约1分钟
esl-seqstat my-genome.fa

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

使用cmscan程序注释基因组的ncRNAs

# Rfam12.2.claninfo 为下载的claninfo文件，需提供所在路径
# Rfam.cm 下载的cm文件
# my-genome.fa 待查询序列
# my-genome.cmscan 输出结果
# my-genome.tblout 有一个输出结果
# 对500M大小的输入序列，48线程，需要7个小时，最好放入后台
cmscan -Z `esl-seqstat my-genome.fa | awk '{if($0~/^Total/) print int($4/2000000);}''` --cut_ga --rfam --nohmmonly --tblout my-genome.tblout --fmt 2 --clanin Rfam12.2.claninfo Rfam.cm my-genome.fa > my-genome.cmscan

参数解释:

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

结果解释

my-genome.cmscan

cmscan命令的标准输出，使用>重定向。

结果的第一部分是运行的命令的参数记录

第二部分是查询的FASTA文件中每个序列的top hits, 根据E-value排序。

Query:       scaffold12  [L=2429009]
Hit scores:
 rank     E-value  score  bias  modelname   start     end   mdl trunc   gc  description
 ----   --------- ------ -----  --------- ------- -------   --- ----- ----  -----------
  (1) !   5.4e-20   86.6   0.0  mir-166    961624  961752 +  cm    no 0.43  -
  (2) !   9.6e-14   70.9   0.0  tRNA      2369877 2369805 -  cm    no 0.59  -
  (3) !   3.7e-13   68.8   0.0  tRNA      1344161 1344232 +  cm    no 0.53  -
  (4) !   2.3e-12   65.9   0.0  tRNA       973203  973274 +  cm    no 0.54  -
  (5) !   5.8e-12   64.5   0.0  tRNA      1241524 1241595 +  cm    no 0.50  -

E-value: 统计显著性，依赖于查询数据库的大小。
score: Log-odds得分，独立于查询序列数据库的大小。在使用了--cut_ga后所有输出的结果都是高于Rfam GA得分的。
modelname: Rfam家族或模型的名字。
start, stop: 查询序列匹配的区域。后面跟着的是链的信息，对于+，起始位置小于终止位置；对于-，其实位置大于终止位置。

互有重叠的查询区域可能会匹配到不同的Rfam家族或模型。推荐保留具有最低的E-value，最高的bit scaore的部分。

结果的下一部分是比对结果。具体可查看文后的参考网址链接的内容。

my-genome.tblout

表格输出包含了cmscan标准输出的大部分内容，并且便于对结果的进一步处理。

#idx target name          accession query name           accession clan name mdl mdl from   mdl to seq from   seq to strand trunc pass   gc  bias  score   E-value inc olp anyidx afrct1 afrct2 winidx wfrct1 wfrct2 description of target
#--- -------------------- --------- -------------------- --------- --------- --- -------- -------- -------- -------- ------ ----- ---- ---- ----- ------ --------- --- --- ------ ------ ------ ------ ------ ------ ---------------------
1    tRNA                 RF00005   scaffold20           -         CL00001    cm        1       71  1399503  1399576      +    no    1 0.57   0.0   58.8   2.4e-10  !   *       -      -      -      -      -      - -
2    tRNA                 RF00005   scaffold20           -         CL00001    cm        1       71   186338   186267      -    no    1 0.54   0.0   55.4     2e-09  !   *       -      -      -      -      -      - -
1    mir-166              RF00075   scaffold12           -         -          cm        1      126   961624   961752      +    no    1 0.43   0.0   86.6   5.4e-20  !   *       -      -      -      -      -      - -
2    tRNA                 RF00005   scaffold12           -         CL00001    cm        1       71  2369877  2369805      -    no    1 0.59   0.0   70.9   9.6e-14  !   *       -      -      -      -      -      - -
3    tRNA                 RF00005   scaffold12           -         CL00001    cm        1       71  1344161  1344232      +    no    1 0.53   0.0   68.8   3.7e-13  !   *       -      -      -      -      -      - -
4    tRNA                 RF00005   scaffold12           -         CL00001    cm        1       71   973203   973274      +    no    1 0.54   0.0   65.9   2.3e-12  !   *       -      -      -      -      -      - -
5    tRNA                 RF00005   scaffold12           -         CL00001    cm        1       71  1241524  1241595      +    no    1 0.50   0.0   64.5   5.8e-12  !   *       -      -      -      -      -      - -
1    Plant_SRP            RF01855   scaffold15           -         CL00003    cm        1      305  1511627  1511325      -    no    1 0.62   0.7  249.5   1.1e-70  !   ^       -      -      -      -      -      - -

每一行有27列，比较关键的是`target name`, `accession`, `query name`, `seq from`, `seq to`, `strand`, `E-value`, `score`。

`olp`列表示查询序列的重叠信息，`*`表示同一条链上，不存在与此查询序列重叠的序列也在Rfam数据库有匹配，这是需要保留的查询序列。`^`表示同一条链上，不存在比此查询序列与Rfam数据库匹配更好的序列，也需要保留。`=`表示同一条链上，存在比此查询序列与Rfam数据库匹配更好的序列，应忽略。

可通过下面的命令获取最终结果。

```bash
grep -v '=' my-genome.tblout >my-genome.deoverlapped.tblout

结果解析

首先转换下结果格式，提取必须得列和非重叠区域或重叠区域中得分高的区域。

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

统计预测出的ncRNA的类型。

首先下载Rfam家族的注释，点击http://rfam.xfam.org/search#tabview=tab5，选择所有复选框，提交，把得到的表格拷贝下来，整理成TAB键分割的格式。并吧第三列拆开，取出类型, 存储为 Rfam_anno.txt。

Accession   ID  Type    Description class
RF00001 5S_rRNA Gene; rRNA  5S ribosomal RNA    rRNA
RF00002 5_8S_rRNA   Gene; rRNA  5.8S ribosomal RNA  rRNA
RF00003 U1  Gene; snRNA; splicing   U1 spliceosomal RNA snRNA

 awk 'BEGIN{OFS=FS="\t"}ARGIND==1{a[$2]=$5;}ARGIND==2{type=a[$1]; if(type=="") type="Others"; count[type]+=1;}END{for(type in count) print type, count[type];}' Rfam_anno.txt my-genome.tblout.final.xls

最终输出

rRNA    61
snRNA   397
sRNA    1
Others  55
antisense   1
tRNA    427
miRNA   95
ribozyme    1
riboswitch  1

Reference

生信宝典，学的更多

Rfam 12.0+本地使用（最新版教程）

RFAMCHENTONG
版权声明：本文为博主原创文章，转载请注明出处。

alipay.png

WeChatPay.png

</footer>