zzh
/
bwa
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

useless script

This commit is contained in:
Heng Li 2011-11-19 21:48:54 -05:00
parent 182cb2e89c
commit 84aa3fa696
1 changed files with 0 additions and 69 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

69
pairsam.pl
View File

@ -1,69 +0,0 @@
#!/usr/bin/perl -w
use strict;
use warnings;
use Getopt::Std;
my %opts = (a=>700, b=>100);
getopts('a:b:', \%opts);
die("Usage: pairsam.pl [-a $opts{a}] <read1.sam> <read2.sam>\n") if (@ARGV < 2);
my ($fh0, $fh1, $l0, $l1);
open($fh0, $ARGV[0] =~ /\.gz$/? "gzip -dc $ARGV[0] |" : $ARGV[0]) || die;
open($fh1, $ARGV[1] =~ /\.gz$/? "gzip -dc $ARGV[1] |" : $ARGV[1]) || die;
while ($l0 = <$fh0>) { last if $l0 !~ /^@/; print $l0 }
while ($l1 = <$fh1>) { last if $l1 !~ /^@/ }
while (defined($l0) && defined($l1)) {
my ($r0, $r1) = &pair_line(\$l0, \$l1, $opts{a}, $opts{b});
while ($l0 = <$fh0>) { last if $l0 !~ /^$r0/; print $l0; }
while ($l1 = <$fh1>) { last if $l1 !~ /^$r1/; print $l1; }
}
close($fh0); close($fh1);
sub pair_line {
my ($l0, $l1, $max_ins, $min_ins) = @_;
my @t0 = split("\t", $$l0);
my @t1 = split("\t", $$l1);
my ($n0, $n1) = ($t0[0], $t1[0]);
my ($cigar, $a0, $a1, $p0, $p1);
# length in alignment
$cigar = $t0[5]; $a0 = 0; $cigar =~ s/(\d+)[MI]/$a0 += $1/eg;
$cigar = $t1[5]; $a1 = 0; $cigar =~ s/(\d+)[MI]/$a1 += $1/eg;
# 5'-end alignment position on the read
$p0 = $t0[1] == 16? $t0[3] + $a0 : $t0[3];
$p1 = $t1[1] == 16? $t1[3] + $a1 : $t1[3];
# adjust mapping quality
if ($t0[2] eq $t1[2] && $t0[1]+$t1[1] == 16) { # on the same chr and forward-reverse
if (abs($p0 - $p1) <= $max_ins && abs($p0 - $p1) >= $min_ins) { # within the right insert size distribution
$t0[1] |= 2; $t1[1] |= 2; # flag as paired
if ($t0[4] < $t1[4]) { # increase mapQ
$t0[4] = $t0[4] + 10 < $t1[4]? $t0[4] + 10 : $t1[4];
} else {
$t1[4] = $t1[4] + 10 < $t0[4]? $t1[4] + 10 : $t0[4];
}
}
}
unless ($t0[1]&2) { # decrease mapQ if unpaired
$t0[4] = $t0[4] > 10? $t0[4] - 10 : 0;
$t1[4] = $t1[4] > 10? $t1[4] - 10 : 0;
}
# strip off /[12]
$t0[0] =~ s/\/[12]$//; $t1[0] =~ s/\/[12]$//;
# update FLAG
$t0[1] |= 0x41 | (($t1[1]&16)? 0x20 : 0) | (($t1[1]&4)? 0x8 : 0);
$t1[1] |= 0x81 | (($t0[1]&16)? 0x20 : 0) | (($t0[1]&4)? 0x8 : 0);
# update mate positions
if ($t0[2] eq $t1[2]) {
$t0[6] = $t1[6] = '=';
$t0[8] = $p1 - $p0; $t1[8] = $p0 - $p1;
} else { $t0[6] = $t1[2]; $t1[6] = $t0[2]; }
$t0[7] = $t1[3]; $t1[7] = $t0[3];
# print out
print join("\t", @t0);
print join("\t", @t1);
return ($n0, $n1);
}