forked from vikas0633/perl
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathGeneID_to_gff3.pl
87 lines (57 loc) · 1.98 KB
/
GeneID_to_gff3.pl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
#!/usr/bin/env perl
use strict;
use warnings;
use FindBin;
use lib ("$FindBin::Bin/../../PerlLib");
use Gene_obj;
my $model_type = "GeneID";
my $usage = "usage: $0 geneID.gff\n\n";
my $input_file = $ARGV[0] or die $usage;
main: {
my %data;
## parse input file
open (my $fh, $input_file) or die "Error, cannot open file $input_file";
while (<$fh>) {
chomp;
my @x = split(/\t/);
my $scaff = $x[0];
my $feat_type = $x[2];
my $lend = $x[3];
my $rend = $x[4];
my $orient = $x[6];
my $model_id = $x[8];
if ($feat_type eq "exon") {
## add feature
my $model_token = join("-", $scaff, $model_id);
my $exon_feature = { scaff => $scaff,
lend => $lend,
rend => $rend,
orient => $orient,
model_id => $model_id,
};
push (@{$data{$model_token}}, $exon_feature);
}
}
close $fh;
## write gff3 output
foreach my $model_id (sort keys %data) {
my @exons = @{$data{$model_id}};
my $scaff = $exons[0]->{scaff};
my $orient = $exons[0]->{orient};
my %coords;
foreach my $exon (@exons) {
my $lend = $exon->{lend};
my $rend = $exon->{rend};
my ($end5, $end3) = ($orient eq '+') ? ($lend, $rend) : ($rend, $lend);
$coords{$end5} = $end3;
}
my $gene_obj = new Gene_obj();
$gene_obj->populate_gene_object(\%coords, \%coords);
$gene_obj->{asmbl_id} = $scaff;
$gene_obj->{TU_feat_name} = "t.$model_id";
$gene_obj->{Model_feat_name} = $model_id;
$gene_obj->{com_name} = "GeneID prediction $model_id";
print $gene_obj->to_GFF3_format(source => "GeneID") . "\n";
}
exit(0);
}