Computational Genomics with R on Your Own Data

Woolf Software
r bioconductor genomics rna-seq variant-annotation
A night forest river splitting into many colored glowing channels that align beneath reeds whose translucent pods flare gold where the light crosses them.

If you already have your own whole-genome VCF, an RNA-seq quantification, a proteomics matrix, and a few months of glucose traces, R is the right tool for roughly the second half of the work. That half covers annotating variants against transcript models and population databases. It also covers normalizing and modeling expression, testing associations, and joining everything onto a common coordinate system or timeline. R is the wrong tool for the first half. Alignment belongs to BWA-MEM2 or minimap2, variant calling to DeepVariant or GATK, and transcript quantification to salmon or kallisto. All of those run outside R.

The existing free textbook on computational genomics with R teaches the language and the statistics well. What it does not do is walk through a single person’s files, and that is where the practical problems live. The sections below follow those files in order. They cover what to run where, how to prepare a VCF, and how to annotate it. They then cover what can and cannot be done with one RNA-seq sample, and how to check that any of it is real.

What R is genuinely good at here, and what to run elsewhere

Before setting up a pipeline it helps to know which parts of it R will make easier and which parts it will make slower. The answer turns mostly on data volume and on Bioconductor’s data model.

R’s advantage is that data model. GenomicRanges provides an interval algebra with sequence names, strand, and genome build attached to every range, so an overlap between your variants and a set of promoter annotations either works or fails loudly on a build mismatch. VariantAnnotation, rtracklayer, Rsamtools, and GenomicFeatures all speak that model, which means a variant set, a bigWig of conservation scores, a BED file of regulatory regions, and a GTF-derived transcript database can be joined without writing coordinate-handling code yourself. That is worth a great deal when you are working alone and every silent off-by-one is a wrong answer you will never notice.

R is poor at anything that streams tens of gigabytes of reads. A 30x whole genome, meaning every position is covered by about thirty reads on average, is roughly 100 GB as a BAM or 50–60 GB as a CRAM. Per-read work in R runs one to two orders of magnitude slower than the equivalent C or Rust tool. The same applies to large-scale statistical genomics: association testing across millions of samples or millions of features has moved to GPU-backed implementations, where the reported speedups over CPU code are in the tens to hundreds of times 1. The practical arrangement is to do the heavy passes with samtools, bcftools, bedtools, and salmon, write intermediate files, and let R start from summaries.

Preparing a personal VCF before R touches it

A VCF, or Variant Call Format file, lists the positions where your genome differs from the reference along with the evidence for each call. Before R reads one, it is worth spending a few minutes putting it into a canonical form, because the failures that follow from skipping this step are quiet rather than loud.

Normalize first. A single-sample whole-genome VCF from DeepVariant has on the order of four to five million variants, and multiallelic records and non-left-aligned indels will silently break joins against gnomAD or ClinVar. The pass we would run is this:

bcftools norm -m -any -f GRCh38.primary_assembly.genome.fa \
  -Ou sample.vcf.gz \
  | bcftools norm -d exact -Oz -o sample.norm.vcf.gz
bcftools index -t sample.norm.vcf.gz

Each flag does one job. -m -any splits multiallelic records into one record per ALT allele, -f left-aligns indels and checks the reference base, and -d exact drops exact duplicates. Read the stderr output afterwards: a nonzero count of “reference allele mismatch” means your VCF and FASTA disagree on build or contig naming, usually chr1 versus 1. Fix that immediately rather than after four hours of annotation.

With a normalized file, read selectively rather than wholesale. readVcf on a full genome VCF will consume well over 20 GB of RAM because it materializes INFO and FORMAT fields as full matrices. Restricting both the fields and the region keeps that under control:

library(VariantAnnotation)
tab <- TabixFile("sample.norm.vcf.gz", yieldSize = 500000)
param <- ScanVcfParam(
  fixed = c("ALT", "FILTER"),
  info  = c("AF"),
  geno  = c("GT", "DP", "GQ")
)
open(tab)
while (nrow(vcf <- readVcf(tab, "hg38", param = param))) {
  # process chunk, keep only what you need
}
close(tab)

