Rfam统计各种RNA
1. 得到Rfam最终结果
1. 代码参考
perl /bin/rfam_scan.pl -blastdb /biosoft/rfam11.0/Rfam.fasta /biosoft/rfam11.0/Rfam.cm {genomeSeq} -o {fileName(genomeSeq)}.rfam.gff3
sum_ncRNA: python {sumNCRNA} -r {rRNAgff2} -t {tRNAout} -s {rfamGff3} -o {fileName(tRNAout)}.ncRNA.result.csv
2. 运行
2.1 确定待查询的genome大小
esl-seqstat my-genome.fa
输出结果
Total # of residues是我们需要的。考虑到基因组为双链和下一步用到的参数的单位为Million，我们使用公式total * 2 / 1000000计算得出结果为11，作为下一步参数-Z的值.
2.2 cmscan程序注释基因组的ncRNAs
cmscan -Z 11 --cut_ga --rfam --nohmmonly --tblout my-genome.tblout --fmt 2 --clanin ../../biosoft/Rfam/Rfam12.2.claninfo ../../biosoft/Rfam/Rfam.cm ../assembly.fasta > my-genome.cmscan
输出结果
得到cmscan和tblout输出文件
2.3 提取需要的信息，得到tblout.final
转换下结果格式，提取必须得列和非重叠区域或重叠区域中得分高的区域。
-
=表示同一条链上，存在比此查询序列与Rfam数据库匹配更好的序列，应忽略。 - 去掉#开头的注释行
- 提取所需要的列。
- 初始将分隔符设定为制表符
- 如果是第一行，打印行名称
- 第3行开始，如果20列不是等号，且不以井号开始，打印2,3,4等列
- 输出
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
2.根据注释信息分组统计
Rfam 注释信息
选中所有复选框，submit后，在文件底部选择unformatted text，复制粘贴下来。（注意看help信息）
cat rfam_anno.txt | awk 'BEGIN {FS=OFS="\t"}{split($3,x,";");class=x[2];print $1,$2,$3,class}' > rfam_anno_class.txt
将第三列的第二个分号信息提取出来，这就是各ncRNA class。
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
先读文件1，将列2和列4也就是名称和class按关系存入字典a。
读文件2，每一行，用自己的列1也就是名称去查询字典，得到class，存入type变量，并计数。
image.png
即rfam-anno-class.txt的第二列名称对应最后一列的RNA类别。xls按第一列查询。
个人认为RF号更精确一些。换RF号查询结果基本一致，ribo多1个，leader多1个。
得到结果:总和178
riboswitch 10
tRNA 77
ribozyme 1
Others 27
leader 2
thermoregulator 1
antisense 23
rRNA 22
sRNA 15
awk 'BEGIN{OFS=FS="\t"}ARGIND==1{a[$2]=$4;}ARGIND==2{type=a[$1]; if(type=="") type="Others"; if($6=="-")len[type]+=$4-$5+1;if($6=="+")len[type]+=$5-$4+1}END{for(type in len) print type, len[type];}' rfam_anno_class.txt my-genome.tblout.final.xls
得到各类型RNA的总长度
riboswitch 1708
tRNA 5782
Others 3446
ribozyme 360
leader 219
thermoregulator 69
antisense 2024
rRNA 31938
sRNA 2926
3.参考
基本来自本文 link
