NCBI gff/gff3 轉(zhuǎn)換為 gtf 格式

適用于 NCBI 數(shù)據(jù)庫中 gff 格式文件轉(zhuǎn)換為 gtf 格式徘熔,使用 perl 代碼實現(xiàn):

#!/usr/bin/perl -w
use strict;
use Getopt::Long;
use Data::Dumper;

my ($help, $infile, $outdir, $gff2gtf, $gene2tr);
GetOptions(
    "infile=s"  => \$infile,
    "outdir:s"  => \$outdir,
    "gene2rt!"  => \$gene2tr,
    "gff2gtf!"  => \$gff2gtf,
    "help|?"    => \$help
);
die &usage if (!defined $infile || defined $help);
$outdir ||= ".";system("mkdir -p $outdir");


if (defined $gff2gtf){
    my ($gene, $gene_info, $mRNA, $exon, $cds) = &gff2gtf($infile);

    foreach my $g (sort keys %$gene){
        my $out_g = $gene_info->{$g};
        my $gene_id = $gene_info->{$g}->[4];
        my $gene_type = $gene_info->{$g}->[5];
        
        print "$out_g->[0]\tGeneBang\tgene\t$out_g->[1]\t$out_g->[2]\t.\t$out_g->[3]\t.\tgene_id \"$gene_id\"; gi \"$g\"; gene_biotype \"$gene_type\";\n";

        foreach my $m (sort keys %{$mRNA->{$g}}){
            if (exists $mRNA->{$g}->{$m}->{'mRNA'}){
                my $out_m = $mRNA->{$g}->{$m}->{'mRNA'};
                my $transcript_id = $out_m->[4];
                print "$out_m->[0]\tGeneBang\ttranscript\t$out_m->[1]\t$out_m->[2]\t.\t$out_m->[3]\t.\tgene_id \"$gene_id\"; transcript_id \"$transcript_id\"; gi \"$g\"; transcript_biotype \"$gene_type\";\n";

                if (exists $exon->{$m}){
                    foreach my $e (sort keys %{$exon->{$m}}){
                        my $out_e = $exon->{$m}->{$e};
                        print "$out_e->[0]\tGeneBang\texon\t$out_e->[1]\t$out_e->[2]\t.\t$out_e->[3]\t.\tgene_id \"$gene_id\"; transcript_id \"$transcript_id\"; gi \"$g\"; transcript_biotype \"$gene_type\"; exon_id \"$e\"; \n";
                    }
                }
                
                if (exists $cds->{$m}){
                    foreach my $c (sort keys %{$cds->{$m}}){
                        my @out_pos = @{$cds->{$m}->{$c}->{pos}};
                        my $protein_id = $cds->{$m}->{$c}->{pep};
                        foreach my $pos (@out_pos){
                            print "$out_m->[0]\tGeneBang\tCDS\t$pos\t.\t$out_m->[3]\t.\tgene_id \"$gene_id\"; transcript_id \"$transcript_id\"; gi \"$g\"; transcript_biotype \"$gene_type\"; cds_id \"$c\"; protein_id \"$protein_id\"; \n";
                        }
                    }
                }


            }elsif(exists $mRNA->{$g}->{$m}->{'ncRNA'} ){
                my $out_n = $mRNA->{$g}->{$m}->{'ncRNA'};
                my $bio_type = $out_n->[4];
                print "$out_n->[0]\tGeneBang\ttranscript\t$out_n->[1]\t$out_n->[2]\t.\t$out_n->[3]\t.\tgene_id \"$gene_id\"; transcript_id \"$m\"; gi \"$g\"; transcript_biotype \"$bio_type\";\n";

                if (exists $exon->{$m}){
                    foreach my $e (sort keys %{$exon->{$m}}){
                        my $out_e = $exon->{$m}->{$e};
                        print "$out_e->[0]\tGeneBang\texon\t$out_e->[1]\t$out_e->[2]\t.\t$out_e->[3]\t.\tgene_id \"$gene_id\"; transcript_id \"$m\"; gi \"$g\"; transcript_biotype \"$bio_type\"; exon_id \"$e\"; \n";
                    }
                }
            }else{
                next;
            }
        }
    }
}


if (defined $gene2tr){

    open OUT, ">$outdir/gene2tr.txt" || die $!;
    my $gene2tr = &gene2tr("$infile");
    foreach my $k (sort keys %{$gene2tr}){
        foreach (@{$gene2tr->{$k}}){
            print OUT "$k\t$_\n";
        }
    }
    close OUT;
}


