| 1 |
#!/usr/bin/perl |
| 2 |
use strict; |
| 3 |
use Getopt::Std; |
| 4 |
use FindBin;use lib $FindBin::Bin; |
| 5 |
|
| 6 |
my $usage = q/Usage: |
| 7 |
rnaseq_mapflt.pl [-x <chromosome>] [-l <seedlen>] [-n <max_seed_mismatches>] \ |
| 8 |
[-t <max_terminal_mismatches>] <bowtie_map_data> |
| 9 |
|
| 10 |
Filters bowtie map output (sorted with best hits first) such that only |
| 11 |
mappings with <max_mismatches> at 3' end (after <seedlen>) are |
| 12 |
allowed. |
| 13 |
|
| 14 |
Use option -x to exclude mappings from a specific chromosome or contig |
| 15 |
|
| 16 |
Defaults are : -l16 -n1 -t3 -xchrM |
| 17 |
(exclude mitochondrial chromosome) |
| 18 |
|
| 19 |
For bowtie output created with -n and --best --strata it will also make sure |
| 20 |
that only the "best strata" is kept based on the total number of mismatches. |
| 21 |
|
| 22 |
/; |
| 23 |
umask 0002; |
| 24 |
getopts('x:t:l:n:o:') || die($usage."\n"); |
| 25 |
my $seedlen=$Getopt::Std::opt_l||16; |
| 26 |
my $maxmmseed=$Getopt::Std::opt_n || 1; |
| 27 |
my $maxmmend=$Getopt::Std::opt_t || 3; |
| 28 |
my $nochr=$Getopt::Std::opt_x || 'chrM'; |
| 29 |
my $outfile=$Getopt::Std::opt_o; |
| 30 |
if ($outfile) { |
| 31 |
open(OUTF, '>'.$outfile) || die("Error creating output file $outfile\n"); |
| 32 |
select(OUTF); |
| 33 |
} |
| 34 |
# -- |
| 35 |
my ($prevq, $prevmmcount, $prevmmseedcount); |
| 36 |
while (<>) { |
| 37 |
my $line=$_; |
| 38 |
chomp; |
| 39 |
next unless $_; |
| 40 |
my ($qry, $strand, $chr, $cstart, $seq, $qv, $d, $mism)=split(/\t/); |
| 41 |
#next if $chr eq $nochr; |
| 42 |
my $len=length($seq); |
| 43 |
next unless $len>8; |
| 44 |
my @mm=map { [split(/\:/)] } (split(/\,/,$mism)) ; |
| 45 |
#mismatches as a list of [mmcoord, nt_change] |
| 46 |
if ($strand eq '-') { #for reverse complement, reverse and adjust |
| 47 |
@mm=map { [$len-$$_[0]-1, $$_[1]] } reverse @mm; |
| 48 |
} |
| 49 |
my @mmseed=grep { $$_[0]<$seedlen } @mm; |
| 50 |
next if @mmseed>$maxmmseed; |
| 51 |
#my @mmbefore_end=grep { $$_[0]<$len-$maxmmend-1 } @mm; |
| 52 |
#next if @mmbefore_end>0; |
| 53 |
next if scalar(@mm)-scalar(@mmseed)>$maxmmend; |
| 54 |
next if ($prevq eq $qry && @mm>$prevmmcount); |
| 55 |
$prevq=$qry; |
| 56 |
$prevmmcount=scalar(@mm); |
| 57 |
$prevmmseedcount=scalar(@mmseed); |
| 58 |
#TODO: trim 5' end if any mismatches are grouped there close to that end |
| 59 |
my $tlen=$len; |
| 60 |
#print STDERR ">> testing $line"; |
| 61 |
foreach my $md (reverse @mm) { |
| 62 |
#print STDERR " .. testing $$md[0] vs $tlen-2\n"; |
| 63 |
if ($$md[0]<$tlen-2) { last; } |
| 64 |
else { $tlen=$$md[0]; } #truncate at mismatch |
| 65 |
} |
| 66 |
if ($tlen<$len) { |
| 67 |
#adjust mismatches and read/qual strings |
| 68 |
my $dlen=$len-$tlen; |
| 69 |
@mm=grep { $$_[0]<$tlen } @mm; |
| 70 |
if ($strand eq '-') { |
| 71 |
@mm = map { [$len-$$_[0]-1, $$_[1]] } reverse @mm; |
| 72 |
$seq=substr($seq,$dlen); |
| 73 |
$qv=substr($qv, $dlen); |
| 74 |
$cstart+=$dlen; |
| 75 |
#print STDERR "$qry\n" if $tlen>=$seedlen; |
| 76 |
} |
| 77 |
else { #truncate sequence |
| 78 |
$seq=substr($seq,0,$tlen); |
| 79 |
$qv=substr($qv,0,$tlen); |
| 80 |
#print STDERR "$qry (+)\n" if @mm>0 && $tlen>=$seedlen; |
| 81 |
} |
| 82 |
next if $tlen<$seedlen; |
| 83 |
my @mmstr = map { $$_[0].':'.$$_[1] } @mm; |
| 84 |
$line=join("\t", $qry, $strand, $chr, $cstart, $seq, $qv, $d, join(',',@mmstr))."\n"; |
| 85 |
} |
| 86 |
print $line unless $chr eq $nochr; |
| 87 |
} |
| 88 |
|
| 89 |
# -- |
| 90 |
if ($outfile) { |
| 91 |
select(STDOUT); |
| 92 |
close(OUTF); |
| 93 |
} |
| 94 |
|
| 95 |
#************ Subroutines ************** |
