G'day bisimen,
You wrote:
I'm stuck.
So, I'm trying to make a program which looks at fasta formated data, and counts the number of words with X length. Say, I'm looking for two words, it will count all two words it can find in each of the sequences. This counting part, kinda works so far, thanks to @toolic who answered my previous post, but it's not really that accurate since it's not actually looking for X length words, it's just splitting the sequence in X parts (ATGCAT -> AT GC AT) and counts all of the different parts it gets...
Anyway, here is an example data of what I want my program to work with:
>ENSCAPT00000000001.1 cdna (...) CACTCCGCCGGCAAGGCCACCAAGTATCTGAAGGATGTCACGTTAAAGAAGCAATGTGTG CCATTCCGGCGTTACAATGGTGGAGTTGGAAGGTGTGCCCAGGCCAAACAATGGGGCTGG ACCCAGGGTCGGTGGCCCAAAAAGAGTGCTGAATTTTTGCTGCAC >ENSCAPT00000000002.1 cdna (...) CAGCTGCTGAAGACGGAGCTGGGGTCCTTCTTCACCGAGTACCTGCAGAACCAGCTGCTG ACCAAAGGCATGGTGATCCTTCGGGACAAGATCCGCTTCTATGAGGGACAGAAACTGCTG GACTCGCTGGCAGAGACCTGGGACTTCTTCTTCAGCGACGTGCTGCCCACGCTGCAGGCC ATCTTCTACCC >ENSCAPT00000000004.1 cdna (...) ATGATGCTGGGAGAAGATGATGAGGAGTTCGTGGTGAAGGTGCGGGGTTTGCCTTGGTCC TGCTCGGCTGACGAGGTGCAGCGGTTCTTCTCCGACTGCAAAATTCAAAATGGTGCTCAA GGTATTTGTTTCATCTACACCAGAGAAGGCAGACCGAGTGGCGAGGCTTTTGTTGAACTT GAATCCGAAGATGAAGTCAAACTATCAAACAACGTTGAAATGGATTGGGTGTTGAAGCAT ACTGGTCCAAATAGTCCTGACACGGCCAATGATGGCTTTGTACGGCTTAGAGAACTCCCC TTTGGATGTAGCAAnd my counting-part of my code so far:
#Open data file; chomp($paths[1]); my ($sfile) = $paths[1].$ARGV[0]; open my $fain,'<', $sfile or die "Unable to open $sfile\nError: $!\n"; my $length = $paths[3]; my $line = 0; my %count; while (my $fasta = <$fain>){ if ($fasta !~ m/^>/){ while ($fasta =~ /(.{$length})/g) { $count{$1}++; } } if ($fasta =~ m/^>/){ #Print ID line print "$_ $count{$_}\n" for sort keys %count; %count = (); $line++; chomp($fasta); print "$line|$fasta\n"; } } print "$_ $count{$_}\n" for sort keys %count;It's weird, I know. But I'm trying... I guess I should count them in a different way, but the main thing that bothers me now, is that I want to get more data, from the data. See, I want to know how long (how many characters) each sequence is. And I also want to manipulate the hash data. I want to detect anomalies. If the data was random, you would expect to find an equal amount of all words. So in a 100 character long sequence, AA is expected to be seen 6.25 (I think) times. But if I count a 100 character sequence, and I find 25 AA's then I could do something like 1-(6.25/25) = 0.75 more AA than expected. Since it really doesn't make sense to compare amount of AA's found in a 50c.long. seq. with AA's found in a 500c.long.seq.
I feel my main obstacle here is that the sequence is on more than one line.
Thanks for any answers!
Then you deleted the entire body of your post! Don't do that; instead, do what it says in "How do I change/delete my post?".
"I feel my main obstacle here is that the sequence is on more than one line."
You can use this technique:
#!/usr/bin/env perl -l use strict; use warnings; { local $/ = "\n>"; while (<DATA>) { $_ = substr $_, 1 if $. == 1; chomp; my ($key, @seq_parts) = split /\n/; print "KEY=|$key|"; my $seq = join '', @seq_parts; print "SEQ=|$seq|"; } } __DATA__ >ENSCAPT00000000001.1 cdna (...) CACTCCGCCGGCAAGGCCACCAAGTATCTGAAGGATGTCACGTTAAAGAAGCAATGTGTG CCATTCCGGCGTTACAATGGTGGAGTTGGAAGGTGTGCCCAGGCCAAACAATGGGGCTGG ACCCAGGGTCGGTGGCCCAAAAAGAGTGCTGAATTTTTGCTGCAC >ENSCAPT00000000002.1 cdna (...) CAGCTGCTGAAGACGGAGCTGGGGTCCTTCTTCACCGAGTACCTGCAGAACCAGCTGCTG ACCAAAGGCATGGTGATCCTTCGGGACAAGATCCGCTTCTATGAGGGACAGAAACTGCTG GACTCGCTGGCAGAGACCTGGGACTTCTTCTTCAGCGACGTGCTGCCCACGCTGCAGGCC ATCTTCTACCC >ENSCAPT00000000004.1 cdna (...) ATGATGCTGGGAGAAGATGATGAGGAGTTCGTGGTGAAGGTGCGGGGTTTGCCTTGGTCC TGCTCGGCTGACGAGGTGCAGCGGTTCTTCTCCGACTGCAAAATTCAAAATGGTGCTCAA GGTATTTGTTTCATCTACACCAGAGAAGGCAGACCGAGTGGCGAGGCTTTTGTTGAACTT GAATCCGAAGATGAAGTCAAACTATCAAACAACGTTGAAATGGATTGGGTGTTGAAGCAT ACTGGTCCAAATAGTCCTGACACGGCCAATGATGGCTTTGTACGGCTTAGAGAACTCCCC TTTGGATGTAGCA
Output:
KEY=|ENSCAPT00000000001.1 cdna (...)| SEQ=|CACTCCGCCGGCAAGGCCACCAAGTATCTGAAGGATGTCACGTTAAAGAAGCAATGTGTGCCATT +CCGGCGTTACAATGGTGGAGTTGGAAGGTGTGCCCAGGCCAAACAATGGGGCTGGACCCAGGGTCGGTG +GCCCAAAAAGAGTGCTGAATTTTTGCTGCAC| KEY=|ENSCAPT00000000002.1 cdna (...)| SEQ=|CAGCTGCTGAAGACGGAGCTGGGGTCCTTCTTCACCGAGTACCTGCAGAACCAGCTGCTGACCAA +AGGCATGGTGATCCTTCGGGACAAGATCCGCTTCTATGAGGGACAGAAACTGCTGGACTCGCTGGCAGA +GACCTGGGACTTCTTCTTCAGCGACGTGCTGCCCACGCTGCAGGCCATCTTCTACCC| KEY=|ENSCAPT00000000004.1 cdna (...)| SEQ=|ATGATGCTGGGAGAAGATGATGAGGAGTTCGTGGTGAAGGTGCGGGGTTTGCCTTGGTCCTGCTC +GGCTGACGAGGTGCAGCGGTTCTTCTCCGACTGCAAAATTCAAAATGGTGCTCAAGGTATTTGTTTCAT +CTACACCAGAGAAGGCAGACCGAGTGGCGAGGCTTTTGTTGAACTTGAATCCGAAGATGAAGTCAAACT +ATCAAACAACGTTGAAATGGATTGGGTGTTGAAGCATACTGGTCCAAATAGTCCTGACACGGCCAATGA +TGGCTTTGTACGGCTTAGAGAACTCCCCTTTGGATGTAGCA|
With respect to the counting part, about 10 days ago I wrote "Re: Identifying Overlapping Matches in Nucleotide Sequence". That also talks about efficiency when using biological data; subsequent discussion compares benchmarks of this and other solutions.
By the way, in the code you originally posted,
%count = ();
removes all data from %count on every iteration of the while loop: not really what you want to do!
— Ken
In reply to Re: Counting & manipulating hashes
by kcott
in thread Counting & manipulating hashes
by bisimen
| For: | Use: | ||
| & | & | ||
| < | < | ||
| > | > | ||
| [ | [ | ||
| ] | ] |