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 ( $record->{ 'S_BEG' } and $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 ( $record->{ 'S_BEG1' } and $record->{ 'S_END1' } and $record->{ 'S_BEG2' } and $record->{ 'S_END2' } )
  78         {
  79             if ( $options->{ 'bridge' } )
  80             {
  81                 ( $record->{ 'S_BEG1' }, $record->{ 'S_END2' } ) = ( $record->{ 'S_BEG2' }, $record->{ 'S_END1' } ) if ( $record->{ 'S_BEG1' } > $record->{ 'S_END2' } );
  82
  83                 for ( $i = $record->{ 'S_BEG1' }; $i <= $record->{ 'S_END2' }; $i++ ) {
  84                     $contig_hash->{ $record->{ 'S_ID' } }->{ $record->{ 'Q_ID' } }->{ $record->{ 'STRAND' } }->[ $i ]++;
  85                 }
  86             }
  87             else
  88             {
  89                 for ( $i = $record->{ 'S_BEG1' }; $i <= $record->{ 'S_END1' }; $i++ ) {
  90                     $contig_hash->{ $record->{ 'S_ID' } }->{ $record->{ 'Q_ID' } }->{ $record->{ 'STRAND' } }->[ $i ]++;
  91                 }
  92
  93                 for ( $i = $record->{ 'S_BEG2' }; $i <= $record->{ 'S_END2' }; $i++ ) {
  94                     $contig_hash->{ $record->{ 'S_ID' } }->{ $record->{ 'Q_ID' } }->{ $record->{ 'STRAND' } }->[ $i ]++;
  95                 }
  96             }
  97         }
  98     }
  99
 100     # Maasha::Biopieces::put_record( $record, $out );
 101 }
 102
 103 $contig = 0;
 104
 105 foreach $s_id ( keys %{ $contig_hash } )
 106 {
 107     foreach $q_id ( keys %{ $contig_hash->{ $s_id } } )
 108     {
 109         foreach $strand ( keys %{ $contig_hash->{ $s_id }->{ $q_id } } )
 110         {
 111             $beg = 0;
 112             $end = 0;
 113
 114             while ( ( $beg, $end ) = Maasha::UCSC::Wiggle::fixedstep_scan( $contig_hash->{ $s_id }->{ $q_id }->{ $strand }, $beg ) )
 115             {
 116                 @array = @{ $contig_hash->{ $s_id }->{ $q_id }->{ $strand } }[ $beg .. $end - 1 ];
 117
 118                 ( $min, $max ) = Maasha::Calc::minmax( \@array );
 119                 $mean          = Maasha::Calc::mean( \@array );
 120
 121                 $new_record->{ 'S_ID' }        = $s_id;
 122                 $new_record->{ 'Q_ID' }        = $q_id   if $options->{ 'q_id' };
 123                 $new_record->{ 'STRAND' }      = $strand if $options->{ 'strand' };
 124                 $new_record->{ 'S_BEG' }       = $beg;
 125                 $new_record->{ 'S_END' }       = $end;
 126                 $new_record->{ 'CONTIG_ID' }   = $contig;
 127                 $new_record->{ 'CONTIG_LEN' }  = $end - $beg + 1;
 128                 $new_record->{ 'CONTIG_MIN' }  = $min;
 129                 $new_record->{ 'CONTIG_MAX' }  = $max;
 130                 $new_record->{ 'CONTIG_MEAN' } = sprintf "%.2f", $mean;
 131                 $new_record->{ 'CONTIG' }      = join ( ";", @array );
 132
 133                 $beg = $end + 1;
 134                 $contig++;
 135
 136                 Maasha::Biopieces::put_record( $new_record, $out );
 137             }
 138         }
 139     }
 140 }
 141
 142 Maasha::Biopieces::close_stream( $in );
 143 Maasha::Biopieces::close_stream( $out );
 144
 145
 146 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 147
 148
 149 BEGIN
 150 {
 151     Maasha::Biopieces::status_set();
 152 }
 153
 154
 155 END
 156 {
 157     Maasha::Biopieces::status_log();
 158 }
 159
 160
 161 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 162
 163
 164 __END__