从gtf到gff

从Ensemble数据库中下载到的基因组坐标文件通常是gtf,在某些处理中有些不便,这里笔者提供一种把gtf转化成gff格式的方法,仅供参考。针对具体的文件,可以对脚本的细节做相关调整。


#!/usr/bin/perl
use strict;
use warnings;
die "use gtf\n" 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[3],$cut[4];
		$info{$id}=[$cut[0],$cut[1],$cut[6],$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;


  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值