How do I calculate tetranucleotide frequencies?

[AndrewRGross]

I'd like to calculate the tetranucleotide frequencies of my contigs and the bin them based on shared tetranucleotide frequencies. Can anyone recommend a program or script to calculate this, and possibly tools for binning?

My data is a low-diversity environmental metagenome. It consists of 150 bp paired end Illumina library generated on a HiSeq 2500, aligned in MIRA.

What I understand so far:
The occurrence of certain four-base sequences (ATAT, GGCA, TAAC) is not random within an organism, but biased toward certain tetranucleotides. By calculating the frequencies of every possible tetranucleotide, one can group contigs which likely came from the same organism, even though the aligner failed to align them. Once binned in their own file, a repeat alignment run may be able to align these contigs to build longer contigs than those produced in the original run.

So far, the only publicly available script I've found is this one, however I don't know if I can install it on the computer I'm running computation on:

https://github.com/wwood/bioruby-kmer_counter

I've been using the following sources for guidance:

Dick, G. J., Andersson, A. F., Baker, B. J., Simmons, S. L., Thomas, B. C., Yelton, a P., & Banfield, J. F. (2009). Community-wide analysis of microbial genome sequence signatures. Genome biology, 10(8), R85. doi:10.1186/gb-2009-10-8-r85

Lesniewski, R. A., Jain, S., Anantharaman, K., Schloss, P. D., & Dick, G. J. (2012). The metatranscriptome of a deep-sea hydrothermal plume is dominated by water column methanotrophs and lithotrophs. The ISME journal, 6(12), 2257–68. doi:10.1038/ismej.2012.63

[gringer]

I remembered a Rosalind problem about this:

http://rosalind.info/problems/kmer/

My "quickie" code to solve that problem was written in Perl:

#!/usr/bin/perl
use warnings;
use strict;
my %merCount = ();
my @merOrder = ();
for(my $i = 0; $i < 256; $i++){
    my $mer4 = sprintf("%s%s%s%s",
		       int($i / 64)%4, int($i / 16)%4, 
		       int($i /  4)%4, int($i /  1)%4);
    $mer4 =~ tr/0123/ACGT/;
    $merCount{$mer4} = 0;
		$merOrder[$i] = $mer4;
}
open(my $f, "< input.txt");
my $fastaLabel = <$f>;
my $seq = "";
while(<$f>){
		s/\s+//;
		$seq .= $_;
}
for(my $i = 0; $i < (length($seq)-3); $i++){
		$merCount{substr($seq, $i, 4)}++;
}
my @counts = ();
for(my $i = 0; $i < 256; $i++){
		printf("%s - %d\n",$merOrder[$i],$merCount{$merOrder[$i]});
		push(@counts, $merCount{$merOrder[$i]});
}
print(join(" ", @counts)."\n");

If you don't care about the order and don't want zero counts, then this code simplifies quite a lot.

[LeightonP]

If you can get your sequences into R, then the oligonucleotideFrequency function (which underlies alphabetFrequency/dinucleotideFrequency etc.) may be helpful to you.

http://svitsrv25.epfl.ch/R-doc/library/Biostrings/html/alphabetFrequency.html

oligonucleotideFrequency(yeast1, 4)

where yeast1 is your DNAString/DNAStringSet.

[yzzhang]

the script "calc.kmerfreq.pl" in the pipeline multi-metagenome assembly: https://github.com/MadsAlbertsen/multi-metagenome can do this job

[AndrewRGross]

Awesome: each of you have been a big help.

My current plan is to employ the R script, because the server I'm working on already has R and the associated command line tools installed.

That Albertsen metagenome guide looks very, very promising too, so I might incorporate whatever tools I can from that.

Again, thank you all.