bp_bin/assemble_contigs

   1 #!/usr/bin/env perl
   2
   3 # Copyright (C) 2007-2009 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 => 'bridge', short => 'b', type => 'flag', mandatory => 'no', default => undef, allowed => undef, disallowed => undef },
  46     ]
  47 );
  48
  49 $in  = Maasha::Biopieces::read_stream( $options->{ "stream_in" } );
  50 $out = Maasha::Biopieces::write_stream( $options->{ "stream_out" } );
  51
  52 while ( $record = Maasha::Biopieces::get_record( $in ) )
  53 {
  54     if ( $record->{ 'S_ID' } and $record->{ 'STRAND' } )
  55     {
  56         if ( $record->{ 'Q_ID' } and $record->{ 'S_BEG' } and $record->{ 'S_END' } )
  57         {
  58             for ( $i = $record->{ 'S_BEG' }; $i <= $record->{ 'S_END' }; $i++ ) {
  59                 $contig_hash->{ $record->{ 'S_ID' } }->{ $record->{ 'STRAND' } }->[ $i ]++;
  60             }
  61         }
  62
  63         if ( $record->{ 'Q_ID1' } and $record->{ 'S_BEG1' } and $record->{ 'S_END1' } and
  64              $record->{ 'Q_ID2' } and $record->{ 'S_BEG2' } and $record->{ 'S_END2' }
  65         )
  66         {
  67             if ( $options->{ 'bridge' } )
  68             {
  69                 ( $record->{ 'S_BEG1' }, $record->{ 'S_END2' } ) = ( $record->{ 'S_BEG2' }, $record->{ 'S_END1' } ) if ( $record->{ 'S_BEG1' } > $record->{ 'S_END2' } );
  70
  71                 for ( $i = $record->{ 'S_BEG1' }; $i <= $record->{ 'S_END2' }; $i++ ) {
  72                     $contig_hash->{ $record->{ 'S_ID' } }->{ $record->{ 'STRAND' } }->[ $i ]++;
  73                 }
  74             }
  75             else
  76             {
  77                 for ( $i = $record->{ 'S_BEG1' }; $i <= $record->{ 'S_END1' }; $i++ ) {
  78                     $contig_hash->{ $record->{ 'S_ID' } }->{ $record->{ 'STRAND' } }->[ $i ]++;
  79                 }
  80
  81                 for ( $i = $record->{ 'S_BEG2' }; $i <= $record->{ 'S_END2' }; $i++ ) {
  82                     $contig_hash->{ $record->{ 'S_ID' } }->{ $record->{ 'STRAND' } }->[ $i ]++;
  83                 }
  84             }
  85         }
  86     }
  87
  88     # Maasha::Biopieces::put_record( $record, $out );
  89 }
  90
  91 $contig = 0;
  92
  93 foreach $s_id ( keys %{ $contig_hash } )
  94 {
  95     foreach $strand ( keys %{ $contig_hash->{ $s_id } } )
  96     {
  97         $beg = 0;
  98         $end = 0;
  99
 100         while ( ( $beg, $end ) = Maasha::UCSC::Wiggle::fixedstep_scan( $contig_hash->{ $s_id }->{ $strand }, $beg ) )
 101         {
 102             @array = @{ $contig_hash->{ $s_id }->{ $strand } }[ $beg .. $end - 1 ];
 103
 104             ( $min, $max ) = Maasha::Calc::minmax( \@array );
 105             $mean          = Maasha::Calc::mean( \@array );
 106
 107             $new_record->{ 'S_ID' }        = $s_id;
 108             $new_record->{ 'STRAND' }      = $strand;
 109             $new_record->{ 'CONTIG_ID' }   = $contig;
 110             $new_record->{ 'CONTIG_BEG' }  = $beg;
 111             $new_record->{ 'CONTIG_END' }  = $end;
 112             $new_record->{ 'CONTIG_LEN' }  = $end - $beg + 1;
 113             $new_record->{ 'CONTIG_MIN' }  = $min;
 114             $new_record->{ 'CONTIG_MAX' }  = $max;
 115             $new_record->{ 'CONTIG_MEAN' } = sprintf "%.2f", $mean;
 116             $new_record->{ 'CONTIG' }      = join ( ";", @array );
 117
 118             $beg = $end + 1;
 119             $contig++;
 120
 121             Maasha::Biopieces::put_record( $new_record, $out );
 122         }
 123     }
 124 }
 125
 126 Maasha::Biopieces::close_stream( $in );
 127 Maasha::Biopieces::close_stream( $out );
 128
 129
 130 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 131
 132
 133 BEGIN
 134 {
 135     Maasha::Biopieces::status_set();
 136 }
 137
 138
 139 END
 140 {
 141     Maasha::Biopieces::status_log();
 142 }
 143
 144
 145 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 146
 147
 148 __END__