PaperlayerView exact revision on arXiv

Source-linked AI summary

A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data

Heng Li

arXiv:1203.6372v3q-bio.GN

TL;DR

Uncertain genotypes limit sequence analysis in low-coverage and somatic-mutation settings. The paper develops a framework that analyzes sequencing data directly, achieving comparable accuracy to alternative methods across allele-count estimation, allele-frequency-spectrum inference, and association mapping.

  • Problem

    Low-coverage sequencing may not cover non-reference alleles, making separate per-sample variant calling and genotype-based analysis unreliable.

  • Method

    The framework uses statistical inference to call variants and estimate population parameters or test associations directly from uncertain sequencing data.

  • Results

    The method achieves comparable accuracy to alternative approaches for allele-count estimation, allele-frequency-spectrum inference, and association mapping on real data.

  • Takeaways & Limitations

    Many population-genetic parameters and statistical tests can be adapted to uncertain sequences, so exact genotypes are not always indispensable.

  • Takeaways & Limitations

    A biallelic assumption can cause false positives when triallelic variants are mismodeled.

Abstract

from arXiv · show

Motivation: Most existing methods for DNA sequence analysis rely on accurate sequences or genotypes. However, in applications of the next-generation sequencing (NGS), accurate genotypes may not be easily obtained (e.g. multi-sample low-coverage sequencing or somatic mutation discovery). These applications press for the development of new methods for analyzing sequence data with uncertainty. Results: We present a statistical framework for calling SNPs, discovering somatic mutations, inferring population genetical parameters and performing association tests directly based on sequencing data without explicit genotyping or linkage-based imputation. On real data, we demonstrate that our method achieves comparable accuracy to alternative methods for estimating site allele count, for inferring allele frequency spectrum and for association mapping. We also highlight the necessity of using symmetric datasets for finding somatic mutations and confirm that for discovering rare events, mismapping is frequently the leading source of errors. Availability: http://samtools.sourceforge.net. Contact: hengli@broadinstitute.org.

Motivation: · 1 INTRODUCTION

Low-coverage, many-sample sequencing improves population sampling but makes variant discovery and accurate individual genotyping difficult. The paper therefore develops methods that compute genetic statistics directly from uncertain sequencing data without precise genotypes or mandatory imputation.

  • 1 INTRODUCTION: At 2–6 fold coverage per sample, non-reference alleles may be missed, making separate variant calls followed by combining calls unreliable.Jointly considering all samples is presented as the preferred strategy for improving variant-discovery power.
  • 1 INTRODUCTION: Joint analysis largely solves variant discovery, but accurate genotypes for each individual remain unavailable for methods such as HWE testing and association mapping.Without accurate genotypes, many methods developed for genotype data cannot be used directly.
  • 1 INTRODUCTION: Genotype imputation reuses genotype-based methods by transferring information across linked loci and between individuals.The motivating example involves low coverage at site A, linked support at site B, and high-coverage samples at A.
  • 1 INTRODUCTION: Imputation cannot infer regional allele frequency spectra, may bias population-genetic inference through LD dependence, and is computationally slow.Its application is limited to candidate variant sites and its effectiveness depends on the pattern of LD.
  • 1 INTRODUCTION: These limitations motivate deriving new methods for medical and population-genetic problems without precise genotypes.The paper explicitly questions whether imputation is always preferable when its main purpose is to reuse genotype-based methods.
  • 1 INTRODUCTION: Somatic and germline mutation discovery also requires genotype information, even when related samples are sequenced at high coverage.Mutation discovery is challenging because sequencing errors can exceed the much lower mutation rate.
  • 1 INTRODUCTION: The paper computes various statistics directly from sequencing data without knowing genotypes and evaluates the methods on real data.This stated scope addresses the preceding motivation for genotype-free analysis.

2 METHODS

The Methods section formalizes equations for estimating genotype frequency and allele frequency spectrum and performing Hardy–Weinberg equilibrium and association tests, while defining common notation and reviewing existing approaches.

  • Methods: The section presents equations for inferring genotype frequency and allele frequency spectrum and conducting Hardy–Weinberg equilibrium and association tests.The equations use the paper’s notation and include formulations drawn from existing literature for theoretical completeness.
  • Methods: The final subsection reviews existing methods and summarizes differences among them and between their formulations and the proposed framework.Detailed derivations are provided in an online document referenced by the paper.
  • Methods: The framework models n individuals with individual-specific ploidy m_i, sequence data d_i, and genotype g_i as the number of reference alleles in [0, m_i].These definitions establish the common notation used throughout the Methods section.
  • Methods: The article considers only biallelic variants.This restriction is stated as an explicit Methods assumption.

