Mercurial > repos > chawhwa > neuma

comparison NEUMA-1.2.1/filter.best.from.bowtieout.2.pl @ 0:c44c43d185ef draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
NEUMA-1.2.1 Uploaded
author chawhwa
date Thu, 08 Aug 2013 00:46:13 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:c44c43d185ef
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 }