NCBI gff/gff3 转换为 gtf 格式

适用于 NCBI 数据库中 gff 格式文件转换为 gtf 格式,使用 perl 代码实现:

#!/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
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,686评论 6 492
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,668评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 158,160评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,736评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,847评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,043评论 1 291
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,129评论 3 410
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,872评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,318评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,645评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,777评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,470评论 4 333
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,126评论 3 317
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,861评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,095评论 1 267
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,589评论 2 362
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,687评论 2 351

推荐阅读更多精彩内容