This is a test version of Biostars. For the public version, visit https://www.biostars.org.
kmer counting for heterozygosity estimation

Hi all,

I want to count kmers in my sequencing reads inorder to be able to estimate heterozygosity of my genome using Genomescope. I have paired end reads (R1.fastq, R2.fastq), I ran Jellyfish to count kmers using the following settings:

jellyfish count -C -m 21 -s 5000000 -t 8 R1.fastq -o reads.jf

but open further thinking, I realized I should maybe incorporate R2.fastq as well. so I ran the command as such:

jellyfish count -C -m 21 -s 5000000 -t 8 R*.fastq -o reads.jf

this works, but the resulting heterozygosity values from Genomescope differ. I was wondering if anyone had some input on the right way to count kmers in paired end sequencing data for the purpose of downstream heterozygosity estimations.

Any ideas are greatly appreciated. Thank you!

jellyfish kmer heterozygosity genome assembly
• 2,561 views • link
updated by Dattatray Mongad Dattatray Mongad • written by el97004 • el97004

1 answer

First of all, if you investigate the quality of your forward and reverse reads, you will find the poor quality of reverse reads as compared to forward. So, when you calculate k-mers frequencies for both reads it will be different (and being just imaginary if you cluster forward and reverse reads on basis of k-mer frequencies, the mates will fall apart from each other due to poor quality in reverse reads). I suggest, if its supports your aim you can merge the pairs (assemble) first and then calculate k-mer frequencies. I had used compSeq from emboss toolkit for calculation of k-mers.

• 0 views • link
written by Dattatray Mongad Dattatray Mongad

Log in to answer this question.

More posts like this
  • jellyfish histo is empty
    written by m.t.lorenc • m.t.lorenc
    I ran Jellyfish ``` > zcat V*_R1.fastp.fq.gz | jellyfish count -t 2 -C -m 21 -s 5G --quality-start=33 -o reads.jf > cat reads.jf 000001535{"alignment":8,"canonical":true,"cmdline":["count","-t","2","-C","-m","21","-s","5G","--quality-start=33","-o","reads.jf"],"counter_len":4,"exe_path":"/mnt/hpccs01/work/x/miniconda2/envs/jellyfish/bin/jellyfish","format":"binary/sorted","hostname":"cl5n014","key_len":42,"matrix1":{"c":42,"columns":[7114389000,2192333556,7375600293,3957177567,1753639420,5791419737,7660672264,3955008356,2531272710,4865888693,6095688361,6974692332,2128614574,2258809230,885814170,1703365744,5945426963,375951490,7217126865,6145592432,7418978358,2071372677,3810558711,7150470371,5456096994,5475478,7949231332,5100019741,6318414384,7141154869,3543742516,6024410566,5605745954,575704984,8124532068,1053004640,2675922405,1674046033,5535382637,1853889752,1425936082,3907867046],"identity":false,"r":33},"max_reprobe":126,"pwd":"/mnt/hpccs01/scratch/x/y/fastpDNA","reprobes":[1,1,3,6,10,15,21,28,36,45,55,66,78,91,105,120,136,153,171,190,210,231,253,276,300,325,351,378,406,435,465,496,528,561,595,630,666,703,741,780,820,861,903,946,990,1035,1081,1128,1176,1225,1275,1326,1378,1431,1485,1540,1596,1653,1711,1770,1830,1891,1953,2016,2080,2145,2211,2278,2346,2415,2485,2556,2628,2701,2775,2850,2926,3003,3081,3160,3240,3321,3403,3486,3570,3655,3741,3828,3916,4005,4095,4186,4278,4371,4465,4560,4656,4753,4851,4950,5050,5151,5253,5356,5460,5565,5671,5778,5886,5995,6105,6216,6328,6441,6555,6670,6786,6903,7021,7140,7260,7381,7503,7626,7750,7875,8001],"size":8589934592,"time":"Wed Nov …
  • Jellyfish problem with Failed to open input file 'reads.jf'
    written by m.t.lorenc • m.t.lorenc
    Hi, I can't get jellyfish running with two different ways: Method 1 jellyfish count -t 10 -C -m 21 -s 10G --quality-start=33 -o reads.jf <(zcat …
  • Questions about using Jellyfish and genomescope with paired end reads?
    written by asilidae • asilidae
    I recently used jellyfish with some paired-end reads to count k-mers, and I saw another post that mentioned you should concatenate the paired-end reads before …
  • Jellyfish for Paired End Reads
    written by hpalk42 • hpalk42
    Hello, I am assembling a whole genome for my masters project. I am trying to run Jellyfish and GenomeScope to get some information about genome …
  • Genome Size Estimation with Jellyfish and Genome Scope is Unexpectedly Small
    written by dk0319 • dk0319
    I used jellyfish to count kmers and the obtained histograms were analyzed using Genome scope (http://qb.cshl.edu/genomescope/) in order to estimate effective genome size. However the …
  • Jellyfish histo produces a blank file
    written by mln.mrt • mln.mrt
    I am trying to estimate genome size of a species with high coverage Illumina data with jellyfish using the following arguments: jellyfish count -t 16 …
  • Heterozygosity in k-mer histogram
    written by jm440 • jm440
    Hi! I ran GenomeScope to try to estimate the level of heterozygosity in my genome, however, the output plot looks quite strange and most alarming, …
  • Does Jellyfish support zipped fastq sequences?
    written by comma1993 • comma1993
    I used jellyfish to counting kmers and my input files were zipped fastq. I used the commond : jellyfish count -m 17 -o 17mer_counts -s …
  • How do I count the number of nucleotides using jellyfish; can I obtain the number of nucleotides th…
    written by Tom • Tom
    I have a comparison that's been puzzling me. Suppose I have a fastq file for a bacteria: rawread.fastq I do: `cat rawread.fastq | awk 'NR%4==2{print}' …
  • Quake Error Corrector
    written by Seigfried Seigfried
    Hello I have been trying to run Quake Error Corrector on my Illumina reads before assembly. I am following the manual at [http://www.cbcb.umd.edu/software/quake/manual.html][1]. I have …