2.1 Assumptions

The framework assumes independence across sites, reads, and individuals, while treating all variants as biallelic. These assumptions can be violated by sequence context and alignment dependencies, although sample independence is essential to the derivations.

  • Assumptions: Data at different sites are assumed independent, although INDEL errors or polymorphisms can create nearby alignment correlations.Post-filtering and explicit modeling of mapping and alignment errors may reduce the effects of site dependency.
  • Assumptions: Sequencing and mapping errors among reads at a site are assumed independent, yielding independent likelihood functions across individuals.
  • Assumptions: Error dependence on sequence context can violate independence, but sample independence is essential to all subsequent derivations.Error dependency within an individual can be modeled, whereas the sample-independence assumption cannot be relaxed in the presented derivations.
  • Assumptions: All variants are assumed biallelic; triallelic SNPs comprise about 0.2% of the human population.This assumption has little impact on SNP modeling but may matter more for INDELs at microsatellites.

2.2 Computing genotype likelihoods

Genotype likelihoods are computed from sequencing-read bases and their base qualities, focusing on the two most evident nucleotide types at each site under a biallelic-variant model.

  • Data representation: For each sample and site, sequencing data comprise read bases and their base qualities.The framework represents the data as an array of bases on sequencing reads plus their base qualities.
  • Biallelic-site simplification: For biallelic variants, the method retains the two most evident nucleotide types and drops less evident types when present.At most two nucleotide types are therefore considered at any site; this simplification is described as sufficient in practice, though not optimal.
  • Likelihood calculation: Assuming independent sequencing errors, genotype likelihoods are derived from the observed read bases and their error probabilities.For k reads, the first l may match the reference and the remaining reads may differ, with ε_j denoting the j-th read-base error probability.

2.3 Inferences from multiple samples

The framework infers allele, genotype, haplotype, and allele-count distributions directly from sequencing data, while providing tests for HWE and group differences and Bayesian methods for variant calling and AFS estimation. It also addresses computational limitations and robustness concerns, including HWE violation, haplotype-scaling costs, and numerical underflow.

  • Allele-frequency inference: Allele-frequency inference uses maximum likelihood under HWE, optimized numerically with Brent’s method or by EM over missing sample genotypes.EM is faster when genotype likelihoods strongly favor one genotype; after 10 EM rounds without convergence, implementation switches to Brent’s method.
  • HWE testing: For diploid samples, HWE testing uses a likelihood-ratio statistic that incorporates genotype uncertainty and approximately follows a 1-degree χ2 distribution.The statistic approaches the standard HWE test statistic for genotype data.
  • Haplotype inference: Haplotype frequencies are estimated iteratively from genotype likelihoods, but joint computation across many loci has O(n · 4^k) complexity and is mainly practical for measuring LD between two loci.When genotypes are certain, the procedure reduces to the standard EM method for genotype data.
  • Group-difference testing: Group-difference testing compares allele-frequency estimates between two sample groups, but HWE violation can cause false positives and the safer alternative is underpowered when both groups satisfy HWE.The proposed practice is to use the P-value from Da1 while filtering candidates with low Da2.
  • Allele-count inference: The framework computes a discrete likelihood for site reference allele count, which is more convenient to maximize or sum over than the allele-frequency likelihood and supports Bayesian inference.An efficient recursive computation is used to obtain the allele-count likelihood from sequencing data.
  • Variant calling and AFS: Bayesian variant calling places a prior on the sample allele-frequency spectrum, defines variant quality from the allele-count likelihood, and is applicable to pooled resequencing with unequal pool sizes.AFS can also be estimated by EM across sites (EM-AFS) or by a normalized histogram of per-site maximum-likelihood allele counts (site-AFS).

2.4 Discovering somatic and germline mutations

The section presents likelihood-ratio frameworks for discovering somatic, de novo, and germline mutations directly from sequencing data, emphasizing genotype support and protection against false positives from triallelic variants.

  • Normal–tumor pairs: Somatic mutations in normal–tumor pairs can be identified using a likelihood ratio computed from the samples’ genotype likelihoods.A candidate is considered certain only when both samples’ genotypes are clearly better supported than alternative genotypes.
  • Family trios: The trio extension compares maximum likelihoods with and without Mendelian inheritance constraints, with larger ratios indicating greater mutation confidence.The constrained configurations satisfy Mendelian inheritance and use child, father, and mother genotype likelihoods.
  • Multiallelic variants: The implementation allows multiallelic variants because assuming biallelic variation can produce false somatic or germline mutations through triallelic-variant mismodeling.The authors specifically observed false calls caused by mismodeling triallelic variants.

