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

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

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

最后確定要在本地使用了，如果是用之前的版本，請(qǐng)?jiān)诰€搜索，有幾篇教程可用。如果有新版本，請(qǐng)參考此教程。

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

GitHub下載

如果不會(huì)或不能使用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

如果當(dāng)前目錄有aclocal.m4，則不需要執(zhí)行

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

ln -s pwd/easel/aclocal.m4 hmmer

如果沒(méi)有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)境變量查看具體配置方法

或微信公眾號(hào) 生信寶典 后臺(tái)回復(fù) 環(huán)境變量 查看

export PATH=${PATH}:pwd
下載數(shù)據(jù)庫(kù)
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
其輸出結(jié)果中有一行，類似于Total # of residues: 3000000是我們需要的。考慮到基因組為雙鏈和下一步用到的參數(shù)的單位為Million，我們使用公式3000000 * 2 / 1000000計(jì)算得出結(jié)果為6，作為下一步參數(shù)-Z的值.

使用cmscan程序注釋基因組的ncRNAs

Rfam12.2.claninfo 為下載的claninfo文件，需提供所在路徑

Rfam.cm 下載的cm文件

my-genome.fa 待查詢序列

my-genome.cmscan 輸出結(jié)果

my-genome.tblout 有一個(gè)輸出結(jié)果

對(duì)500M大小的輸入序列，48線程，需要7個(gè)小時(shí)，最好放入后臺(tái)

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算出或自己寫程序計(jì)算，記得乘以2，除以10^6.

–cut_ga: 輸出不小于Rfam GA閾值的結(jié)果。這是Rfam認(rèn)證RNA家族的閾值，不低于這個(gè)閾值的序列得分被認(rèn)為是真同源序列。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.

結(jié)果解釋
my-genome.cmscan

cmscan命令的標(biāo)準(zhǔn)輸出，使用>重定向。

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

結(jié)果的下一部分是比對(duì)結(jié)果。具體可查看文后的參考網(wǎng)址鏈接的內(nèi)容。
my-genome.tblout

表格輸出包含了cmscan標(biāo)準(zhǔn)輸出的大部分內(nèi)容，并且便于對(duì)結(jié)果的進(jìn)一步處理。

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列，比較關(guān)鍵的是target name, accession, query name, seq from, seq to, strand, E-value, score。

olp列表示查詢序列的重疊信息，*表示同一條鏈上，不存在與此查詢序列重疊的序列也在Rfam數(shù)據(jù)庫(kù)有匹配，這是需要保留的查詢序列。^表示同一條鏈上，不存在比此查詢序列與Rfam數(shù)據(jù)庫(kù)匹配更好的序列，也需要保留。=表示同一條鏈上，存在比此查詢序列與Rfam數(shù)據(jù)庫(kù)匹配更好的序列，應(yīng)忽略。

可通過(guò)下面的命令獲取最終結(jié)果。

“`bash
grep -v ‘=’ my-genome.tblout >my-genome.deoverlapped.tblout
結(jié)果解析
首先轉(zhuǎn)換下結(jié)果格式，提取必須得列和非重疊區(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)計(jì)預(yù)測(cè)出的ncRNA的類型。

首先下載Rfam家族的注釋，點(diǎn)擊http://rfam.xfam.org/search#tabview=tab5，選擇所有復(fù)選框，提交，把得到的表格拷貝下來(lái)，整理成TAB鍵分割的格式。并吧第三列拆開，取出類型, 存儲(chǔ)為 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/
生信寶典，學(xué)的更多
Rfam 12.0+本地使用 （最新版教程）
RFAM
CHENTONG
版權(quán)聲明：本文為博主原創(chuàng)文章，轉(zhuǎn)載請(qǐng)注明出處。
alipay.png WeChatPay.png

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

總結(jié)

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

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