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