2.5 Working with diploid multi-allelic sites

The framework models diploid multi-allelic sites through allele frequencies under Hardy–Weinberg equilibrium, using EM estimation and likelihood-ratio tests for multi-allelic status and variant calling.

  • 2.5 Working with diploid multi-allelic sites: At a site with p alleles, allele frequencies ψ_h sum to 1, and Hardy–Weinberg equilibrium defines the likelihood for estimating them.The model parameterizes the site by allele frequencies {ψ_1, …, ψ_p}.
  • 2.5 Working with diploid multi-allelic sites: Diploid genotypes are represented as unordered allele pairs ⟨h, h′⟩, from which the EM iteration equation is derived.The unordered-pair representation corresponds to diploid genotypes.
  • 2.5 Working with diploid multi-allelic sites: Multi-allelic status is tested with a likelihood ratio that approximately follows a 1-degree χ2 distribution.The same framework also supports variant calling through a corresponding hypothesis test, as proposed by Kim et al. (2010).

2.6 Related works

Related methods estimate site allele frequency by numerical likelihood maximization or EM, while this framework adds procedures for HWE, haplotype frequency, association testing, somatic mutation discovery, allele-count estimation, and AFS inference.

  • Related methods: Thunder and glfMultiples estimate site allele frequency by numerically maximizing the likelihood, whereas GATK uses EM and Kim et al. use both algorithms.Li et al. (2010b) proposed an alternative estimator not covered in the article.
  • Contributions: The framework tests Hardy–Weinberg equilibrium and estimates haplotype frequency from sequencing data.
  • Contributions: It proposes a two-degree association test and a simple but effective model for discovering somatic mutations.
  • Contributions: It rigorously derives and numerically stabilizes a discrete allele count estimator and introduces an EM algorithm for inferring the allele frequency spectrum.

3 RESULTS

The framework is implemented in SAMtools and GATK and performs competitively across allele-count, allele-frequency-spectrum, association-testing, and somatic-difference analyses. Results show that performance depends on linkage, coverage, HWE control, and mapping-error mitigation.

  • Implementation: Most diploid-sample equations are implemented in SAMtools, while exact Eqs. (17)–(19) are the default SNP-calling model in GATK.SAMtools is freely available under the MIT open-source license.
  • Allele-count estimation: Allele-count RMSD was 15.0 without imputation versus 12.7 with imputation in 670 European samples, whereas target resequencing favored imputation-free estimation, 42.5 versus 54.8.The authors attribute the target-resequencing result to a lack of linked SNPs around fragmented regions and caution that Omni sites may not represent all SNPs.
  • Allele frequency spectrum: The EM-AFS method inferred a fairly accurate allele frequency spectrum with 3-fold coverage per sample, while site-AFS estimates were less stable despite the correct overall trend.Accuracy was evaluated by comparing low-coverage 1000 Genomes data with high-coverage Complete Genomics data.
  • Association testing: The Da1 association statistic largely followed the 1-degree χ2 distribution but produced one false positive SNP at P < 10^-6 when HWE was violated.Using Da2 suppressed the false positive to P > 0.001, so HWE control is important for the 1-degree likelihood-ratio test.
  • Somatic mutation discovery: Mapping-error mitigation reduced 667 apparent differences between HiSeq and older Illumina data to 33 after remapping and 13 after clustered-SNP filtering.The remaining 13 were potential differences, while SomaticSniper identified 1,826 differences at a comparable score cutoff.

4 DISCUSSIONS

The framework analyzes sequencing data directly, supporting mutation discovery, population-genetic estimation, and association testing without explicit SNP calling or imputation. Evaluations show comparable accuracy to alternative methods and emphasize when uncertain sequence data are preferable to inferred genotypes.

  • The framework supports SNP and mutation discovery, allele-frequency and linkage-disequilibrium estimation, Hardy–Weinberg testing, and association analysis without explicitly calling SNPs or genotypes.It can discover somatic and germline mutations with appropriate input data.
  • On real data, the method achieved comparable accuracy to the best alternative methods and produced sensible results on several unpublished datasets.The evaluation supports obtaining useful information directly from sequencing data without SNP calling or imputation.
  • Imputation can transfer genotyping knowledge to low-coverage sequencing data and may be more accurate with homogeneous whole-genome data from many samples.Its accuracy depends on nearby linkage disequilibrium, a relationship supported here by direct evidence from real data.
  • When sequencing or mapping makes sequences uncertain, directly analyzing uncertain sequence data may preserve more information than using arbitrarily ascertained sequences.Exact DNA sequences or genotypes are often needed only as intermediate results for parameter estimation or statistical computation.
Loading 1203.6372v3…