Beefy Boxes and Bandwidth Generously Provided by pair Networks
"be consistent"

Re: unique sequences

by Cristoforo (Curate)
on Dec 12, 2017 at 21:35 UTC ( #1205371=note: print w/replies, xml ) Need Help??

in reply to unique sequences


I used BioPerl, Bio::SeqIO to load the fasta sequences all into a long string. Then I tested for kmers that only appeared once and printed them out.

The results agree with your desired output.

#!/usr/bin/perl use strict; use warnings; use Bio::SeqIO; my $in = Bio::SeqIO->new( -file => "dmel-all-chromosome-r6.18.fasta +" , -format => 'fasta'); my $fasta; while ( my $seq = $in->next_seq() ) { $fasta .= $seq->seq; } my %seen; for my $i (0 .. length($fasta) - 21) { my $kmer = substr $fasta, $i, 21; next unless substr($kmer, -2) eq 'GG'; my $match = substr($kmer, -12); $seen{$match}{count}++; $seen{$match}{kmer} = $kmer; } my $crispr; for my $key (keys %seen) { next unless $seen{$key}{count} == 1; print "crispr_", ++$crispr, "\n"; print $seen{$key}{kmer}, "\n"; } __END__ *** ouput crispr_1 TTTAGACTCCCCTTGTACAGG crispr_2 TCTTCAGTCTCCAGTCTCCGG crispr_3 TTGCGTTGCGGAGCATACTGG crispr_4 TGCCACCAGTGGTTCCAAGGG crispr_5 TTATGTTTGTACGAGGGGGGG crispr_6 TCTCTTTGGTTTACGGATGGG crispr_7 TTGGCAAGGAGACGGTCCTGG crispr_8 TGAATTAAAGCTTGCGCGAGG crispr_9 GGAAGAGGCATCAACGAGGGG crispr_10 TGCAGCGGCCTAACAAGGCGG crispr_11 CTGCCCGATCCTAACTCCAGG crispr_12 ATATATGTTTGACCGTCGGGG ... crispr_126892 TTGCTTGGCACTAAGGCGGGG crispr_126893 CACCAAAAAGGACTTGCGTGG crispr_126894 GTGCCCCTCACTCATGCGGGG crispr_126895 TAAAAAGCGACGCAGTATTGG crispr_126896 CCGGATTTCTTCGTACAGGGG crispr_126897 GGTGGCTATGCTATGGTACGG crispr_126898 CTGCGTTGATGTTAGGTAGGG crispr_126899 GCTGGGACCCGAATACGTAGG

This program runs in under 2 minutes with the string of fasta characters about 145M. From your specs it looks like you want to combine the sequences into one string to test for kmers that only occur 1 time.

Update: Just realized how odd it was for my results to agree with the ones you posted. The order of my results were in the (un)ordered keys from the hash. It is curious why the order I got agreed with the order in your post!

Also, I think the reason you were getting 1 million results rather than 250,000 is you are testing for uniqueness on the whole 21 char windowsize instead of testing just the last 12 chars (including the 'GG' ending). There would be more unique 21 char kmers than unique 12 char kmers.

Log In?

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

How do I use this? | Other CB clients
Other Users?
Others avoiding work at the Monastery: (6)
As of 2021-04-15 15:27 GMT
Find Nodes?
    Voting Booth?

    No recent polls found