DOCSEARCHbyAlgolia
Submit Docs

GenomeScope 2.0

repository·master·Indexed 41 days ago

https://github.com/tbenavi1/genomescope2.0

A tool for reference-free genome profiling of polyploid genomes (up to hexaploid). It uses a polyploid-aware mixture model to estimate genome size, heterozygosity, and repetitiveness from k-mer frequency distributions generated by KMC or Jellyfish.

Tokens
2.5K
Snippets
5
Records
11
Agent score
50%

What's inside GenomeScope 2.0

  1. Check ploidy limits and estimation capabilities

    master

    GenomeScope 2.0 is designed for reference-free profiling of polyploid genomes with the following constraints:

    • Ploidy Range: Supports genomes up to hexaploid (6) ploidy levels.
    • Ploidy Estimation: GenomeScope does not estimate ploidy itself; use Smudgeplot for this task.
    • Aneuploidy/Unequal Chromosomes: GenomeScope does not support genomes with uneven copy numbers of chromosomes (e.g., aneuploid cancer genomes or unequal sex chromosomes). In these cases, the reported heterozygosity rate will reflect the fraction of bases that are haploid (copy number 1) versus diploid (copy number 2), alongside any heterozygous positions in other chromosomes.

  2. Understand GenomeScope genome size estimation logic

    master

    GenomeScope estimates haploid genome size by dividing the total number of k-mers (excluding presumptive sequencing errors) by the estimated homozygous coverage ($2\lambda$).

    Key Variables and Concepts

    • $\lambda$ (lambda) and kcov: These refer to the same value in the code and documentation. It represents the estimated coverage for heterozygous k-mers.
    • Homozygous Coverage: The model identifies peaks at $1\lambda$, $2\lambda$, $3\lambda$, and $4\lambda$. The $2\lambda$ peak corresponds to the mean coverage of unique homozygous sequences.
    • Haploid vs. Diploid:
      • Haploid genome size: The span of one complete set of chromosomes.
      • Diploid genome size: The total DNA content in one diploid cell (roughly $2 \times$ haploid size).
    • Scaling: K-mers with higher coverage (e.g., $200\times$ or $300\times$) are automatically scaled by the homozygous coverage ($2\lambda$) to contribute the correct number of copies to the estimate.
    • Unequal Haplotypes: If the two haplotypes have significantly different lengths, the reported haploid genome size is the average of the two.

  3. Manage high-frequency k-mers and contamination

    master

    To handle artifacts or highly repetitive genomes, use the following strategies:

    Filtering Contamination

    Use the --max_kmercov parameter to set a cutoff for high-frequency k-mers. This is useful for excluding artifacts like phiX sequencing or organelle/contamination sequences that may appear as extremely high-frequency k-mers (e.g., $10,000\times$ or $100,000\times$).

    Handling Repetitive/Large Genomes

    If you are working with highly repetitive or very large genomes, the default jellyfish histogram (which truncates at $10,000$) might be insufficient. You may need to:

    1. Regenerate the histogram from jellyfish or KMC with a higher upper bound (e.g., $100,000$, $1,000,000$, or higher).
    2. Note that while this includes ultra-high frequency k-mers in the analysis (potentially increasing the estimated genome size), it also increases sensitivity to data artifacts.

  4. Prepare k-mer frequency histograms using Jellyfish

    master

    To use Jellyfish for k-mer profiling, follow these two steps:

    1. Count k-mers: Use jellyfish count. Always use the -C flag for "canonical k-mers" to account for both DNA strands. Adjust -m (k-mer length) and -s (size/memory) according to your data.
    2. Generate histogram: Use jellyfish histo to redirect the output to a .histo file.

    Requirements: High-quality sequencing (e.g., Illumina) with at least 15x haploid genome coverage is recommended.

    # 1. Count k-mers
jellyfish count -C -m 21 -s 1000000000 -t 10 *.fastq -o reads.jf

# 2. Generate histogram
jellyfish histo -t 10 reads.jf > reads.histo

  5. Prepare k-mer frequency histograms using KMC

    master

    Before using GenomeScope 2.0, you must generate a k-mer count histogram from your sequencing reads. Using KMC involves two steps: counting the k-mers and then transforming the output into a histogram.

    1. Count k-mers: Use kmc to generate a k-mer database. Adjust -m (memory) and -t (threads) based on your hardware. Use -k for k-mer length (21 is recommended for most genomes). -ci and -cs set the lower and upper count bounds.
    2. Transform to histogram: Use kmc_tools transform to create the .histo file. The -cx flag sets the upper cutoff for the histogram.

    Requirements: High-quality sequencing (e.g., Illumina) with at least 15x haploid genome coverage is recommended. Long-read data (Nanopore/PacBio) is currently not supported due to high error rates.

    # 1. Count k-mers