The same discipline applies to what each iteration returns. Keep per-chunk output small, ideally a GRanges of the variants that pass your filters plus a handful of annotation columns. On a laptop this brings peak memory down to a few gigabytes.

Annotation: transcript consequences, frequency, and constraint

With a filtered variant set in hand, three annotation layers carry most of the information. They are what a variant does to a transcript, how common it is in the population, and how conserved its position is across species.

The first layer is transcript consequence. GenomicFeatures together with a TxDb package gives you the gene model, and predictCoding returns reference and variant amino acids for coding changes:

library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(BSgenome.Hsapiens.UCSC.hg38)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
loc  <- locateVariants(vcf, txdb, AllVariants())
cod  <- predictCoding(vcf, txdb, BSgenome.Hsapiens.UCSC.hg38)

One caveat matters here. This gives you consequences against the UCSC known-gene set, which is not identical to the Ensembl or RefSeq transcript sets used by VEP and ClinVar. For anything you intend to read carefully, run VEP or bcftools csq with an explicit transcript set and import the annotated VCF back into R. Disagreements between annotators are common, and they usually come down to which transcript is treated as canonical.

The second layer is population frequency, which is the single most effective filter for a personal genome. Most of your four million variants are common and uninteresting. Restricting to a gnomAD allele frequency below roughly 0.1 percent in your inferred ancestry, and to predicted loss-of-function or missense consequences, takes the list from millions to hundreds. This is the standard interpretive funnel, and its limits are well described. Even in a healthy genome, the number of candidate variants far exceeds the number that can be interpreted with confidence, and most rare variants remain of uncertain significance 2. Statistical and machine-learning predictors add ranking rather than certainty 3.

The third layer is evolutionary constraint, imported as a bigWig track with rtracklayer and restricted to your intervals so that you never load a genome-wide track into memory:

library(rtracklayer)
gerp <- import("gerp_hg38.bw", which = my_variants, as = "NumericList")

GERP++ estimates the number of substitutions rejected at each position relative to a neutral rate. Its authors estimated that a large fraction of the human genome, well beyond coding exons, shows measurable constraint 4. A high constraint score on a noncoding variant is a reason to keep it in a shortlist, though it is not a reason to believe it does anything in you.

It bears stating plainly that nothing produced by this pipeline is a clinical result. A research-grade VCF has per-site error rates high enough that any variant you would act on needs orthogonal confirmation in a clinical laboratory and interpretation by a genetics clinician. Use R to build the question and leave the answer to that process.

RNA-seq and the n = 1 problem

Expression data raises a different difficulty. The import is straightforward, and the statistics are where a single person’s data runs into a wall.

Start by importing salmon output with tximport, which handles transcript-to-gene summarization and the effective-length offsets that DESeq2 expects:

library(tximport); library(DESeq2)
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM")

Here the constraint bites. With one sample and one time point there is no dispersion to estimate and therefore no differential expression to test. Two approaches work, depending on what you have.

  • A longitudinal series from the same person. Build a DESeqDataSet with time as the design and use vst(dds, blind = FALSE) for visualization, treating within-person variation as your error term.
  • A single sample. Do not test at all. Rank your genes instead and compare those ranks against a reference distribution such as GTEx whole blood, looking at pathway-level aggregates rather than individual genes.

The reason for that second recommendation is that absolute TPM values in one sample are dominated by blood cell composition, time of day, and library preparation. Any single gene’s value therefore carries little information on its own.

Sanity checks, negative controls, and where R ends

Every enrichment you compute on your own genome needs a null to compare against, and a few identity checks will catch the errors that no statistical test will. This section covers both, and then the point at which R stops being the right environment.

