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
|
#!/usr/bin/env perl
use strict;
use warnings;
use lib ("/usr/lib/trinityrnaseq/PerlLib");
use Gene_obj;
use Fasta_reader;
use GFF3_utils;
use Carp;
my $usage = "usage: $0 genes.gff3 genome.fasta\n\n";
my $genes_gff3 = $ARGV[0] or die $usage;
my $genome_fasta_file = $ARGV[1] or die $usage;
my $fasta_reader = new Fasta_reader($genome_fasta_file);
my %genome = $fasta_reader->retrieve_all_seqs_hash();
my $gene_obj_indexer_href = {};
## associate gene identifiers with contig id's.
my $contig_to_gene_list_href = &GFF3_utils::index_GFF3_gene_objs($genes_gff3, $gene_obj_indexer_href);
foreach my $asmbl_id (sort keys %$contig_to_gene_list_href) {
my $genome_seq = $genome{$asmbl_id} or die "Error, cannot find sequence for $asmbl_id"; #cdbyank_linear($asmbl_id, $fasta_db);
my @gene_ids = @{$contig_to_gene_list_href->{$asmbl_id}};
foreach my $gene_id (@gene_ids) {
my $gene_obj_ref = $gene_obj_indexer_href->{$gene_id};
my $orientation = $gene_obj_ref->get_orientation();
my $intron_text = "";
foreach my $isoform ($gene_obj_ref, $gene_obj_ref->get_additional_isoforms()) {
my @intron_coords = $isoform->get_intron_coordinates();
my $trans_id = $isoform->{Model_feat_name};
foreach my $intron (@intron_coords) {
my ($end5, $end3) = @$intron;
my ($lend, $rend) = sort {$a<=>$b} ($end5, $end3);
my $left_splice_dinuc = substr($genome_seq, $lend-1, 2);
my $right_splice_dinuc = substr($genome_seq, $rend-1-1, 2);
$intron_text .= join("\t", $gene_id, $trans_id, $asmbl_id, "$lend-$rend", $orientation, "$left_splice_dinuc..$right_splice_dinuc") . "\n";
}
}
if ($intron_text) {
print "$intron_text\n";
}
}
}
exit(0);
|
