scripts/utilities/mask_by_search.PLS

   1 #!perl -w
   2 # Author:  Jason Stajich <jason-at-bioperl-dot-org>
   3
   4
   5 =head1 NAME
   6
   7 mask_by_search - mask sequence(s) based on its alignment results
   8
   9 =head1 SYNOPSIS
  10
  11   mask_by_search.pl -f blast genomefile blastfile.bls > maskedgenome.fa
  12
  13 =head1 DESCRIPTION
  14
  15 Mask sequence based on significant alignments of another sequence.
  16 You need to provide the report file and the entire sequence data which
  17 you want to mask.  By default this will assume you have done a TBLASTN
  18 (or TFASTY) and try and mask the hit sequence assuming you've provided
  19 the sequence file for the hit database.  If you would like to do the
  20 reverse and mask the query sequence specify the -t/--type query flag.
  21
  22 This is going to read in the whole sequence file into memory so for
  23 large genomes this may fall over.  I'm using DB_File to prevent
  24 keeping everything in memory, one solution is to split the genome into
  25 pieces (BEFORE you run the DB search though, you want to use the exact
  26 file you BLASTed with as input to this program).
  27
  28 Below the double dash (--) options are of the form
  29 --format=fasta or --format fasta
  30 or you can just say
  31 -f fasta
  32
  33 By -f/--format I mean either are acceptable options.  The =s or =n
  34 or =c specify these arguments expect a 'string'
  35
  36 Options:
  37     -f/--format=s    Search report format (fasta,blast,axt,hmmer,etc)
  38     -sf/--sformat=s  Sequence format (fasta,genbank,embl,swissprot)
  39     --hardmask       (booelean) Hard mask the sequence
  40                      with the maskchar [default is lowercase mask]
  41     --maskchar=c     Character to mask with [default is N], change
  42                      to 'X' for protein sequences
  43     -e/--evalue=n    Evalue cutoff for HSPs and Hits, only
  44                      mask sequence if alignment has specified evalue
  45                      or better
  46     -o/--out/
  47     --outfile=file   Output file to save the masked sequence to.
  48     -t/--type=s      Alignment seq type you want to mask, the
  49                      'hit' or the 'query' sequence. [default is 'hit']
  50     --minlen=n       Minimum length of an HSP for it to be used
  51                      in masking [default 0]
  52     -h/--help        See this help information
  53
  54 =head1 AUTHOR - Jason Stajich
  55
  56 Jason Stajich, jason-at-bioperl-dot-org.
  57
  58 =cut
  59
  60
  61 use strict;
  62 use Bio::SeqIO;
  63 use Bio::SearchIO;
  64 use Getopt::Long;
  65 use Bio::Seq;
  66 use DB_File;
  67 # assuming tblastn or tfasty type alignment
  68
  69 my $format = 'blast';
  70 my $sformat= undef;
  71 my $evalue = undef;
  72 my $type   = 'hit';
  73 my $minlen = 50;
  74 my $hardmask = 0; # mask with $maskchar instead of lowercase
  75 my $maskchar = 'N'; # if we hard mask, mask with this cahr
  76 my $outfile;
  77 GetOptions(
  78            'f|format:s'  => \$format,
  79            'sf|sformat:s'=> \$sformat,
  80            'hardmask'    => \$hardmask,
  81            'maskchar:s'  => \$maskchar,
  82            'e|evalue:s'  => \$evalue,
  83            'o|out|outfile:s' => \$outfile,
  84            't|type:s'    => \$type,
  85            'minlen:s'    => \$minlen,
  86            'h|help'      => sub { system('perldoc', $0);
  87                                   exit; },
  88            );
  89 if( $type !~ /^(hit|query)/i ) {
  90     die("type must be query or hit[default] not $type") ;
  91 }
  92 $type = lc($type);
  93
  94 if(length($maskchar) > 1 ) {
  95     die("expected a mask character, not a string (you gave $maskchar)");
  96 }
  97 my $genomefile = shift || die('need a file containing the genome');
  98 my $reportfile = shift;
  99
 100 # this could be problem for large genomes, figure out a
 101 # better way to do this later on
 102 # or force people to split it up
 103 my $genomeparser = new Bio::SeqIO(-file  => $genomefile,
 104                                   -format=> $sformat);
 105 my %seqs;
 106 unlink('/tmp/genome.idx');
 107 tie(%seqs,'DB_File','/tmp/genome.idx');
 108 while( my $seq = $genomeparser->next_seq ) {
 109     # should we pre-force to upper case?
 110     $seqs{$seq->display_id} = $seq->seq();
 111 }
 112
 113 my $parser = new Bio::SearchIO(-file   => $reportfile,
 114                                -format => $format);
 115
 116 while( my $r = $parser->next_result ) {
 117     while( my $h = $r->next_hit ) {
 118         last if( defined $evalue && $h->significance > $evalue );
 119         my $hname = $h->name;
 120         if( ! $seqs{$hname} ) {
 121             die("Cannot find sequence $hname in genome seq");
 122         }
 123         while( my $hsp = $h->next_hsp ) {
 124             last if( defined $evalue && $hsp->evalue > $evalue );
 125             next if( $hsp->length('total') < $minlen);
 126             my ($s,$len) = ( $hsp->$type()->start,
 127                              $hsp->$type()->length);
 128
 129             if( $hardmask ) {
 130                 substr($seqs{$hname}, $s,$len, $maskchar x $len);
 131             } else {
 132                 substr($seqs{$hname}, $s,$len,
 133                        lc(substr($seqs{$hname}, $s,$len)));
 134             }
 135         }
 136     }
 137 }
 138
 139 my $out;
 140 if( $outfile ) {
 141     $out = new Bio::SeqIO(-file   => ">$outfile",
 142                           -format => $sformat);
 143 } else {
 144     $out = new Bio::SeqIO(-format => $sformat);
 145 }
 146
 147 while( my ($seqname,$seq) = each %seqs ) {
 148     $out->write_seq(Bio::Seq->new(-seq        => $seq,
 149                                   -display_id => $seqname,
 150                                   -description=> 'MASKED'));
 151 }
 152 END {
 153     unlink('/tmp/genome.idx');
 154 }