text: user180: 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?
  
  user776: 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.
  
  user939: 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),]
  
  
  user964: 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
  
  
original_text: 180: <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>
   776: <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>
   939: <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>
   964: <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>
  