Skip to content

Extracting variants

Often you just want a few SNPs for every person: to check a genotype, fit a quick model in R, or make a plot. c4 dosages and get_dosages() pull them out of a cohort's VCFs and link each sample to its patient.

Example output

The demo cohort's VCFs are empty placeholders, so the dosage values on this page come from canned bcftools output. Everything else is what c4 prints.

From the command line

c4 dosages --dataset demo --variants 6:1000:G>A --variants 6:5:A:C -o snps.tsv

Samples:  3
Variants: 1 of 2 extracted
  ! 6:1000:G>A: Alleles flipped relative to VCF (VCF: A>G); dosage inverted.
  ✗ 6:5:A:C: Variant not found in VCF
✓ Dosages saved: snps.tsv (+ snps.variants.tsv)

snps.tsv has one row per VCF sample and one column per variant:

patient_id  imputed_id  6:1000:G>A  6:5:A:C
P001        ARR-0001    2.0         NA
P003        ARR-0003    1.0         NA
NA          CTRL-1      1.0         NA

snps.variants.tsv says how each variant was matched:

variant     vcf_id         status  warning                                                          error
6:1000:G>A  chr6:1000:A:G  ok      Alleles flipped relative to VCF (VCF: A>G); dosage inverted.     NA
6:5:A:C     NA             error   NA                                                               Variant not found in VCF

Options:

Option Meaning
--variants, -v A variant; repeat for several. Quote anything with >.
--variant-list A file with one variant per line.
--variant-build Build of your variants, e.g. b37. They are lifted over to the VCFs' build with a chain from the catalog.
--dna Which genetic data set to read (a dna: entry, default imputed).
--no-dosage Use hard calls (GT) instead of the DS dosage field.
--keep-dropped Keep samples of patients marked drop=1 in the ID map.
--chain-file Use this liftover chain instead of the catalog's.

In Python

from c4.genotypes import get_dosages, write_dosages

dos, info = get_dosages("demo", ["6:1000:G>A", "6:5:A:C"])
write_dosages(dos, info, "snps.tsv")

dos and info are the two tables above. Merge dos with anything else you loaded on patient_id:

from c4.cohort import cohort

demo = cohort("demo", show=False)
df = demo.load("pheno").merge(dos, on="patient_id", how="left")

How variants are matched

  • Formats: chr:pos:ref:alt, chr:pos:ref>alt, or chr:pos for whatever is at that position. The chr prefix is optional.
  • The dosage counts copies of the ALT allele as you wrote it. If the VCF has the alleles the other way round, the dosage is flipped (2 − DS) and the variant gets a warning saying so.
  • Alleles that match neither way give an error and an empty column.
  • Multi-allelic positions: the VCF record with your alleles is used. If none matches, or you gave no alleles, the first record is used, with a warning.
  • Another genome build: with --variant-build b37 (or build="b37"), positions are lifted over first, in Python with pyliftover (no R needed). The chain comes from the catalog's resources.
    • Positions are 1-based throughout, as in VCFs and dbSNP.
    • Where the chain maps a position onto the reverse strand of the new build, the alleles are complemented (A↔T, C↔G) before matching, and the variant gets a warning saying so. Indels there are reported as errors, since their position would also shift.
    • A position that falls in a gap of the chain, maps to more than one place, or maps to another chromosome is an error for that variant; the others still run.
  • Samples: every sample in the VCF, linked to patient_id through the ID map. Samples that aren't patients (controls) are kept with an empty patient_id. Patients marked drop=1 are left out unless you keep them.

Saving the dosages behind survival and PRS results

The survival and PRS commands use exactly the same extraction. Add --save-dosages to keep what they used:

c4 survival --dataset demo -v 6:1000:A:G --survival OS --save-dosages surv_dosages.tsv
c4 prs --dataset demo -w weights.tsv --survival OS -o prs_out --save-dosages

survival writes the file you name (plus .variants.tsv); prs writes prs_dosages.tsv and prs_dosages.variants.tsv into its output directory. For PRS, the dosages are coded to the ALT allele in the weights file's variant IDs, before flipping to the effect allele.