Least present k-mers in the human genome

Question

What are the least present k-mers in the human genome at different sizes?

Starting with k=4 and going up in size until k=10, what are the k-mers least seen (or not at all) in the human genome? I am only interested in the reference human genome, so I am discounting SNPs/Indels in the population.

If this is not pre-calculated somewhere, what tool is recommended to use that starts with the GRCh38 reference as input?

Answer 1

You can use the Jellyfish software to calculate the k-mer profiles up to length 31.

From the instructions in the user guide:

The basic command to count all k-mers is as follows:

jellyfish count -m 21 -s 100M -t 10 -C reads.fasta

To compute the histogram of the k-mer occurrences, use the histo subcommand (see section 3.1):

jellyfish histo mer_counts.jf

To query the counts of a particular k-mer, use the query subcommand (see section 3.3):

jellyfish query mer_counts.jf AACGTTG

To output all the counts for all the k-mers in the le, use the dump subcommand (see section 3.2):

jellyfish dump mer_counts.jf > mer_counts_dumps.fa

Answer 2

You can also use R. I give you an example of only chr1 and only kmer=4.

library(BSgenome.Hsapiens.UCSC.hg38)
library(Biostrings)

genome <- BSgenome.Hsapiens.UCSC.hg38
kmers <- oligonucleotideFrequency(genome$chr1, 4)
kmers

m <- as.matrix(kmers)

m[order(m),]

Answer 3

Grabbing chromosomes for hg38:

$ wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/*.fa.gz
$ for fn in `ls *.fa.gz`; do gunzip $fn; done

Via kmer-counter and Python, here's how to search for kmers of length 7 from chromosome chrY:

#!/usr/bin/env python                                                                                                                                                                                                                                                                                                                                                                                                                    

import sys
import subprocess
import itertools

k = 7
chr = 'chrY'
fastaFile = '%s.fa' % (chr)
kmerCmd = 'kmer-counter --fasta --no-rc --k=%d %s' % (k, fastaFile)

try:
    output = subprocess.check_output(kmerCmd, shell=True)
    result = {}
    for line in output.splitlines(): 
        (header, counts) = line.strip().split('\t')
        header = header[1:]
        kmers = dict((key,int(val)) for (key,val) in [d.split(':') for d in counts.split(' ')])
        result[header] = kmers
except subprocess.CalledProcessError as error:
    sys.stderr.write("%s\n" % (str(error)))

kmers = result[chr]
comp = {'A':'T', 'C':'G', 'G':'C', 'T':'A'}
for kmerList in itertools.product('ACGT', repeat=k):
    kmerKey = ''.join(kmerList)
    kmerCompKey = ''.join(reversed([comp.get(b,b) for b in kmerList]))
    if kmerKey not in kmers and kmerCompKey not in kmers:
        kmers[kmerKey] = 0

for key, val in sorted(kmers.iteritems(), key=lambda (key,val):(val,key)):
    sys.stdout.write("%s\t%s\n" % (key, val))

This script will print a two-column tab-delimited text file to standard output, where the first column is the 7mer and the second column is the count of that 7mer and its reverse complement over chrY (including zero-counts):

CGACGCG 20
CGTCGCG 20
CGCGATA 23
TACGCGC 25
...
AAATAAA 33521
TTTCTTT 34014
GAATGGA 35361
AATGGAA 36906
TTTTTTT 103093

With a couple tweaks, this script could run over a range of k and chromosomes for your reference genome.