Beefy Boxes and Bandwidth Generously Provided by pair Networks
Problems? Is your data what you think it is?
 
PerlMonks  

Re^2: Mustering Reads

by neversaint (Deacon)
on Nov 02, 2010 at 07:30 UTC ( [id://868944]=note: print w/replies, xml ) Need Help??


in reply to Re: Mustering Reads
in thread Mustering Reads

Dear BrowserUK,
I was trying the following input to your code above:
my $k = 15; my %repo = ( "readA" => "AAACATGAAAAGAAATGATGAAACAGA", "readB" => "AAATCATGAAACAGAGCCTCATCTCCC ", ); my @order = ( "readB_1","readA_2");
Why it gives this output
AAACATGAAAAGGAGCCTCATCTCCC GCCTCATCTCCC
Instead of:
AAACATGAAAAGAAAT(G,C)ATGAAACAGAGCCTCATCTCCC
Following this alignment:
readA AAACATGAAAAGAAATGATGAAACAGA readB AAATCATGAAACAGAGCCTCATCTCCC
How can I fix it?

---
neversaint and everlastingly indebted.......

Replies are listed 'Best First'.
Re^3: Mustering Reads
by BrowserUk (Patriarch) on Nov 02, 2010 at 08:17 UTC
    How can I fix it?

    By giving better examples :) I had assumed on the basis of the sample supplied, that the overlapping sections were identical as in both cases supplied. I think javafans make that same assumption.

    You'd need to accumulate the overlaps in a third array in the while, and then add another for loop to process them is the same way as @heads & @tails. (There obviously scope there for making the for loop a subroutine.

    I'm done for today, but I'll be back tonight.

Re^3: Mustering Reads
by BrowserUk (Patriarch) on Nov 02, 2010 at 14:34 UTC

    As I suggested:

    #! perl -slw use strict; use Data::Dump qw[ pp ]; my $k = 15; my %repo = ( "readA" => "AAACATGAAAAGAAATGATGAAACAGA", "readB" => "AAATCATGAAACAGAGCCTCATCTCCC ", ); my @order = ( "readB_1","readA_2"); #my $k = 25; #my %repo = ( # "readA" => "GCTGAGGCAGGAGAATTGCTTGAACCTGGGAGGCA", # "readB" => "TACTCAGGAGGCTGAGGCAGGAGAATTGCTTGAAC", # "readC" => "GCTGAGGCAGGAGAATTGCTTGAACTTAGGGGATG", # "readD" => "TACTCGGGAGGCTGAGGCAGGAGAATTGCTTGAAC", #); #my @order = ( "readA_1", "readB_2", "readC_1", "readD_2"); sub muster { local $^W; my( $bits ) = @_; my $muster = ''; for my $p ( 0 .. length( $bits->[0] )-1 ) { my %uniq; ++$uniq{ substr $bits->[ $_ ], $p, 1 } for 0 .. $#$bits; if( keys %uniq > 1 ) { $muster .= '(' . join( ',', keys %uniq ) . ')'; } else { $muster .= each %uniq; } } return $muster; } my( @heads, @bodys, @tails ); while( @order ) { my( $s1, $p1, $s2, $p2 ) = map split( '_', shift @order ), 1 .. 2; ( $s1, $s2 ) = ( $s2, $s1 ) if $p1 > $p2; push @heads, substr $repo{ $s2 }, 0, length( $repo{ $s2 } ) - $k; push @bodys, substr $repo{ $s2 }, -$k; push @tails, substr $repo{ $s1 }, $k; push @bodys, substr $repo{ $s1 }, 0, $k; } my $head = muster( \@heads ); my $body = muster( \@bodys ); my $tail = muster( \@tails ); print $head, $body, $tail; __END__ TACTC(A,G)GGAGGAGAATTGCTTGAACCTGGGAGGCA(T,C)T(A,G)GG(A,G)G(A,G)(T,C)(A +,G) TACTC(A,G)GGAGGCTGAGGCAGGAGAATTGCTTGAAC(T,C)T(A,G)GG(A,G)G(A,G)(T,C)(A +,G) AAACATGAAAAGAAAT(C,G)ATGAAACAGAGCCTCATCTCCC AAACATGAAAAGAAAT(G,C)ATGAAACAGAGCCTCATCTCCC

    Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
    "Science is about questioning the status quo. Questioning authority".
    In the absence of evidence, opinion is indistinguishable from prejudice.

Log In?
Username:
Password:

What's my password?
Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: note [id://868944]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this?Last hourOther CB clients
Other Users?
Others exploiting the Monastery: (4)
As of 2024-03-29 10:28 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    No recent polls found