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
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
|
#!/usr/bin/perl
#
# parse the output of sim4, fiter the alignments and
# output them to stdout
#
# Mario Stanke, August 17th, 2005
use SplicedAlignment;
use strict;
my $usage = "$0 -- filter the output of sim4\n";
$usage .= "Usage: $0 sim4.out seqfile.fa\n";
$usage .= "output to stdout. first messages. then complete genes. then genes complete only at the 5' end\n";
if ($#ARGV != 1) {
die "Unknown option\n\n$usage";
}
my $simfilename = $ARGV[0];
my $seqfilename = $ARGV[1];
my @alignmentList=();
###########################################################################################
#
# read in the sim4 alignments
#
###########################################################################################
open(SIM, "<$simfilename") or die ("Could not open $simfilename");
my ($sa, $estName, $estLength, $genomicName, $genomicLength, $estExonBegin, $estExonEnd, $exonBegin, $exonEnd, $intronFlag, $quality, $complement);
my @genomicMatches;
$/="--------\n";
while (<SIM>) {
next unless /seq1 =/;
/>(.*)\n/;
$estName = $1;
/seq1 = .*, (\d+) bp/;
$estLength = $1;
@genomicMatches = split(/seq2 = /, $_);
shift @genomicMatches;
foreach my $match (@genomicMatches) {
$match =~ /\((.*)\), (\d+) bp/;
$genomicName = $1;
$genomicLength = $2;
$complement = ($match =~ /complement/);
$sa = new SplicedAlignment($estName, $estLength, $genomicName, $genomicLength, $complement);
foreach (split (/\n/, $match)){
if (/(\d+)-(\d+) +\((\d+)-(\d+)\) +(\d+)%(.*)/) {
$estExonBegin = $1;
$estExonEnd = $2;
$exonBegin = $3;
$exonEnd = $4;
$quality = $5 / 100;
$intronFlag = $6;
$intronFlag =~ s/ //;
$sa->addExon($estExonBegin-1, $estExonEnd-1, $exonBegin-1, $exonEnd-1 , $quality, $intronFlag);
}
}
push @alignmentList, $sa;
}
}
my @geneList=();
###########################################################################################
#
# filter and enrich the data
#
###########################################################################################
foreach $sa (@alignmentList){
$sa->checkAlignment();
print STDERR $sa->get_status(), "\n";
if ($sa->get_status() eq "OK"){
my $seq = getSequence($sa->get_contigname);
$sa->makeGene($seq);
$sa->findSplicedUTR();
print STDERR "\t\t", $sa->get_status(), "\n";
if ($sa->get_status() eq "OK"){
push @geneList, $sa;
}
}
}
###########################################################################################
#
# output the data
#
###########################################################################################
print "### complete genes\n";
foreach $sa (@geneList){
if ($sa->get_complete5prime() && $sa->get_complete3prime()){
print $sa->output;
}
}
print "### genes incomplete at the 3' end and complete at the 5' end\n";
foreach $sa (@geneList){
if ($sa->get_complete5prime() && !$sa->get_complete3prime()){
print $sa->output;
}
}
###########################################################################################
#
# subroutines
#
###########################################################################################
sub getSequence {
my $seqname = shift;
my $seq;
# let cdbyank use the index to search the sequence
system("echo '$seqname' | cdbyank ${seqfilename}.cidx -d $seqfilename > genomic.fa");
system ("perl -i.orig -p -e 's/^Incorrectly.*fixed.\n//' genomic.fa");
open (SEQ, "<genomic.fa") or die ("Could not open genomic.fa");
$/="\n";
my @lines=<SEQ>;
shift @lines;
$seq = join ("", @lines);
$seq =~ s/\n//g;
return $seq;
}
|
