bp_bin/assemble_contigs

   1 #!/usr/bin/env perl
   2
   3 # Copyright (C) 2007-2010 Martin A. Hansen.
   4
   5 # This program is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU General Public License
   7 # as published by the Free Software Foundation; either version 2
   8 # of the License, or (at your option) any later version.
   9
  10 # This program is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 # GNU General Public License for more details.
  14
  15 # You should have received a copy of the GNU General Public License
  16 # along with this program; if not, write to the Free Software
  17 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18
  19 # http://www.gnu.org/copyleft/gpl.html
  20
  21
  22 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  23
  24 # Calculate coverage of reads mapped to a backbone sequence.
  25
  26 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  27
  28
  29 use warnings;
  30 use strict;
  31
  32 use Data::Dumper;
  33 use Maasha::Biopieces;
  34 use Maasha::Filesys;
  35 use Maasha::UCSC::Wiggle;
  36
  37
  38 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  39
  40
  41 my ( $options, $in, $out, $record, $contig_hash, $i, $s_id, $strand, $new_record, $contig, $beg, $end, @array, $min, $max, $mean );
  42
  43 $options = Maasha::Biopieces::parse_options(
  44     [
  45         { long => 'strand', short => 's', type => 'flag', mandatory => 'no', default => undef, allowed => undef, disallowed => undef },
  46         { long => 'bridge', short => 'b', type => 'flag', mandatory => 'no', default => undef, allowed => undef, disallowed => undef },
  47     ]
  48 );
  49
  50 $in  = Maasha::Biopieces::read_stream( $options->{ "stream_in" } );
  51 $out = Maasha::Biopieces::write_stream( $options->{ "stream_out" } );
  52
  53 while ( $record = Maasha::Biopieces::get_record( $in ) )
  54 {
  55     if ( $record->{ 'S_ID' } and $record->{ 'STRAND' } )
  56     {
  57         if ( $options->{ 'strand' } ) {
  58             next if not $record->{ 'STRAND' };
  59         } else {
  60             $record->{ 'STRAND' } = 'X';
  61         }
  62
  63         if ( $record->{ 'Q_ID' } and $record->{ 'S_BEG' } and $record->{ 'S_END' } )
  64         {
  65             for ( $i = $record->{ 'S_BEG' }; $i <= $record->{ 'S_END' }; $i++ ) {
  66                 $contig_hash->{ $record->{ 'S_ID' } }->{ $record->{ 'STRAND' } }->[ $i ]++;
  67             }
  68         }
  69
  70         if ( $record->{ 'Q_ID1' } and $record->{ 'S_BEG1' } and $record->{ 'S_END1' } and
  71              $record->{ 'Q_ID2' } and $record->{ 'S_BEG2' } and $record->{ 'S_END2' }
  72         )
  73         {
  74             if ( $options->{ 'bridge' } )
  75             {
  76                 ( $record->{ 'S_BEG1' }, $record->{ 'S_END2' } ) = ( $record->{ 'S_BEG2' }, $record->{ 'S_END1' } ) if ( $record->{ 'S_BEG1' } > $record->{ 'S_END2' } );
  77
  78                 for ( $i = $record->{ 'S_BEG1' }; $i <= $record->{ 'S_END2' }; $i++ ) {
  79                     $contig_hash->{ $record->{ 'S_ID' } }->{ $record->{ 'STRAND' } }->[ $i ]++;
  80                 }
  81             }
  82             else
  83             {
  84                 for ( $i = $record->{ 'S_BEG1' }; $i <= $record->{ 'S_END1' }; $i++ ) {
  85                     $contig_hash->{ $record->{ 'S_ID' } }->{ $record->{ 'STRAND' } }->[ $i ]++;
  86                 }
  87
  88                 for ( $i = $record->{ 'S_BEG2' }; $i <= $record->{ 'S_END2' }; $i++ ) {
  89                     $contig_hash->{ $record->{ 'S_ID' } }->{ $record->{ 'STRAND' } }->[ $i ]++;
  90                 }
  91             }
  92         }
  93     }
  94
  95     # Maasha::Biopieces::put_record( $record, $out );
  96 }
  97
  98 $contig = 0;
  99
 100 foreach $s_id ( keys %{ $contig_hash } )
 101 {
 102     foreach $strand ( keys %{ $contig_hash->{ $s_id } } )
 103     {
 104         $beg = 0;
 105         $end = 0;
 106
 107         while ( ( $beg, $end ) = Maasha::UCSC::Wiggle::fixedstep_scan( $contig_hash->{ $s_id }->{ $strand }, $beg ) )
 108         {
 109             @array = @{ $contig_hash->{ $s_id }->{ $strand } }[ $beg .. $end - 1 ];
 110
 111             ( $min, $max ) = Maasha::Calc::minmax( \@array );
 112             $mean          = Maasha::Calc::mean( \@array );
 113
 114             $new_record->{ 'S_ID' }        = $s_id;
 115             $new_record->{ 'STRAND' }      = $strand if $options->{ 'strand' };
 116             $new_record->{ 'S_BEG' }       = $beg;
 117             $new_record->{ 'S_END' }       = $end;
 118             $new_record->{ 'CONTIG_ID' }   = $contig;
 119             $new_record->{ 'CONTIG_LEN' }  = $end - $beg + 1;
 120             $new_record->{ 'CONTIG_MIN' }  = $min;
 121             $new_record->{ 'CONTIG_MAX' }  = $max;
 122             $new_record->{ 'CONTIG_MEAN' } = sprintf "%.2f", $mean;
 123             $new_record->{ 'CONTIG' }      = join ( ";", @array );
 124
 125             $beg = $end + 1;
 126             $contig++;
 127
 128             Maasha::Biopieces::put_record( $new_record, $out );
 129         }
 130     }
 131 }
 132
 133 Maasha::Biopieces::close_stream( $in );
 134 Maasha::Biopieces::close_stream( $out );
 135
 136
 137 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 138
 139
 140 BEGIN
 141 {
 142     Maasha::Biopieces::status_set();
 143 }
 144
 145
 146 END
 147 {
 148     Maasha::Biopieces::status_log();
 149 }
 150
 151
 152 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 153
 154
 155 __END__