If you ask whether your rare variants fall disproportionately in enhancers, the comparison should be against matched variants rather than uniformly sampled positions. Match them on allele frequency, GC content, and repeat context. Purpose-built synthetic sequence with matched composition and repeat structure exists precisely because naive nulls inflate apparent signal in computational genomics 5. In practice we generate matched control sets with regioneR::randomizeRegions constrained to a mappability mask, and use regioneR::permTest for the p-value.

Two further checks are cheap and catch real errors. Compare your VCF’s genotype calls at a few hundred common SNPs against an independent source, such as an array file or pileups from the BAM via Rsamtools::pileup, to confirm sample identity. Then verify that inferred sex chromosome coverage matches expectation, which catches file mix-ups faster than anything else.

When the R version of an analysis takes more than a few minutes per iteration, move it. data.table will carry you a long way for tabular joins, and plyranges keeps the interval work readable. Once your matrices exceed memory, arrow with Parquet is a better intermediate format than RDS. Beyond that, the work belongs in the command-line tools that were built for read-level throughput and short-read structural inference 6.

Questions people also ask

Do I need the Computational Genomics with R book to do this? It is a good grounding in R programming, statistics, and the standard genomics workflows. It is also free. Because it is organized around methods rather than around one person’s files, pair it with the Bioconductor vignettes for VariantAnnotation, GenomicRanges, and DESeq2, which are where the specifics live.

Can I do variant calling in R? Practically, no. VariantTools and gmapR exist and work for small targeted regions, but for a whole genome you want DeepVariant or GATK HaplotypeCaller. Start from the VCF.

How much RAM do I need? For per-chromosome VCF processing with restricted ScanVcfParam fields, 16 GB is workable and 32 GB is comfortable. A full-genome readVcf with all INFO and FORMAT fields is the usual cause of an out-of-memory crash.

R or Python? Use R for variant annotation, interval arithmetic, and count-based statistics. Bioconductor’s coverage of genome annotation resources has no real equivalent. Use Python for pipeline orchestration, machine learning, and anything GPU-accelerated 1. Mixing the two through Parquet files is less painful than either reticulate or rpy2.

What is the single most useful first analysis? Build one annotated table of your rare, predicted-deleterious variants with gene, consequence, gnomAD frequency, and constraint score, sorted and under a thousand rows. It becomes the object that every later question refers back to.

Woolf Software builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. Amaro Taylor-Weiner, François Aguet, Nicholas J. Haradhvala, et al. Scaling computational genomics to millions of individuals with GPUs. Genome Biology, 2019. https://doi.org/10.1186/s13059-019-1836-7 ↩ ↩2

  2. Tuuli Lappalainen, Alexandra J. Scott, Margot Brandt, et al. Genomic Analysis in the Age of Human Genome Sequencing. Cell, 2019. https://doi.org/10.1016/j.cell.2019.02.032 ↩

  3. Chengsheng Zhu, Maximilian Miller, Zishuo Zeng, et al. Computational Approaches for Unraveling the Effects of Variation in the Human Genome and Microbiome. Annual Review of Biomedical Data Science, 2020. https://doi.org/10.1146/annurev-biodatasci-030320-041014 ↩

  4. Eugene V. Davydov, David L. Goode, Marina Sirota, et al. Identifying a High Fraction of the Human Genome to be under Selective Constraint Using GERP++. PLoS Computational Biology, 2010. https://doi.org/10.1371/journal.pcbi.1001025 ↩

  5. Juan Caballero, Arian F. A. Smit, Leroy Hood, et al. Realistic artificial DNA sequences as negative controls for computational genomics. Nucleic Acids Research, 2014. https://doi.org/10.1093/nar/gku356 ↩

  6. Paolo Carnevali, Jonathan Baccash, Aaron L. Halpern, et al. Computational Techniques for Human Genome Resequencing Using Mated Gapped Reads. Journal of Computational Biology, 2011. https://doi.org/10.1089/cmb.2011.0201 ↩