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 ************** |