mkdir tmp
ls *.fastq > FILES
kmc -k21 -t10 -m64 -ci1 -cs10000 @FILES reads tmp/

# 2. Transform to histogram
kmc_tools transform reads histogram reads.histo -cx10000

  6. Install GenomeScope 2.0 via Conda and R

    master

    To run GenomeScope 2.0 on the command line, you need an R environment with specific packages. Follow these steps:

    1. Install R dependencies: Use conda to install r-base, r-minpack.lm, and r-argparse.
    2. Clone the repository: Download the source from GitHub.
    3. Install the R package: Run the provided install.R script within the repository.
    4. Setup PATH: Ensure genomescope.R is in your system PATH to run it from any directory.
    # Install R dependencies
conda install -c conda-forge r-base r-minpack.lm r-argparse

# Clone and enter directory
git clone https://github.com/tbenavi1/genomescope2.0.git
cd genomescope2.0/

# Install the package
Rscript install.R

  8. Prepare k-mer histograms for GenomeScope Online

    master

    GenomeScope Online requires a k-mer count histogram file as input. You can generate this using either Jellyfish or KMC.

    Using Jellyfish

    1. Count k-mers: Use jellyfish count. Ensure you use the -C flag for canonical k-mers.
      $ jellyfish count -C -m 21 -s 1000000000 -t 10 *.fastq -o reads.jf
      Note: Adjust -s (memory) and -t (threads) to your server capacity. -m is the k-mer length.
    2. Export histogram: Use jellyfish histo to create the file for upload.
      $ jellyfish histo -t 10 reads.jf > reads.histo

    Using KMC

    1. Count k-mers: Use kmc with a file containing your input fastq paths.
      $ mkdir tmp
$ ls *.fastq > FILES
$ kmc -k21 -t10 -m64 -ci1 -cs10000 @FILES reads tmp/
      Note: -k is k-mer length, -ci and -cs are lower and upper count bounds.
    2. Export histogram: Use kmc_tools transform to generate the histogram.
      $ kmc_tools transform reads histogram reads.histo -cx10000
      Note: -cx sets the upper bound cutoff for the histogram.

    Once you have the .histo file, upload it to the GenomeScope web interface.

    # Jellyfish workflow
jellyfish count -C -m 21 -s 1000000000 -t 10 *.fastq -o reads.jf
jellyfish histo -t 10 reads.jf > reads.histo

# KMC workflow
kmc -k21 -t10 -m64 -ci1 -cs10000 @FILES reads tmp/
kmc_tools transform reads histogram reads.histo -cx10000

  9. Troubleshoot model convergence and unexpected results

    master

    If the GenomeScope model fails to converge or produces results that deviate significantly from expectations (e.g., lack of peaks in k-mer plots or a poor model fit), check the following:

    1. Coverage: Ensure you have sufficient sequencing coverage. Low coverage prevents the model from confidently identifying homozygous k-mer peaks.
    2. K-mer Counting Mode: Verify that you used the canonical k-mer counting mode. If using jellyfish, use the -C flag.
    3. K-mer Size: If coverage is sufficient but the model still fails, try decreasing the k-mer size to 19 or 17.
    4. Data Quality: If the above steps fail, additional sequencing data may be required.

  10. Run GenomeScope 2.0 via command line

    master

    Once installed, use the genomescope.R script to perform genome profiling. The script requires an input histogram, an output directory, and the k-mer length used during the counting step.

    Required Arguments:

    • -i: Path to the input histogram file (generated by KMC or Jellyfish).
    • -o: Directory where output plots and results will be saved.
    • -k: The k-mer length used to generate the histogram.

    Optional Arguments:

    • -p ploidy: Sets the ploidy of the model.
    • -l lambda: Initial guess for the average k-mer coverage.
    • -n 'name_prefix': Prefix for output files.
    • -m max_kmercov: Cutoff for excluding high-frequency k-mers.
    # Basic usage
genscope.R -i histogram_file -o output_dir -k k-mer_length

# Example with Arabidopsis F1 data
genscope.R -i ara_F1_21.hist -o output -k 21

  11. Configure GenomeScope analysis parameters

    master

    When uploading a histogram to GenomeScope, you can provide the following parameters to refine the polyploid genome profiling:

    ParameterDescription
    descriptionA label for your sample (default: my sample)
    kmer_lengthThe length of the k-mers used in the histogram (e.g., 21)
    ploidyThe expected ploidy of the genome (range: 1 to 6, default: 2)
    max_kmer_covMaximum k-mer coverage to filter out high copy-number DNA (e.g., chloroplasts). Set to -1 for no filter.
    lambdaAverage k-mer coverage for polyploid genome. Set to -1 for default.