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

<pre><code>library(BSgenome.Hsapiens.UCSC.hg38)
library(Biostrings)

genome &lt;- BSgenome.Hsapiens.UCSC.hg38
kmers &lt;- oligonucleotideFrequency(genome$chr1, 4)
kmers

m &lt;- as.matrix(kmers)

m[order(m),]
</code></pre>

CC BY-SA 3.0
2017-11-15T12:23:23.280Z
2017-11-15T12:45:32.263Z
2017-11-15T12:45:32.263Z
benn
<p>Grabbing chromosomes for <code>hg38</code>:</p>

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

<p>Via <a href="https://github.com/alexpreynolds/kmer-counter" rel="nofollow noreferrer"><code>kmer-counter</code></a> and Python, here's how to search for kmers of length 7 from chromosome <code>chrY</code>:</p>

<pre><code>#!/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))
</code></pre>

<p>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 <code>chrY</code> (including zero-counts):</p>

<pre><code>CGACGCG 20
CGTCGCG 20
CGCGATA 23
TACGCGC 25
...
AAATAAA 33521
TTTCTTT 34014
GAATGGA 35361
AATGGAA 36906
TTTTTTT 103093
</code></pre>

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

CC BY-SA 3.0
2017-11-15T21:57:50.937Z
2017-11-15T21:57:50.937Z
Alex Reynolds
<p>You can use the <a href="https://www.ncbi.nlm.nih.gov/pubmed/21217122/" rel="noreferrer" title="Link to the article">Jellyfish</a> <a href="http://www.cbcb.umd.edu/software/jellyfish" rel="noreferrer" title="Link to the help page">software</a> to calculate the k-mer profiles up to length 31.</p>

<p>From the instructions in the user <a href="http://www.genome.umd.edu/docs/JellyfishUserGuide.pdf" rel="noreferrer">guide</a>: </p>

<p>The basic command to count all k-mers is as follows:</p>

<pre><code>jellyfish count -m 21 -s 100M -t 10 -C reads.fasta
</code></pre>

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

<pre><code>jellyfish histo mer_counts.jf
</code></pre>

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

<pre><code>jellyfish query mer_counts.jf AACGTTG
</code></pre>

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

<pre><code>jellyfish dump mer_counts.jf &gt; mer_counts_dumps.fa
</code></pre>

CC BY-SA 3.0
2017-11-15T12:55:20.883Z
I am trying jellyfish now.
719016
CC BY-SA 3.0
2017-11-15T16:06:07.653Z
I tried jellyfish, and got the list of kmers. But I couldn't get the list of kmers of count=0.
719016
CC BY-SA 3.0
2017-11-15T11:56:29.283Z
2017-11-15T16:21:05.233Z
2017-11-15T16:21:05.233Z
Bioathlete
<p>What are the least present <em>k</em>-mers in the human genome at different sizes?</p>

<p>Starting with k=4 and going up in size until k=10, what are the <em>k</em>-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.</p>

<p>If this is not pre-calculated somewhere, what tool is recommended to use that starts with the <em><a href="https://en.wikipedia.org/wiki/Reference_genome#Human_reference_genome" rel="noreferrer">GRCh38</a></em> reference as input?</p>

CC BY-SA 3.0
2017-11-16T16:09:53.623Z
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.: http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0029344. Another phrase to google is "word avoidance".
Karel Břinda
<human-genome><sequence-analysis>
Least present k-mers in the human genome
# Least present k-mers in the human genome

(Asked by: username_0 on 2017-11-15)

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](https://en.wikipedia.org/wiki/Reference_genome#Human_reference_genome)* reference as input?

**Tags:** human-genome, sequence-analysis
> **Comments:**
> **username_1**: 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.: http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0029344. Another phrase to google is "word avoidance".  

---

## 3 answers

### Answer by username_4 on 2017-11-15. Score: 12 | Accepted answer

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

> **Comments:**
> **username_0**: I am trying jellyfish now.  
> **username_0**: I tried jellyfish, and got the list of kmers. But I couldn't get the list of kmers of count=0.  

---

### Answer by username_2 on 2017-11-15. Score: 10

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 by username_3 on 2017-11-15. Score: 4

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.


---

