From a3a66d2f110793bff4a8feb7053c1b6b5c524c8b Mon Sep 17 00:00:00 2001 From: martinahansen Date: Mon, 21 Sep 2009 15:32:40 +0000 Subject: [PATCH] added assemble_contigs git-svn-id: http://biopieces.googlecode.com/svn/trunk@682 74ccb610-7750-0410-82ae-013aeee3265d --- bp_bin/assemble_contigs | 148 +++++++++++++++++++++++++++++++++++++ code_perl/Maasha/Matrix.pm | 25 +++++-- 2 files changed, 165 insertions(+), 8 deletions(-) create mode 100755 bp_bin/assemble_contigs diff --git a/bp_bin/assemble_contigs b/bp_bin/assemble_contigs new file mode 100755 index 0000000..3e73a76 --- /dev/null +++ b/bp_bin/assemble_contigs @@ -0,0 +1,148 @@ +#!/usr/bin/env perl + +# Copyright (C) 2007-2009 Martin A. Hansen. + +# This program is free software; you can redistribute it and/or +# modify it under the terms of the GNU General Public License +# as published by the Free Software Foundation; either version 2 +# of the License, or (at your option) any later version. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. + +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + +# http://www.gnu.org/copyleft/gpl.html + + +# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + +# Calculate coverage of reads mapped to a backbone sequence. + +# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +use warnings; +use strict; + +use Data::Dumper; +use Maasha::Biopieces; +use Maasha::Filesys; +use Maasha::UCSC::Wiggle; + + +# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +my ( $options, $in, $out, $record, $contig_hash, $i, $s_id, $strand, $new_record, $contig, $beg, $end, @array, $min, $max, $mean ); + +$options = Maasha::Biopieces::parse_options( + [ + { long => 'bridge', short => 'b', type => 'flag', mandatory => 'no', default => undef, allowed => undef, disallowed => undef }, + ] +); + +$in = Maasha::Biopieces::read_stream( $options->{ "stream_in" } ); +$out = Maasha::Biopieces::write_stream( $options->{ "stream_out" } ); + +while ( $record = Maasha::Biopieces::get_record( $in ) ) +{ + if ( $record->{ 'S_ID' } and $record->{ 'STRAND' } ) + { + if ( $record->{ 'Q_ID' } and $record->{ 'S_BEG' } and $record->{ 'S_END' } ) + { + for ( $i = $record->{ 'S_BEG' }; $i <= $record->{ 'S_END' }; $i++ ) { + $contig_hash->{ $record->{ 'S_ID' } }->{ $record->{ 'STRAND' } }->[ $i ]++; + } + } + + if ( $record->{ 'Q_ID1' } and $record->{ 'S_BEG1' } and $record->{ 'S_END1' } and + $record->{ 'Q_ID2' } and $record->{ 'S_BEG2' } and $record->{ 'S_END2' } + ) + { + if ( $options->{ 'bridge' } ) + { + ( $record->{ 'S_BEG1' }, $record->{ 'S_END2' } ) = ( $record->{ 'S_BEG2' }, $record->{ 'S_END1' } ) if ( $record->{ 'S_BEG1' } > $record->{ 'S_END2' } ); + + for ( $i = $record->{ 'S_BEG1' }; $i <= $record->{ 'S_END2' }; $i++ ) { + $contig_hash->{ $record->{ 'S_ID' } }->{ $record->{ 'STRAND' } }->[ $i ]++; + } + } + else + { + for ( $i = $record->{ 'S_BEG1' }; $i <= $record->{ 'S_END1' }; $i++ ) { + $contig_hash->{ $record->{ 'S_ID' } }->{ $record->{ 'STRAND' } }->[ $i ]++; + } + + for ( $i = $record->{ 'S_BEG2' }; $i <= $record->{ 'S_END2' }; $i++ ) { + $contig_hash->{ $record->{ 'S_ID' } }->{ $record->{ 'STRAND' } }->[ $i ]++; + } + } + } + } + + # Maasha::Biopieces::put_record( $record, $out ); +} + +$contig = 0; + +foreach $s_id ( keys %{ $contig_hash } ) +{ + foreach $strand ( keys %{ $contig_hash->{ $s_id } } ) + { + $beg = 0; + $end = 0; + + while ( ( $beg, $end ) = Maasha::UCSC::Wiggle::fixedstep_scan( $contig_hash->{ $s_id }->{ $strand }, $beg ) ) + { + @array = @{ $contig_hash->{ $s_id }->{ $strand } }[ $beg .. $end - 1 ]; + + ( $min, $max ) = Maasha::Calc::minmax( \@array ); + $mean = Maasha::Calc::mean( \@array ); + + $new_record->{ 'S_ID' } = $s_id; + $new_record->{ 'STRAND' } = $strand; + $new_record->{ 'CONTIG_ID' } = $contig; + $new_record->{ 'CONTIG_BEG' } = $beg; + $new_record->{ 'CONTIG_END' } = $end; + $new_record->{ 'CONTIG_LEN' } = $end - $beg + 1; + $new_record->{ 'CONTIG_MIN' } = $min; + $new_record->{ 'CONTIG_MAX' } = $max; + $new_record->{ 'CONTIG_MEAN' } = sprintf "%.2f", $mean; + $new_record->{ 'CONTIG' } = join ( ";", @array ); + + $beg = $end + 1; + $contig++; + + Maasha::Biopieces::put_record( $new_record, $out ); + } + } +} + +Maasha::Biopieces::close_stream( $in ); +Maasha::Biopieces::close_stream( $out ); + + +# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +BEGIN +{ + Maasha::Biopieces::status_set(); +} + + +END +{ + Maasha::Biopieces::status_log(); +} + + +# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +__END__ diff --git a/code_perl/Maasha/Matrix.pm b/code_perl/Maasha/Matrix.pm index 4a29c4d..ed90a25 100644 --- a/code_perl/Maasha/Matrix.pm +++ b/code_perl/Maasha/Matrix.pm @@ -906,11 +906,11 @@ sub list_check_sort } -sub list_shrink +sub list_deflate { # Martin A. Hansen, September 2009. - # Shrinks a list of values to a specified size + # Defaltes a list of values to a specified size # and at the same time average the values. my ( $list, @@ -919,19 +919,28 @@ sub list_shrink # Returns nothing. - my ( $old_size, $bucket_size, $i, @new_list ); + my ( $old_size, $bucket_size, $bucket_rest, $i, @new_list ); $old_size = scalar @{ $list }; - Maasha::Common::error( qq(can't shrink to a bigger list: $old_size < $new_size ) ) if $old_size < $new_size; + Maasha::Common::error( qq(Can't shrink to a bigger list: $old_size < $new_size ) ) if $old_size < $new_size; - $bucket_size = int( $old_size / $new_size ); + $bucket_size = int( $old_size / $new_size ); + $bucket_rest = $old_size - ( $new_size * $bucket_size ); - for ( $i = 0; $i < $old_size - $bucket_size + 1; $i += $bucket_size ) { - push @new_list, Maasha::Calc::mean( [ splice @{ $list }, $i, $bucket_size ] ); + print STDERR "old_size: $old_size new_size: $new_size bucket_size: $bucket_size bucket_rest: $bucket_rest\n"; + + $i = 0; + + while ( $i < $new_size ) + { + # push @new_list, [ @{ $list }[ $i * $bucket_size .. $i * $bucket_size + $bucket_size - 1 ] ]; + push @new_list, Maasha::Calc::mean( [ @{ $list }[ $i * $bucket_size .. $i * $bucket_size + $bucket_size - 1 ] ] ); + + $i ++; } - $list = \@new_list; + @{ $list } = @new_list; } -- 2.39.5