歡迎關注天下博客：http://blog.genesino.com/2017/06/Rfam/
Jump to…
下載infernal軟件
下載數(shù)據(jù)庫
確定待查詢的核苷酸序列或基因組的大小，作為后續(xù)命令的參數(shù)
使用cmscan程序注釋基因組的ncRNAs
結果解釋
my-genome.cmscan
my-genome.tblout
結果解析
Reference
生信寶典，學的更多
Rfam是用來鑒定non-coding RNAs的數(shù)據(jù)庫，常用于注釋新的核酸序列或者基因組序列。

最新版本的Rfam (12.2)，包含2588個RNA家族，其在線網(wǎng)站提供了便捷的查詢使用功能，http://rfam.xfam.org/，尤其是對小批量數(shù)據(jù)。

對已經(jīng)注釋好的物種，建議在ENSEMBLE或UCSC直接下載官方的注釋文件，直接下載GTF或使用bioMart或TableBrowser都可。具體可在本博客或微信公眾號后臺回復關鍵字獲取使用方法。

最后確定要在本地使用了，如果是用之前的版本，請在線搜索，有幾篇教程可用。如果有新版本，請參考此教程。

下載infernal軟件
官網(wǎng)下載：infernal軟件本來應該在http://eddylab.org/infernal/下載，但這個網(wǎng)站最近打不開了。

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，則不需要執(zhí)行

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

就可以看到安裝的程序了，加入環(huán)境變量即可，本站搜索環(huán)境變量查看具體配置方法

或微信公眾號 生信寶典 后臺回復 環(huán)境變量 查看

export PATH=${PATH}:pwd
下載數(shù)據(jù)庫
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
確定待查詢的核苷酸序列或基因組的大小，作為后續(xù)命令的參數(shù)

大約1分鐘

esl-seqstat my-genome.fa
其輸出結果中有一行，類似于Total # of residues: 3000000是我們需要的。考慮到基因組為雙鏈和下一步用到的參數(shù)的單位為Million，我們使用公式3000000 * 2 / 1000000計算得出結果為6，作為下一步參數(shù)-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
參數(shù)解釋:

-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命令的標準輸出，使用>重定向。

結果的第一部分是運行的命令的參數(shù)記錄
第二部分是查詢的FASTA文件中每個序列的top hits, 根據(jù)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: 統(tǒng)計顯著性，依賴于查詢數(shù)據(jù)庫的大小。
score: Log-odds得分，獨立于查詢序列數(shù)據(jù)庫的大小。在使用了–cut_ga后所有輸出的結果都是高于Rfam GA得分的。
modelname: Rfam家族或模型的名字。
start, stop: 查詢序列匹配的區(qū)域。后面跟著的是鏈的信息，對于+，起始位置小于終止位置；對于-，其實位置大于終止位置。
互有重疊的查詢區(qū)域可能會匹配到不同的Rfam家族或模型。推薦保留具有最低的E-value，最高的bit scaore的部分。

結果的下一部分是比對結果。具體可查看文后的參考網(wǎng)址鏈接的內容。
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數(shù)據(jù)庫有匹配，這是需要保留的查詢序列。^表示同一條鏈上，不存在比此查詢序列與Rfam數(shù)據(jù)庫匹配更好的序列，也需要保留。=表示同一條鏈上，存在比此查詢序列與Rfam數(shù)據(jù)庫匹配更好的序列，應忽略。

可通過下面的命令獲取最終結果。

“`bash
grep -v ‘=’ my-genome.tblout >my-genome.deoverlapped.tblout
結果解析
首先轉換下結果格式，提取必須得列和非重疊區(qū)域或重疊區(qū)域中得分高的區(qū)域。

awk ‘BEGIN{OFS=”\t”;}{if(FNR==1) print “target_name\taccession\tquery_name\tquery_start\tquery_end\tstrand\tscore\tEvalue”; if(FNR>2 && 20!="="?&&20!="="?&&0!~/^#/) print 2,2,3,4,4,10,11,11,12,17,17,18; }’ my-genome.tblout >my-genome.tblout.final.xls
統(tǒng)計預測出的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]=2]=5;}ARGIND==2{type=a[$1]; if(type==”“) type=”O(jiān)thers”; 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
http://rfam.readthedocs.io/en/latest/genome-annotation.html
http://rfam.xfam.org/
生信寶典，學的更多
Rfam 12.0+本地使用 （最新版教程）
RFAM
CHENTONG
版權聲明：本文為博主原創(chuàng)文章，轉載請注明出處。
alipay.png WeChatPay.png

創(chuàng)作挑戰(zhàn)賽新人創(chuàng)作獎勵來咯，堅持創(chuàng)作打卡瓜分現(xiàn)金大獎

總結

以上是生活随笔為你收集整理的本地使用Rfam 12.0+的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內容還不錯，歡迎將生活随笔推薦給好友。