|
1
|
1 use strict;
|
|
|
2
|
|
|
3 my $bowtiefilename=$ARGV[0];
|
|
|
4 my $readsfilename=$ARGV[1];
|
|
|
5 my $elandfilename=$ARGV[2];
|
|
|
6
|
|
|
7 open(readsfile, $readsfilename);
|
|
|
8
|
|
|
9 my $detectformat=`head -c 1 $readsfilename`;
|
|
|
10
|
|
|
11 #my $firstletter=$detectformat;
|
|
|
12 #my $looplinenumbers=4;
|
|
|
13
|
|
|
14 #$looplinenumbers=2 if ($detectformat eq ">");
|
|
|
15 open(bowtiefile, $bowtiefilename);
|
|
|
16 open(elandfile, ">".$elandfilename);
|
|
|
17 my $readfilelinenum=0;
|
|
|
18 # hash the positions of the alignments for each read id
|
|
|
19 my %readposhash;
|
|
|
20 my $bowtiepos = tell (bowtiefile);
|
|
|
21 while (my $bowtieline=<bowtiefile>)
|
|
|
22 {
|
|
|
23 my ($bowtiereadname, $strand, $mapped_id, $pos, $seq, $qt,$num, $mapinfo)=split("\t",$bowtieline);
|
|
|
24 if (not exists $readposhash{$bowtiereadname} )
|
|
|
25 {
|
|
|
26 $readposhash{$bowtiereadname} = $bowtiepos;
|
|
|
27 }
|
|
|
28 $bowtiepos = tell (bowtiefile);
|
|
|
29 }
|
|
|
30
|
|
|
31 while(my $readline=<readsfile>)
|
|
|
32 {
|
|
|
33 $readfilelinenum++;
|
|
|
34 if(($readline=~/^$detectformat/) && ($readfilelinenum%2 ==1))
|
|
|
35 {
|
|
|
36 chomp($readline);
|
|
|
37 my $readname=substr($readline, 1, length($readline)-1);
|
|
|
38 if( not exists $readposhash{$readname} )
|
|
|
39 {
|
|
|
40 print elandfile $readname,"\tNA\tNM\n";
|
|
|
41 next;
|
|
|
42 }
|
|
|
43 else
|
|
|
44 {
|
|
|
45 my @mapped_ids=();
|
|
|
46 my @mapped_pos=();
|
|
|
47 my @mapped_strand=();
|
|
|
48 seek(bowtiefile, $readposhash{$readname}, 0);
|
|
|
49 while (my $bowtieline=<bowtiefile>)
|
|
|
50 {
|
|
|
51 my ($bowtiereadname, $strand, $mapped_id, $pos, $seq, $qt,$num, $mapinfo)=split("\t",$bowtieline);
|
|
|
52 if($readname eq $bowtiereadname)
|
|
|
53 {
|
|
|
54 push(@mapped_ids, $mapped_id);
|
|
|
55 push(@mapped_pos,$pos);
|
|
|
56 push(@mapped_strand,$strand);
|
|
|
57 }
|
|
|
58 else
|
|
|
59 {
|
|
|
60 last;
|
|
|
61 }
|
|
|
62
|
|
|
63 }
|
|
|
64 print elandfile $readname,"\t";
|
|
|
65 print elandfile "NA\t";
|
|
|
66 print elandfile scalar(@mapped_ids),":0:0\t";
|
|
|
67 for(my $i=0;$i<@mapped_ids;$i++)
|
|
|
68 {
|
|
|
69 print elandfile "/",$mapped_ids[$i];
|
|
|
70 print elandfile ":",$mapped_pos[$i]+1;
|
|
|
71 if($mapped_strand[$i] eq "+")
|
|
|
72 {
|
|
|
73 print elandfile "F0,";
|
|
|
74 }
|
|
|
75 else
|
|
|
76 {
|
|
|
77 print elandfile "R0,";
|
|
|
78 }
|
|
|
79
|
|
|
80 }
|
|
|
81 print elandfile "\n";
|
|
|
82
|
|
|
83 }
|
|
|
84 }
|
|
|
85 }
|
|
|
86
|
|
|
87 close(elandfile);
|
|
|
88 close(bowtiefile);
|
|
|
89 close(readsfile);
|
|
|
90
|
|
|
91 exit;
|
|
|
92 while(my $bowtieline=<bowtiefile>)
|
|
|
93 {
|
|
|
94 my ($bowtiereadname, $strand, $mapped_id, $pos, $seq, $qt,$num, $mapinfo)=split("\t",$bowtieline);
|
|
|
95 while(my $readline=<readsfile>)
|
|
|
96 {
|
|
|
97 $readfilelinenum++;
|
|
|
98 if(($readline=~/^$detectformat/) && ($readfilelinenum%2 ==1))
|
|
|
99 #if($readline=~/^$detectformat/)
|
|
|
100 {
|
|
|
101 chomp($readline);
|
|
|
102 my $readname=substr($readline, 1, length($readline)-1);
|
|
|
103
|
|
|
104
|
|
|
105 if($readname ne $bowtiereadname)
|
|
|
106 {
|
|
|
107 print elandfile $readname,"\tNA\tNM\n";
|
|
|
108 next;
|
|
|
109 }
|
|
|
110 else
|
|
|
111 {
|
|
|
112 my @mapped_ids=();
|
|
|
113 my @mapped_pos=();
|
|
|
114 my @mapped_strand=();
|
|
|
115 push(@mapped_ids, $mapped_id);
|
|
|
116 push(@mapped_pos,$pos);
|
|
|
117 push(@mapped_strand,$strand);
|
|
|
118 while(1)
|
|
|
119 {
|
|
|
120 $bowtieline=<bowtiefile>;
|
|
|
121 my ($bowtiereadname, $strand, $mapped_id, $pos, $seq, $qt,$num, $mapinfo)=split("\t",$bowtieline);
|
|
|
122 if( $bowtiereadname eq $readname )
|
|
|
123 {
|
|
|
124 push(@mapped_ids, $mapped_id);
|
|
|
125 push(@mapped_pos,$pos);
|
|
|
126 push(@mapped_strand,$strand);
|
|
|
127 }
|
|
|
128 else
|
|
|
129 {
|
|
|
130 seek(bowtiefile, -1*length($bowtieline),1);
|
|
|
131 print elandfile $readname,"\t";
|
|
|
132 print elandfile "NA\t";
|
|
|
133 print elandfile scalar(@mapped_ids),":0:0\t";
|
|
|
134 for(my $i=0;$i<@mapped_ids;$i++)
|
|
|
135 {
|
|
|
136 print elandfile "/",$mapped_ids[$i];
|
|
|
137 print elandfile ":",$mapped_pos[$i]+1;
|
|
|
138 if($mapped_strand[$i] eq "+")
|
|
|
139 {
|
|
|
140 print elandfile "F0,";
|
|
|
141 }
|
|
|
142 else
|
|
|
143 {
|
|
|
144 print elandfile "R0,";
|
|
|
145 }
|
|
|
146
|
|
|
147 }
|
|
|
148 print elandfile "\n";
|
|
|
149 last;
|
|
|
150 }
|
|
|
151 }
|
|
|
152 last;
|
|
|
153
|
|
|
154 }
|
|
|
155 }
|
|
|
156 }
|
|
|
157 }
|
|
|
158
|
|
|
159 while(my $readline=<readsfile>)
|
|
|
160 {
|
|
|
161 $readfilelinenum++;
|
|
|
162 if(($readline=~/^$detectformat/) && ($readfilelinenum%2 ==1))
|
|
|
163 #if($readline=~/^$detectformat/)
|
|
|
164 {
|
|
|
165 chomp($readline);
|
|
|
166 my $readname=substr($readline, 1, length($readline)-1);
|
|
|
167 print elandfile $readname,"\tNA\tNM\n";
|
|
|
168 }
|
|
|
169 }
|
|
|
170
|
|
|
171 close(elandfile);
|
|
|
172 close(bowtiefile);
|
|
|
173
|
|
|
174
|
|
|
175 close(readsfile);
|
