基因组注释从0开始(3):ncRNA注释

06 Jan 2025 3733字 13分
CC BY 4.0 （除特别声明或转载文章外）
如果这篇博客帮助到你，可以请我喝一杯咖啡~

1. ncRNA

非编码RNA(Non-coding RNA, ncRNA) 包括rRNA，tRNA，snRNA，snoRNA 和microRNA 等不编码蛋白质的RNA，它们转录后直接在RNA 水平上就能行使各自的生物学功能，并不需要翻译成蛋白质。

2. 注释软件

  • 非编码RNA种类繁多，且结构特征各不相同，所以开发出了许多注释特定某一类RNA的软件，比如tRNAScan-SE预测tRNA，rnammer预测rRNA，snoScan 搜索带C/D盒的snoRNAs，SnoGps 搜索带H/ACA盒的snoRNAs，mirScan搜索microRNA等。
  • Sanger实验室开发了Infernal软件，建立了1600多个RNA家族，并对每个家族建立了一致性二级结构和协方差模型，形成了Rfam数据库。采用Rfam数据库中的每个RNA的协方差模型，结合Infernal软件可以预测出已有RNA家族的新成员，只是特异性比较差。

如果不是专门研究ncRNA，可以用Infernal注释所有ncRNA。如果需要更精细的注释，则可以选择特定软件注释特定RNA。

这篇博客是介绍用Infernal程序与Rfam数据库一起用来注释与数据库中已知ncRNA同源的序列（这里用来注释完整的基因组）。注释结果包括tRNA，rRNA，snRNA，snoRNA和miRNA等。

3. tRNAscan-SE

tRNAscan-SE的安装需要依赖Infernal软件，因此可以用conda直接安装tRNAscan-SE顺带解决依赖问题：

conda install -c bioconda trnascan-se

写一个脚本运行tRNAscan-SE：

mkdir 3.ncRNA
cd 3.ncRNA
mkdir 1.tRANscan-se
cd 1.tRANscan-se

conda activate transcan-se
tRNAscan-SE --thread 50 -E -I -m tRNA_luobuma.stats -o tRNA_luobuma.out -f tRNA_luobuma.ss ../../sample.fa

参数解释：

-E 搜寻真核生物tRNA

-I 使用Infernal软件进行搜索

-m 保存结果统计文件

-o 输出tRNA预测结果

-f 保存tRNA的二级结构

这里只是粗略统计一下tRNA数量，可根据stat查看：

4. Infernal + Rfam

Infernal（INFERence of RNA ALignment）是Sanger实验室开发的ncRNA预测软件，他们建立了1600多个RNA家族，每个家族建立了一致性二级结构和协方差模型，也就是Rfam数据库。总体的注释思路是基因组与 Rfam数据库进行比对，Rfam是一个RNA分类数据库，其比对方法是调用软件Infernal中的程序cmscan，将提交的序列在Rfam.cm数据库中进行检索，从而得到其比对的结果。

cmscan（search sequence(s) against a covariance model database, 针对协方差模型数据库的序列搜索），主要参考官方手册Userguide.pdf (eddylab.org)中的Searching the Rfam CM database with a query sequence步骤。

前面已经安装了Infernal，这里需要再下载一个Rfam数据库。

mkdir 2.infernal
cd 2.infernal

# 下载Rfam数据库，注意两个文件版本必须一致
wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/CURRENT/Rfam.cm.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/CURRENT/Rfam.clanin

# 解压建库
gunzip Rfam.cm.gz
cmpress Rfam.cm

4.1 运行

4.1.1 序列索引

推荐的参数：

nohup cmscan -Z 1024 --cut_ga --rfam --nohmmonly --fmt 2 --tblout sample.tblout -o sample.result --clanin Rfam.clanin Rfam.cm ../..sample.fa &
  • -Z：根据基因组大小来定，基因组大小的2倍，Mb单位，选一个整数。比如512Mb的基因组，-Z 1024。
  • --cut_ga --rfam --nohmmonly --fmt 2：推荐使用
  • –tblout sample.tblout：指定table格式输出文件
  • -o sample.result：指定比对结果输出文件
  • –clanin Rfam.clanin：指定clanin文件
  • Rfam.cm genome.fa：指定数据库Rfam.cm和基因组genome.fa

note：-o sample.result要放在Rfam.cm genome.fa前面，否则报错。

此步骤耗时参考：585Mb基因组，默认线程，耗时4h左右。

4.1.2结果文件

  • sample.result：比对结果
  • sample.tblout：table格式结果

5. 整理结果

5.1. 将注释结果整理成gff3文件

gff3文件可用于提交注释到数据库。

用perl脚本infernal-tblout2gff.pl实现，脚本来自https://www.cnblogs.com/jessepeng/p/15392809.html。

perl infernal-tblout2gff.pl --cmscan --fmt2 sample.tblout >sample.infernal.ncRNA.gff3

5.2. 统计各类ncRNA总数

  1. 整理注释结果文件sample.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; }' sample.tblout >sample.tblout.xls

  2. 下载rfam注释

  • 在rfam官网，选择【SEARCH】-【Entry type】
  • 然后选中所有的Entry types（包括Gene，Intron，Cis-regulatory element），点击【Submit】，会列出所有RNA family的注释信息。
  • 手动选择所有注释信息，复制，粘贴到新建的空白文本文件rfam.txt并保存。
  • 把rfam.txt传输到服务器，最好用dos2unix rfam.txt转换文件格式为unix版本。
  • 拆分第三列cat rfam.txt | awk 'BEGIN {FS=OFS="\t"}{split($3,x,";");class=x[2];print $1,$2,$3,$4,class}' > rfam_anno.txt

rfam注释文件rfam_anno.txt包含了所有rfam的类型type和功能描述description信息。

  1. 统计ncRNA注释结果 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 sample.tblout.xls >sample.ncRNA.statistic

sample.ncRNA.statistic输出示例：

 riboswitch     1
 ribozyme       3
 tRNA   600
Others  30
 miRNA  267
 antisense      1
 rRNA   524
 sRNA   1
 snRNA  6730
  1. 统计细分分类
  • 也可以根据细分分类分别统计，细分分类参考rfam官网，【SEARCH】-【Entry type】。
  • 可参考的统计值包括每个细分ncRNA的数量(copy)，平均长度(average length)，总长(total length)，总长占基因组的比例（Percent of the genome）等
  • 统计后整理成发表文章用的表格。

比如：snRNA包括snoRNA和splicing，snoRNA包括CD-box，HACA-box和scaRNA。下面用统计CD-box这个细分分类的ncRNA举例。

  • 提取CD-box的Accession(RF00000格式)：grep "CD-box" rfam_anno.txt |cut -f1 >cdbox.tem
  • 提取注释到的CD-box信息：grep -f cdbox.tem mc.tblout.xls >cdbox.txt
  • cdbox.txt的行数就是CD-box的数量；利用第四五列的位置信息，即可统计平均长度(average length)，总长(total length)，注意正反链。