8
$\begingroup$

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?

$\endgroup$
1
  • 1
    $\begingroup$ 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". $\endgroup$ Nov 16, 2017 at 16:09

3 Answers 3

12
$\begingroup$

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
$\endgroup$
2
  • $\begingroup$ I am trying jellyfish now. $\endgroup$
    – 719016
    Nov 15, 2017 at 12:55
  • $\begingroup$ I tried jellyfish, and got the list of kmers. But I couldn't get the list of kmers of count=0. $\endgroup$
    – 719016
    Nov 15, 2017 at 16:06
10
$\begingroup$

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),]
$\endgroup$
4
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.