sub gff2gtf{
    my $gff = shift @_;
    my (%gene, %gene_info, %mRNA, %exon, %cds);
    open GFF, "< $gff" || die $!;
    while(<GFF>){
        chomp;
        next if /^$|^#/;
        my @a = split /\t/, $_;

        if ($a[2] =~ /gene/i){
            my ($id) = $a[8] =~ /ID=(gene\d+)/;
            my ($gene) = $a[8] =~ /Dbxref=GeneID:(\d+)/;
            my ($gene_type) = $a[8] =~ /gene_biotype=([^;]+)/;
            $gene_info{$gene} = [$a[0], $a[3], $a[4], $a[6], $id, $gene_type];
            $gene{$gene} = $id;

        }elsif($a[2] =~ /mRNA|lnc_RNA|trascript/i){
            my ($rna_id) = $a[8] =~ /ID=(rna\d+)/;
            my ($gene) = $a[8] =~ /Dbxref=GeneID:(\d+)/;
            my ($transcript) = $a[8] =~ /transcript_id=([^;]+)/;
            $mRNA{$gene}{$rna_id}{mRNA} = [$a[0], $a[3], $a[4], $a[6], $transcript];

        }elsif($a[2] =~ /exon/){
            my ($exon_id) = $a[8] =~ /ID=(id\d+)/;
            my ($rna_id) = $a[8] =~ /Parent=([^;]+)/;
            $exon{$rna_id}{$exon_id} = [$a[0], $a[3], $a[4], $a[6]];

        }elsif($a[2] =~ /CDS/i){
            my ($cds_id) = $a[8] =~ /ID=(cds\d+)/;
            my ($rna_id) = $a[8] =~ /Parent=([^;]+)/;
            my ($protein)= $a[8] =~ /protein_id=([^;]+)/;
            $cds{$rna_id}{$cds_id}{pep} = $protein;
            push @{$cds{$rna_id}{$cds_id}{pos}},"$a[3]\t$a[4]";

        }elsif($a[2] =~ /rRNA/i){
            my ($rna_id) = $a[8] =~ /ID=(rna\d+)/;
            my ($gene) = $a[8] =~ /Dbxref=GeneID:(\d+)/;
            $mRNA{$gene}{$rna_id}{ncRNA} = [$a[0], $a[3], $a[4], $a[6], "rRNA"];
            
        }elsif($a[2] =~ /tRNA/i){
            my ($rna_id) = $a[8] =~ /ID=(rna\d+)/;
            my ($gene) = $a[8] =~ /Dbxref=GeneID:(\d+)/;
            $mRNA{$gene}{$rna_id}{ncRNA} = [$a[0], $a[3], $a[4], $a[6], "tRNA"];
            
        }else{
            next;
        }
    }
    return (\%gene, \%gene_info, \%mRNA, \%exon, \%cds);
}


