Rfam数据库的Small nuclear RNA分类统计

常见的Non-coding RNA基因主要分四类:rRNA基因、tRNA基因、snRNA基因和miRNA基因。tRNA基因一般使用tRNAScan-SE软件进行分析;rRNA基因一般使用RNAmmer软件或BLAST软件进行分析;对于snRNA基因和miRNA基因,则使用Infernal软件将基因组序列比对到Rfam数据库进行分析得到。但Rfam数据库分析能得到很多Non-conding RNA信息,提取其中的snRNA信息则有点麻烦。具体分析方法如下:

1. snRNA简介

snRNA一类较小的存在细胞核中的RNA分子,通常~150个核苷酸长度。可以分成四类:

1. Sm-class snRNA: 由RNA聚合酶II转录,包含U1, U2, U4, U4atac, U5, U7, U11,和U12几类。该类snRNA能从核孔转移到细胞质中,最后形成splicesome,用于对mRNA前体进行加工,去除introns。

2. Lsm-class snRNA:和第一类snRNA类似,但由RNA聚合酶III转录,不会离开细胞核。

3. C/D box snoRNA: 存在核仁中的RNA分子,包含两个保守的序列motifs,C (RUGAUGA) and D (CUGA),通常能对rRNA进行甲基化加工。

4. H/ACA box snoRNA:存在核仁中的RNA分子,包含两个保守的序列motifs,H box (consensus ANANNA) and the ACA box (ACA),通常能对tRNA进行假尿苷化加工。

2. snoRNA分类和RF编号的获取

相比于snRNA,snoRNA家族庞大。若需要统计snoRNA信息,首先需要了解那些Rfam中的编号属于那种类型的snoRNA。

在Rfam数据库文件Rfam.cm中以关键词Small nucleolar RNA搜索,从而确定属于snoRNA的Rfam编号。

perl -e '$/ = "\n//"; while (<>) { if (m/Small nucleolar RNA/i) { my $name = $1 if m/NAME\s+(\S+)/; my $acc = $1 if m/ACC\s+(\S+)/; print "$acc\t$name\n"; } }' Rfam.cm > snoRNA_1.list

然后,根据Rfam编号在Rfam网站中鉴定所属的snoRNA分类,C/D box或 H/ACA box。此步骤需要编写如下程序snoRNA_type.pl,来自动下载相应Rfam编号的网页信息,并解析其所属分类。

#!/usr/bin/perl
use strict;

my $info = `curl http://rfam.xfam.org/family/$ARGV[0]`;

if ($info =~ m/Gene; snRNA; snoRNA; ([^\s;]+)/) {
    print "$ARGV[0]\t$1\n";
}
else {
    print "$ARGV[0]\tother\n";
}

批量化运行解析程序,并合并结果得到snoRNA的Rfam编号和分类信息。

for i in `cut -f 1 snoRNA_1.list`
do
    echo "perl snoRNA_type.pl $i > $i.type"
done > command.snoRNA_type.list
ParaFly -c command.snoRNA_type.list -CPU 50

cat *.type > snoRNA_2.list
rm *.type
cut -f 2 snoRNA_1.list > aa
paste snoRNA_2.list aa > snoRNA.list
rm aa snoRNA_1.list snoRNA_2.list snoRNA_type.pl

3. 编写程序统计snRNA信息

程序名称:snRNA_stats_from_Rfam.pl

#!/usr/bin/perl
use strict;
use Getopt::Long;

my $usage = <<USAGE;
Usage:
    $0 [options] rfam_out.tab snoRNA_type.txt

    --out-prefix <string>    default: snRNA_out
    设置输出文件前缀。默认输出文件:
    (1)snRNA_out.stats 对各类snRNA的统计信息。
    (2)snRNA_out.txt 对Rfam数据库进行Infernal分析中属于snRNA的表格结果。
    (3)snRNA_out.gff3 将第2个结果转换为GFF3格式的结果。

    程序通过读取Rfam的infernal分析结果,统计出small nuclear RNA的信息。snRNA_out.stat结果文件包含多列:
    第1列:snRNA分类
    第2列:snRNA名称
    第2列:RF编号
    第3列:基因家族数量
    第4列:序列平均长度
    第5列:序列总长度

    small nuclear RNA分四大类:
    Sm-class snRNA: (U1, U2, U4, U4atac, U5, U6, U7, U11, U12)
    Lsm-class snRNA: (U6, U6atac)
    C/D box snoRNA
    H/ACA box snoRNA

    snoRNA是最复杂的一类,在Rfam v14.1中包含565个家族,其中388个属于C/D box,177个属于H/ACA box。通过读取snoRNA_type.txt文件中RF编号和snoRNA分类对应的信息,来统计两类snoRNA的数量。

USAGE
if (@ARGV==0){die $usage}

my ($out_prefix);
GetOptions(
    "out-prefix:s" => \$out_prefix,
);
$out_prefix ||= "snRNA_out";

