从gtf到gff
发布时间:2020-12-15 23:54:48 所属栏目:大数据 来源:网络整理
导读:从Ensemble数据库中下载到的基因组坐标文件通常是gtf,在某些处理中有些不便,这里笔者提供一种把gtf转化成gff格式的方法,仅供参考。针对具体的文件,可以对脚本的细节做相关调整。 #!/usr/bin/perluse strict;use warnings;die "use gtfn" unless @ARGV==
从Ensemble数据库中下载到的基因组坐标文件通常是gtf,在某些处理中有些不便,这里笔者提供一种把gtf转化成gff格式的方法,仅供参考。针对具体的文件,可以对脚本的细节做相关调整。 #!/usr/bin/perl use strict; use warnings; die "use gtfn" unless @ARGV==1; my $gtf=shift; my %gene; my %mRNA; my %info; $gtf=~/.gz$/ ? (open IN,"gzip -cd $gtf|"||die) : (open IN,$gtf||die); while (<IN>){ chomp; next if /^#/; my @cut=split/t/,$_; my $id=(split/s+/,$cut[8])[4]; $id=~s/"//g; $id=~s/;//; my $gene=(split/s+/,$cut[8])[2]; $gene=~s/"//g; $gene=~s/;//; if($cut[2] eq "CDS" || $cut[2] eq "exon"){ my $join=join("t",($cut[0],$cut[1],"CDS",$cut[3],$cut[4],$cut[5],$cut[6],$cut[7],"Parent=$id;")) if $cut[2] eq "CDS"; push @{$gene{$id}},$join if $join; push @{$mRNA{$id}},$cut[4]; $info{$id}=[$cut[0],$gene]; } } close IN; my @output; foreach my $key (keys %gene){ @{$mRNA{$key}}=sort{$a<=>$b} @{$mRNA{$key}}; my $out=join "t",$info{$key}[0],$info{$key}[1],"mRNA",$mRNA{$key}[0],$mRNA{$key}[-1],".",$info{$key}[2],"ID=$key;Gene=$info{$key}[3]"; push @output,$out; for(my $i=0;$i<@{$gene{$key}};$i++){ push @output,$gene{$key}[$i]; } } my @output_sort=map{$_->[3]} sort{$a->[0] cmp $b->[0] || $a->[1] <=> $b->[1]||$a->[2] cmp $b->[2]} map{[(split/t/,$_)[0,3,2],$_]}@output; print "$_n" for@output_sort; (编辑:李大同) 【声明】本站内容均来自网络,其相关言论仅代表作者个人观点,不代表本站立场。若无意侵犯到您的权利,请及时与联系站长删除相关内容! |