# Least present k-mers in the human genome

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?

• If you are interested in least present k-mers, you might be interested also in the concept of minimum absent words (these two things are related). There are many papers on this topic, see, e.g.: journals.plos.org/plosone/article?id=10.1371/…. Another phrase to google is "word avoidance". Nov 16, 2017 at 16:09

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

• I am trying jellyfish now. Nov 15, 2017 at 12:55
• I tried jellyfish, and got the list of kmers. But I couldn't get the list of kmers of count=0. Nov 15, 2017 at 16:06

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),]  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():
kmers = dict((key,int(val)) for (key,val) in [d.split(':') for d in counts.split(' ')])
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.