# Sm-class或Lsm-class类型的snRNA在Rfam数据库中名称和RF编号获取:
# perl -e '$/ = "\n//"; while (<>) { if (m/small nuclear RNA/i or m/spliceosomal RNA/i) { my $name = $1 if m/NAME\s+(\S+)/; my $acc = $1 if m/ACC\s+(\S+)/; print "$acc\t$name\n"; } }' Rfam.cm | uniq

my %snRNA_RF2Name = (
    "RF00003" => "U1",
    "RF00004" => "U2",
    "RF00007" => "U12",
    "RF00015" => "U4",
    "RF00020" => "U5",
    "RF00026" => "U6",
    "RF00066" => "U7",
    "RF00488" => "U1_yeast",
    "RF00548" => "U11",
    "RF00618" => "U4atac",
    "RF00619" => "U6atac",
    "RF01844" => "SmY",
    "RF02491" => "Gl_U1",
    "RF02492" => "Gl_U2",
    "RF02493" => "Gl_U4",
    "RF02494" => "Gl_U6"
);

my %snRNA_RF2Class = (
    "RF00003" => "Sm-class snRNA",
    "RF00004" => "Sm-class snRNA",
    "RF00007" => "Sm-class snRNA",
    "RF00015" => "Sm-class snRNA",
    "RF00020" => "Sm-class snRNA",
    "RF00026" => "Lsm-class snRNA",
    "RF00066" => "Sm-class snRNA",
    "RF00488" => "Sm-class snRNA",
    "RF00548" => "Sm-class snRNA",
    "RF00618" => "Sm-class snRNA",
    "RF00619" => "Lsm-class snRNA",
    "RF01844" => "Sm-class snRNA",
    "RF02491" => "Sm-class snRNA",
    "RF02492" => "Sm-class snRNA",
    "RF02493" => "Sm-class snRNA",
    "RF02494" => "Sm-class snRNA"
);

open IN, $ARGV[1] or die "Can not open file $ARGV[1], $!\n";
while (<IN>) {
    next if m/^#/;
    next if m/^\s*$/;
    chomp;
    @_ = split /\t/;
    $snRNA_RF2Name{$_[0]} = $_[2];

    if ($_[1] eq "CD-box") {
        $snRNA_RF2Class{$_[0]} = 'C/D box snoRNA';
    }
    elsif ($_[1] eq "HACA-box") {
        $snRNA_RF2Class{$_[0]} = 'H/ACA box snoRNA';
    }
    else {
        $snRNA_RF2Class{$_[0]} = 'other';
    }
}
close IN;

my %class;
foreach (keys %snRNA_RF2Class) {
    $class{$snRNA_RF2Class{$_}}{$_} = 1;
}

