0
|
1 #!/usr/bin/perl
|
|
2 use Getopt::Long;
|
|
3
|
|
4 # This version fixed the problem of not being able to merge when the order of the two mates were flipped. (version 2)
|
|
5
|
|
6
|
|
7 my $delete_original;
|
|
8 my $newbowtiefile; ## output bowtiefile
|
|
9 my $mmcol=4;
|
|
10 &GetOptions ( 'd|delete_original' => \$delete_original,
|
|
11 'o=s' => \$newbowtiefile,
|
|
12 'mmcol=i' => \$mmcol ## column for mismatch info (0-based)
|
|
13 );
|
|
14
|
|
15 if(@ARGV<2) { &help; exit; }
|
|
16
|
|
17 my $bowtiefile = shift @ARGV;
|
|
18 my $PE = shift @ARGV; # 1 for paired-end, 0 for single-end
|
|
19
|
|
20 if(!defined($newbowtiefile)){
|
|
21 ($newbowtiefile) = $bowtiefile=~/(.+)\.bowtieout/;
|
|
22 $newbowtiefile .= ".best.bowtieout";
|
|
23 }
|
|
24
|
|
25
|
|
26
|
|
27
|
|
28 my $prev_head1=''; my $prev_head2='';
|
|
29 my $totmm=100000; # just some huge number
|
|
30 my @mapping=();
|
|
31 open $IN, $bowtiefile or die "Can't open bowtiefile $bowtiefile\n";
|
|
32 open $OUT, ">$newbowtiefile" or die "Can't write to output bowtiefile $newbowtiefile\n";
|
|
33 while(<$IN>){
|
|
34 chomp;
|
|
35 my ($head1,$strand1,$tr,$pos1,$mmstr1) = (split/\t/)[0,1,2,3,$mmcol];
|
|
36 my $mm1 = &parse_mmstr($mmstr1);
|
|
37 if($PE==1) {
|
|
38 chomp(my $secondline = <$IN>);
|
|
39 ($head2,$strand2,$pos2,$mmstr2) = (split/\t/,$secondline)[0,1,3,$mmcol];
|
|
40 my $mm2 = &parse_mmstr($mmstr2);
|
|
41 }
|
|
42
|
|
43 if($PE==1){
|
|
44 if($prev_head1 eq $head1 || $prev_head2 eq $head1) {
|
|
45 if($totmm>$mm1+$mm2) { @mapping=(); push @mapping,[$tr,$strand1,$pos1,$mmstr1,$strand2,$pos2,$mmstr2]; $totmm=$mm1+$mm2; }
|
|
46 elsif($totma==$mm1+$mm2) { push @mapping,[$tr,$strand1,$pos1,$mmstr1,$strand2,$pos2,$mmstr2]; }
|
|
47 }
|
|
48 else {
|
|
49 &print_mapping_PE(\@mapping,$OUT, $prev_head1,$prev_head2);
|
|
50 @mapping=();
|
|
51 push @mapping,[$tr,$strand1,$pos1,$mmstr1,$strand2,$pos2,$mmstr2];
|
|
52 $totmm=$mm1+$mm2;
|
|
53 }
|
|
54 }
|
|
55 else {
|
|
56 if($prev_head1 eq $head1) {
|
|
57 if($totmm>$mm1) { @mapping=(); push @mapping,[$tr,$strand1,$pos1,$mmstr1]; $totmm=$mm1; }
|
|
58 elsif($totmm==$mm1) { push @mapping,[$tr,$strand1,$pos1,$mmstr1]; }
|
|
59 }
|
|
60 else {
|
|
61 &print_mapping_SE(\@mapping,$OUT, $prev_head1);
|
|
62 @mapping=();
|
|
63 push @mapping,[$tr,$strand1,$pos1,$mmstr1];
|
|
64 $totmm=$mm1;
|
|
65 }
|
|
66 }
|
|
67
|
|
68 $prev_head1= $head1;
|
|
69 if($PE==1) { $prev_head2 = $head2; }
|
|
70 }
|
|
71
|
|
72
|
|
73 ## last line
|
|
74 if($PE==1){
|
|
75 &print_mapping_PE(\@mapping,$OUT, $prev_head1, $prev_head2);
|
|
76 }
|
|
77 else {
|
|
78 &print_mapping_SE(\@mapping,$OUT, $prev_head1);
|
|
79 }
|
|
80
|
|
81 close $OUT;
|
|
82 close $IN;
|
|
83
|
|
84
|
|
85
|
|
86 ## delete original bowtieout file
|
|
87 if(defined $delete_original) { `rm -f $bowtiefile`; }
|
|
88
|
|
89
|
|
90
|
|
91
|
|
92
|
|
93 ## subroutines @@
|
|
94
|
|
95
|
|
96 sub parse_mmstr {
|
|
97 if(!defined $_[0]){ return 0; }
|
|
98 else {
|
|
99 return(scalar(split/,/,$_[0]));
|
|
100 }
|
|
101 }
|
|
102
|
|
103
|
|
104 sub print_mapping_PE {
|
|
105 my $pMapping = shift @_;
|
|
106 my $OUT = shift @_;
|
|
107 my $prev_head1 = shift @_;
|
|
108 my $prev_head2 = shift @_;
|
|
109
|
|
110 if($prev_head1 ne ''){
|
|
111 for my $i (0..$#$pMapping){
|
|
112 my ($tr,$strand1,$pos1,$mm1,$strand2,$pos2,$mm2) = @{$pMapping->[$i]};
|
|
113 print $OUT "$prev_head1\t$strand1\t$tr\t$pos1\t$mm1\n";
|
|
114 print $OUT "$prev_head2\t$strand2\t$tr\t$pos2\t$mm2\n";
|
|
115 }
|
|
116 }
|
|
117 }
|
|
118
|
|
119
|
|
120
|
|
121 sub print_mapping_SE {
|
|
122 my $pMapping = shift @_;
|
|
123 my $OUT = shift @_;
|
|
124 my $prev_head1 = shift @_;
|
|
125
|
|
126 if($prev_head1 ne ''){
|
|
127 for my $i (0..$#$pMapping){
|
|
128 my ($tr,$strand1,$pos1,$mm1) = @{$pMapping->[$i]};
|
|
129 print $OUT "$prev_head1\t$strand1\t$tr\t$pos1\t$mm1\n";
|
|
130 }
|
|
131 }
|
|
132 }
|
|
133
|
|
134
|
|
135
|
|
136 sub help {
|
|
137 print STDERR<<EOF
|
|
138
|
|
139 usage: $0 <options> bowtieout_file PE(1/0)
|
|
140
|
|
141 Options
|
|
142 -d|--delete_original : delete the original bowtieout file
|
|
143 -o <outfilename> : output bowtieout file name (default : originalfiletitle.best.bowtieout )
|
|
144
|
|
145 EOF
|
|
146 }
|