sub gene2tr{
    my $gff = shift @_;
    my %out;
    open GFF, "<$gff" || die $!;
    while(<GFF>){
        chomp;
        next if (/^$|^#/);
        my @a = split /\t/, $_;
        if ($a[2] =~ /mRNA|transcript|lnc_RNA/i){
            my ($gene) = $a[8] =~ /Parent=(gene\d+)/;
            my ($transcript) = $a[8] =~ /transcript_id=(\w+_\d+\.\d+)/;
            push @{$out{$gene}}, $transcript;
        }else{
            next;
        }
    }
    close GFF;
    
    return \%out;
}

sub usage{
    print STDERR<<EOF;
===========================================================
Name:
    $0
Usage:
    perl $0 [options]
Options:
    -infile*    input file [gff/gff3/gtf]
    -outdir     outdir for ouputs files
    -fasta      fasta sequence [-ncr/-cds requared]
    -gene2tr    out gene2tr file
    -gff2gtf    convert gff/gff3 to gtf format
    -ncrna      fetch ncRNA tpye [defuall: all, rRNA/tRNA/ncRNA]
    -cds        fetch cds sequence 
    -help|?     print this help information
e.g:
    perl $0 -infile boluo.gff -outdir ./ -gene2tr -gff2gtf
===========================================================
EOF
    exit 1;
}

__END__
Author : liupeng@genebang.com
Date   : Fri Dec 15 09:55:07 CST 2017
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末荒椭,一起剝皮案震驚了整個濱河市出嘹,隨后出現(xiàn)的幾起案子唯咬,更是在濱河造成了極大的恐慌纱注,老刑警劉巖,帶你破解...
    沈念sama閱讀 217,185評論 6 503
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件胆胰,死亡現(xiàn)場離奇詭異狞贱,居然都是意外死亡,警方通過查閱死者的電腦和手機蜀涨,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 92,652評論 3 393
  • 文/潘曉璐 我一進(jìn)店門瞎嬉,熙熙樓的掌柜王于貴愁眉苦臉地迎上來,“玉大人厚柳,你說我怎么就攤上這事氧枣。” “怎么了别垮?”我有些...
    開封第一講書人閱讀 163,524評論 0 353
  • 文/不壞的土叔 我叫張陵便监,是天一觀的道長。 經(jīng)常有香客問我宰闰,道長茬贵,這世上最難降的妖魔是什么? 我笑而不...
    開封第一講書人閱讀 58,339評論 1 293
  • 正文 為了忘掉前任移袍,我火速辦了婚禮解藻,結(jié)果婚禮上,老公的妹妹穿的比我還像新娘葡盗。我一直安慰自己螟左,他們只是感情好,可當(dāng)我...
    茶點故事閱讀 67,387評論 6 391
  • 文/花漫 我一把揭開白布觅够。 她就那樣靜靜地躺著胶背,像睡著了一般。 火紅的嫁衣襯著肌膚如雪喘先。 梳的紋絲不亂的頭發(fā)上钳吟,一...
    開封第一講書人閱讀 51,287評論 1 301
  • 那天,我揣著相機與錄音窘拯,去河邊找鬼红且。 笑死,一個胖子當(dāng)著我的面吹牛涤姊,可吹牛的內(nèi)容都是我干的暇番。 我是一名探鬼主播,決...
    沈念sama閱讀 40,130評論 3 418
  • 文/蒼蘭香墨 我猛地睜開眼思喊,長吁一口氣:“原來是場噩夢啊……” “哼壁酬!你這毒婦竟也來了?” 一聲冷哼從身側(cè)響起,我...
    開封第一講書人閱讀 38,985評論 0 275
  • 序言:老撾萬榮一對情侶失蹤舆乔,失蹤者是張志新(化名)和其女友劉穎岳服,沒想到半個月后,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體蜕煌,經(jīng)...
    沈念sama閱讀 45,420評論 1 313
  • 正文 獨居荒郊野嶺守林人離奇死亡派阱,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 37,617評論 3 334
  • 正文 我和宋清朗相戀三年,在試婚紗的時候發(fā)現(xiàn)自己被綠了斜纪。 大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片贫母。...
    茶點故事閱讀 39,779評論 1 348
  • 序言:一個原本活蹦亂跳的男人離奇死亡,死狀恐怖盒刚,靈堂內(nèi)的尸體忽然破棺而出腺劣,到底是詐尸還是另有隱情,我是刑警寧澤因块,帶...
    沈念sama閱讀 35,477評論 5 345
  • 正文 年R本政府宣布橘原,位于F島的核電站,受9級特大地震影響涡上,放射性物質(zhì)發(fā)生泄漏趾断。R本人自食惡果不足惜,卻給世界環(huán)境...
    茶點故事閱讀 41,088評論 3 328
  • 文/蒙蒙 一吩愧、第九天 我趴在偏房一處隱蔽的房頂上張望芋酌。 院中可真熱鬧,春花似錦雁佳、人聲如沸脐帝。這莊子的主人今日做“春日...
    開封第一講書人閱讀 31,716評論 0 22
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽堵腹。三九已至,卻和暖如春星澳,著一層夾襖步出監(jiān)牢的瞬間疚顷,已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 32,857評論 1 269
  • 我被黑心中介騙來泰國打工禁偎, 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留腿堤,地道東北人。 一個月前我還...
    沈念sama閱讀 47,876評論 2 370
  • 正文 我出身青樓届垫,卻偏偏與公主長得像释液,于是被迫代替她去往敵國和親全释。 傳聞我的和親對象是個殘疾皇子装处,可洞房花燭夜當(dāng)晚...
    茶點故事閱讀 44,700評論 2 354

推薦閱讀更多精彩內(nèi)容