open OUT1, ">", "$out_prefix.txt" or die "Can not create file $out_prefix.txt, $!\n";
open OUT2, ">", "$out_prefix.gff3" or die "Can not create file $out_prefix.gff3, $!\n";
open OUT3, ">", "$out_prefix.stats" or die "Can not create file $out_prefix.stats, $!\n";
my (%stats, %gff3, $total_snRNA_number);
open IN, $ARGV[0] or die "Can not open file $ARGV[0], $!\n";
while (<IN>) {
    if (m/^#/) { print OUT1; next; }
    next if m/^\s*$/;

    @_ = split /\s+/;
    if (exists $snRNA_RF2Name{$_[3]}) {
        $total_snRNA_number ++;
        print OUT1;
        $gff3{$_} = 1;

        $stats{$_[3]}{"num"} ++;
        $stats{$_[3]}{"total_length"} += (abs($_[8] - $_[7]) + 1);
    }
}
close IN;

my @class = ("Sm-class snRNA", "Lsm-class snRNA", 'C/D box snoRNA', 'H/ACA box snoRNA');
push @class, "other" if exists $class{"other"};
foreach my $class (@class) {
    my @rf = keys %{$class{$class}};
    my (%out1 , %out2);
    foreach (keys %{$class{$class}}) {
        my $num = $stats{$_}{"num"};
        next unless $num;
        my $total_length = $stats{$_}{"total_length"};
        my $avg = int(($total_length / $num) * 100 + 0.5) / 100;
        my $out = "$class\t$snRNA_RF2Name{$_}\t$_\t$num\t$avg\t$total_length";
        $out1{$out} = $num;
        $out2{$out} = $total_length;
    }
    foreach (sort {$out1{$b} <=> $out1{$a} or $out2{$b} <=> $out2{$a}} keys %out1) {
        print OUT3 "$_\n";
    }
}

my (%sortGFF1, %sortGFF2, %sortGFF3, %sortGFF4, %gff3_name, %gff3_RF, $number);
foreach (keys %gff3) {
    @_ = split /\s+/, $_;
    my ($seqID, $start, $end, $score, $strand, $name, $rf_code) = ($_[0], $_[7], $_[8], $_[14], $_[9], $_[2], $_[3]);
    if ($start > $end) {
        my $tmp_vaule = $start;
        $start = $end;
        $end = $tmp_vaule;
    }

    my $out = "$seqID\t.\tsnRNA\t$start\t$end\t$score\t$strand\t\.";
    $sortGFF1{$out} = $seqID;
    $sortGFF2{$out} = $start;
    $sortGFF3{$out} = $end;
    $sortGFF4{$out} = $strand;
    $gff3_name{$out} = $name;
    $gff3_RF{$out} = $rf_code;
}

foreach (sort {$sortGFF1{$a} cmp $sortGFF1{$b} or $sortGFF2{$a} <=> $sortGFF2{$b} or $sortGFF3{$a} <=> $sortGFF3{$b} or $sortGFF4{$a} cmp $sortGFF4{$b}} keys %sortGFF1) {
    $number ++;
    my $id = '0' x (length($total_snRNA_number) - length($number)) . $number;
    print OUT2 "$_\tID=snRNA_$id;Name=$gff3_name{$_};RF=$gff3_RF{$_}\n";
}
close OUT1; close OUT2; close OUT3;

4. 编写程序统计miRNA信息

程序名称:miRNA_stats_from_Rfam.pl

#!/usr/bin/perl
use strict;
use Getopt::Long;

my $usage = <<USAGE;
Usage:
    $0 [options] rfam_out.tab Rfam.cm

    --out-prefix <string>    default: miRNA_out
    设置输出文件前缀。默认输出文件:
    (1)miRNA_out.stats 对各类microRNA的统计信息。
    (2)miRNA_out.txt 对Rfam数据库进行Infernal分析中属于microRNA的表格结果.
    (3)miRNA_out.gff3 将第2个结果转换为GFF3格式的结果。

    程序通过读取Rfam数据库信息,使用关键词microRNA搜索,找到属于miRNA的编号,再分析infernal结果文件,统计miRNA信息。

USAGE
if (@ARGV==0){die $usage}

my ($out_prefix);
GetOptions(
    "out-prefix:s" => \$out_prefix,
);
$out_prefix ||= "miRNA_out";

open IN, $ARGV[1] or die "Can not open file $ARGV[1], $!\n";
my %miRNA_RF2Name;
$/ = "\n//";
while (<IN>) {
    if (m/microRNA/i) {
        my $name = $1 if m/NAME\s+(\S+)/;
        my $acc = $1 if m/ACC\s+(\S+)/;
        $miRNA_RF2Name{$acc} = $name;
    }
}
close IN;
$/ = "\n";

open OUT1, ">", "$out_prefix.txt" or die "Can not create file $out_prefix.txt, $!\n";
open OUT2, ">", "$out_prefix.gff3" or die "Can not create file $out_prefix.gff3, $!\n";
open OUT3, ">", "$out_prefix.stats" or die "Can not create file $out_prefix.stats, $!\n";
my (%stats, %gff3, $total_miRNA_number);
open IN, $ARGV[0] or die "Can not open file $ARGV[0], $!\n";
while (<IN>) {
    if (m/^#/) { print OUT1; next; }
    next if m/^\s*$/;

    @_ = split /\s+/;
    if (exists $miRNA_RF2Name{$_[3]}) {
        $total_miRNA_number ++;
        print OUT1;
        $gff3{$_} = 1;

        $stats{$_[3]}{"num"} ++;
        $stats{$_[3]}{"total_length"} += (abs($_[8] - $_[7]) + 1);
    }
}
close IN;

foreach (sort {$stats{$b}{"num"} <=> $stats{$a}{"num"}} keys %stats) {
    my $num = $stats{$_}{"num"};
    next unless $num;
    my $total_length = $stats{$_}{"total_length"};
    my $avg = int(($total_length / $num) * 100 + 0.5) / 100;
    print OUT3 "$miRNA_RF2Name{$_}\t$_\t$num\t$avg\t$total_length\n";
}

my (%sortGFF1, %sortGFF2, %sortGFF3, %sortGFF4, %gff3_name, %gff3_RF, $number);
foreach (keys %gff3) {
    @_ = split /\s+/, $_;
    my ($seqID, $start, $end, $score, $strand, $name, $rf_code) = ($_[0], $_[7], $_[8], $_[14], $_[9], $_[2], $_[3]);
    if ($start > $end) {
        my $tmp_vaule = $start;
        $start = $end;
        $end = $tmp_vaule;
    }

    my $out = "$seqID\t.\tmiRNA\t$start\t$end\t$score\t$strand\t\.";
    $sortGFF1{$out} = $seqID;
    $sortGFF2{$out} = $start;
    $sortGFF3{$out} = $end;
    $sortGFF4{$out} = $strand;
    $gff3_name{$out} = $name;
    $gff3_RF{$out} = $rf_code;
}

foreach (sort {$sortGFF1{$a} cmp $sortGFF1{$b} or $sortGFF2{$a} <=> $sortGFF2{$b} or $sortGFF3{$a} <=> $sortGFF3{$b} or $sortGFF4{$a} cmp $sortGFF4{$b}} keys %sortGFF1) {
    $number ++;
    my $id = '0' x (length($total_miRNA_number) - length($number)) . $number;
    print OUT2 "$_\tID=miRNA_$id;Name=$gff3_name{$_};RF=$gff3_RF{$_}\n";
}
close OUT1; close OUT2; close OUT3;

发表评论

电子邮件地址不会被公开。 必填项已用*标注

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据