使用Rfam和Infernal注释基因组的非编码RNA

1. 全过程

官网：https://rfam.org/
官网给的Rfam/Infernal基因组注释教程： https://docs.rfam.org/en/latest/genome-annotation.html
简书文章：https://www.jianshu.com/p/ca3f382c1941

conda create -n rna -y
conda activate rna
conda install -y infernal
cd /mnt/caigui/41_sk_genome_129M/36_infernal

# 我当前下载的是Rfam14.9（2022年11月版本，4108 families，点开https://rfam.org/就可以看到最新版本号了）
wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/CURRENT/Rfam.cm.gz
gunzip Rfam.cm.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/Rfam/CURRENT/Rfam.clanin

# 建立索引文件
cmpress Rfam.cm

# 获得基因组大小
# method1
python ~/scripts/genomesize.py xiaocuiyun_genome_20211217.fa 
# method2
locate esl-seqstat
/home/caigui/miniconda3/envs/SubPhaser/bin/esl-seqstat xiaocuiyun_genome_20211217.fa
# method3
conda install -y easel
esl-seqstat xiaocuiyun_genome_20211217.fa
# 获得基因组大小后，计算Z值 = 132374076*2/1000000=264.748152

# 运行
nohup cmscan -Z 264.748152 --cut_ga --rfam --nohmmonly --tblout sk-genome.tblout --fmt 2 --clanin Rfam.clanin Rfam.cm xiaocuiyun_genome_20211217.fa > sk-genome.cmscan &

# 关于这个Z值如何计算，有两种不同的说法。Rfam的方法就是以上方法。但是infernal的教程中还要乘以一个模型数量(模型数量获得方法：cat Rfam.cm | grep 'NAME'|sort|wc -l)。
# 为了解决这个问题，我写了封邮件，如下
# rfam-help@ebi.ac.uk
# Hello,
# I hope this message finds you well. I have a question regarding the calculation of Z-values on the webpage "https://docs.rfam.org/en/latest/genome-annotation.html". On page 39 of "http://eddylab.org/infernal/Userguide.pdf", it is mentioned that when calculating Z-values using cmscan, it should be multiplied by "the number of models in the target CM database". I am a bit confused about this and would appreciate it if you could clarify how the Z-value is actually calculated.
# Thank you very much for your help. Have a great day!

# 获得两个结果文件
$ ls sk-genome.*
sk-genome.cmscan sk-genome.tblout

# 从tblout文中删除得分较低的重叠
grep -v " = " sk-genome.tblout > sk-genome.deoverlapped.tblout

# Rfam/Infernal可以预测所有类型的RNA
# 以下是另外几种用于搜索特定类型RNA的工具
# tRNAscan-SE for tRNAs
# RNAMMER for rRNA
# snoscan for snoRNAs
# SRPscan for SRP RNA

# 提取需要保存的信息
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; }' sk-genome.deoverlapped.tblout > sk-genome.deoverlapped.clean.tblout

# 获得注释，https://rfam.org/search，Entry type，勾选全部类型，Submit，划到页面最下面，保存unformatted list格式
vim rfam.annot.txt #粘贴以上部分

# 提取注释文件第三列信息
cat rfam.annot.txt | awk 'BEGIN {FS=OFS="\t"}{split($3,x,";");class=x[2];print $1,$2,$3,class}' > rfam.annot.modif.txt

# 统计小RNA种类
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.annot.modif.txt sk-genome.deoverlapped.clean.tblout
#  riboswitch    2
#  ribozyme    3
#  tRNA    599
#  Others    10
#  miRNA    14
#  rRNA    364
#  sRNA    1
#  snRNA    336

# 统计小RNA长度
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.annot.modif.txt sk-genome.deoverlapped.clean.tblout
#  riboswitch    204
#  Others    1694
#  ribozyme    712
#  tRNA    44286
#  miRNA    1571
#  rRNA    215788
#  sRNA    379
#  snRNA    46451

