in reply to Re^9: mismatching characters in dna sequence
in thread mismatching characters in dna sequence
#!/opt/local/bin/perl use Inline C => Config => BUILD_NOISY => 1; use Inline C => <<'END_C', NAME => 'dnacmp', CLEAN_AFTER_BUILD => 0; void dnacmp( char *a, char *b ) { int n = 0; Inline_Stack_Vars; Inline_Stack_Reset; while( *a && *b ) { if( *a != *b ) { Inline_Stack_Push( sv_2mortal( newSVpvf( "%d:%c:%c", n, *a +, *b ) ) ); } a++; b++; n++; } Inline_Stack_Done; return; } SV *dnacmp2( char *a, char *b ) { SV *res = newSVpv( "", 0 ); int n = 0; while( *a && *b ) { if( *a != *b ) sv_catpvf( res, "%d>%c:%c ", n, *a, *b ); ++a, ++b, ++n; } return res; } END_C use Bio::Seq; use Bio::DB::Sam; use Data::Dump qw[ pp ]; my $bamfile = Bio::DB::Sam->new( -bam => "eg1.bam" ); my @allReads = $bamfile->features( -type => 'match' ); my %stats; for my $read (@allReads) { my ( $ref, $matches, $query ) = $read->padded_alignment; m[\d+:(.:.)] and ++$stats{ $1 } for dnacmp2( $ref, $query ); print $read->qname, "\n", dnacmp2( $ref, $query ), "\n"; } pp \%stats;
|
|---|
| Replies are listed 'Best First'. | |
|---|---|
|
Re^11: mismatching characters in dna sequence
by BrowserUk (Patriarch) on Dec 30, 2011 at 22:00 UTC | |
by prbndr (Acolyte) on Dec 30, 2011 at 22:24 UTC | |
by BrowserUk (Patriarch) on Dec 30, 2011 at 22:36 UTC | |
by prbndr (Acolyte) on Dec 30, 2011 at 22:42 UTC |