# 附上以上命令的注释，出自chatGPT
awk 'BEGIN{OFS=FS="\t"} # 设置输入输出字段分隔符为tab
ARGIND==1{a[$2]=$4;} # 处理第一个文件，将第二列作为关键字，第四列作为值存储在一个名为a的数组中
ARGIND==2{ # 处理第二个文件
    type=a[$1]; # 根据第一列的值从a数组中获取类型
    if(type=="") type="Others"; # 如果类型为空，则将其设置为Others
    if($6=="-")len[type]+=$4-$5+1; # 如果第6列的值为负号，则使用第四列的值减去第五列的值加1来更新该类型的长度
    if($6=="+")len[type]+=$5-$4+1 # 如果第6列的值为正号，则使用第五列的值减去第四列的值加1来更新该类型的长度
}
END{ # 处理完两个文件后，进行最终的输出
    for(type in len) # 遍历len数组中的每个类型
        print type, len[type]; # 输出类型和对应的长度
}' rfam.annot.modif.txt sk-genome.deoverlapped.clean.tblout

2. 修改Z值重新预测一下，比较一下什么结果

万一人家不回邮件呢，咱再试试。理论上来说，这个Z值好像只影响E值的计算结果吧？

# cd /mnt/caigui/41_sk_genome_129M/36_infernal/Z_big
cat Rfam.cm | grep 'NAME'|sort|wc -l
# 8216
# 132374076*2*8216/1000000 = 2175171
nohup cmscan -Z 2175171 --cut_ga --rfam --nohmmonly --tblout sk-genome.tblout --fmt 2 --clanin Rfam.clanin Rfam.cm xiaocuiyun_genome_20211217.fa > sk-genome.cmscan &

# 经过比较，确实只有E值改变了，而检测出来的基因种类数量都完全没有变。
# 再比较一下，基因E值，在Z值乘了8216倍后，变大了很多。但这个玩意是越小表明显著性越高。
# 随便算了算，Z值变大的倍数也基本上就是变大了8216倍左右。

3. 试一试对拟南芥进行注释

因为注释出来的ncRNA数量好像不是很多，所以对拟南芥试试看，使用big Z。

nohup cmscan -Z 1966380 --cut_ga --rfam --nohmmonly --tblout at-genome.tblout --fmt 2 --clanin Rfam.clanin Rfam.cm Athaliana_447_TAIR10.fa > at-genome.cmscan &

# 统计小RNA种类
#  riboswitch    1
#  ribozyme    2
#  tRNA    658
#  Others    42
#  miRNA    398
#  antisense    1
#  rRNA    561
#  snRNA    654
#  sRNA    9

# 统计小RNA长度
#  riboswitch    128
#  ribozyme    500
#  tRNA    49842
#  Others    6132
#  miRNA    61678
#  antisense    187
#  rRNA    85627
#  snRNA    71522
#  sRNA    1788

4. 比较一下SK和AT注释出来的关于miRNA的内容

# 小翠云miRNA个数
grep -E "MIR|mir-" sk-genome.deoverlapped.clean.tblout | wc -l
# 14个
# 小翠云miRNA种类
grep -E "MIR|mir-"  sk-genome.deoverlapped.clean.tblout | awk '{print $1}' | sort | uniq| wc -l
# 5种

# 拟南芥miRNA个数
grep -E "MIR|mir-"  at-genome.deoverlapped.clean.tblout| wc -l
# 398个
# 拟南芥miRNA种类
grep -E "MIR|mir-"  at-genome.deoverlapped.clean.tblout | awk '{print $1}' | sort | uniq| wc -l
# 74种

在Rfam中大概看了下，小翠云有的这几个miRNA确实在很多植物中都是特别保守的，但是种类和个数确实是很少啊。

5. 邮件后续

简单来说，Rfam的维护者Blake表示也对这个问题感到迷惑，并帮我咨询了Infernal的作者Eric，Eric表示就按Rfam文档中的来就好，不用乘以模型数量。
附上Eric的详细解释如下：

这么说，这是全世界第一份关于这个问题的正确解释了？