Analyzing Your Own Molecular Data in R

Oak
r bioconductor rnaseq genomics reproducibility
Five colored glowing streams braid into one pool in a living stone basin, with an identical basin repeating behind it in misty forest.

By the end of this guide you will have a single R project that loads five kinds of data about one person and holds them in objects you can query: a VRanges of your filtered germline variants annotated against a transcript database, a SummarizedExperiment of longitudinal RNA-seq counts imported from salmon, a log2 protein abundance matrix with missingness handled explicitly, and a tidy time-indexed table joining blood biomarkers to continuous glucose data. You will also have an renv.lock file that pins every package version, so the analysis reruns identically in two years. You need R 4.4 or later, a matching Bioconductor release (3.20 or 3.21), roughly 32 GB of RAM if you plan to load a whole-genome VCF in memory, and your own files: a bgzipped and tabix-indexed VCF, one salmon or kallisto quantification directory per RNA-seq timepoint, a proteomics panel export (Olink NPX or SomaScan RFU, usually wide CSV), lab results as CSV, and a CGM export as CSV. Everything below is measurement and interpretation. Nothing here tells you what any result means for your health, and several steps end at a point where you should take the output to a clinician.

On the R-versus-Python question that brings most people to this topic: we use both, and the split is not arbitrary. Read alignment, variant calling, and quantification are command-line tools (bwa-mem2, DRAGEN, GATK, salmon) orchestrated by Snakemake or Nextflow. Everything downstream of a count matrix or a variant table, meaning normalization, differential testing, correlation networks, and figures, is where R earns its place, because Bioconductor’s core data structures encode genomic coordinates and sample metadata together and the statistical methods you want (limma’s empirical Bayes moderation, DESeq2’s shrinkage estimators) have their reference implementations there. Use Python for glue, machine learning, and anything touching raw sequence at scale.

1. Set up a project that reruns

Start with dependency pinning, because Bioconductor packages depend on the R version in a way that breaks silently. Bioconductor releases twice a year and each release is valid for exactly one R minor version, so a script written against Bioconductor 3.18 and run under R 4.5 will install different package versions without warning you.

install.packages("renv")
renv::init(bioconductor = "3.21")

install.packages("BiocManager")
BiocManager::install(c(
  "VariantAnnotation", "GenomicRanges", "GenomicFeatures",
  "TxDb.Hsapiens.UCSC.hg38.knownGene", "org.Hs.eg.db",
  "tximport", "DESeq2", "limma", "WGCNA", "SummarizedExperiment",
  "AnnotationHub"
))
install.packages(c("data.table", "dplyr", "lubridate", "iglu", "targets"))

renv::snapshot()
BiocManager::valid()   # flags packages out of sync with the release

Commit renv.lock alongside your scripts. If you expect to hand the project to someone else or to a compute environment you do not control, build on top of the Bioconductor Docker image (bioconductor/bioconductor_docker:RELEASE_3_21) and mount your data directory read-only, which removes the system library dependencies (libxml2, GSL, HDF5) that cause most first-run install failures. Containerized, version-pinned pipelines are the practical answer to reproducibility failures in bioinformatics, and the tooling for it is mature enough that there is no reason to skip it 1.

Keep raw data read-only and write every derived object to a derived/ directory as Parquet or RDS with a timestamped filename. You will rerun step 4 a dozen times and you want to know which version of the count matrix produced which figure.

2. Load and filter your variants

A whole-genome VCF for one person is typically 4 to 8 million variant records and 1 to 3 GB compressed. Do not read it whole unless you have the memory. Read by region with ScanVcfParam, or stream it in chunks.

library(VariantAnnotation)

vcf_file <- "raw/sample.dragen.hard-filtered.vcf.gz"
tbx <- TabixFile(vcf_file, yieldSize = 100000)

# Region-restricted read: a single gene's locus
roi <- GRanges("chr17", IRanges(43044295, 43170245))  # BRCA1, hg38
param <- ScanVcfParam(which = roi,
                      info = c("AC", "AF"),
                      geno = c("GT", "DP", "GQ"))
vcf <- readVcf(vcf_file, genome = "hg38", param = param)

Filter on the caller’s own quality fields before you do anything interpretive. For a germline single-sample callset we keep records where FILTER == "PASS", genotype quality GQ >= 20 (a phred-scaled 99% confidence in the genotype call), and read depth DP >= 10. Depth matters more than people expect: at DP of 8, a true heterozygote has roughly a 1-in-100 chance of being called homozygous by sampling alone, and low-depth regions cluster in exactly the GC-rich promoters and segmental duplications you might care about.

gt <- geno(vcf)$GT[, 1]
dp <- geno(vcf)$DP[, 1]
gq <- geno(vcf)$GQ[, 1]

keep <- fixed(vcf)$FILTER == "PASS" & !is.na(dp) & dp >= 10 & gq >= 20
vcf  <- vcf[keep, ]

Then annotate consequence against a transcript database. locateVariants assigns each variant to a region class (coding, intron, promoter, five/threeUTR, intergenic) and predictCoding gives you the amino acid change for coding variants.

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, seqSource = Hsapiens)
table(cod$CONSEQUENCE)   # synonymous / nonsynonymous / nonsense / frameshift

predictCoding will emit warnings for variants whose reference allele does not match the genome. That is a real signal, almost always meaning your VCF is on GRCh37 while your TxDb is hg38, or the seqlevels are mismatched. Fix it with seqlevelsStyle(vcf) <- "UCSC" for the naming problem and with rtracklayer::liftOver plus a UCSC chain file for the build problem, checking how many variants fail to map (a few percent is normal, twenty percent means something is wrong).

One thing to be clear about: annotating a variant as nonsynonymous, or finding it in ClinVar with a pathogenic assertion, is not a clinical result. Germline variant interpretation depends on allele frequency in matched ancestry, functional evidence, segregation, and the specific ACMG criteria a laboratory director applies. If a variant in your file looks consequential, the next step is a clinical-grade confirmatory test ordered through a physician or genetic counselor, not a downstream R script.

3. Ask interval questions with GenomicRanges

Most of the useful genomic queries on a personal VCF are interval arithmetic: which of my variants fall in a liver enhancer, in a pharmacogene, in a region a published GWAS implicated. GenomicRanges makes this one function call, and AnnotationHub gives you the interval sets.

library(AnnotationHub)
ah <- AnnotationHub()
query(ah, c("EncodeDCC", "hg38", "cCREs"))   # browse, then pull by ID
ccre <- ah[["AH104048"]]                      # example accession

hits <- findOverlaps(rowRanges(vcf), ccre)
length(unique(queryHits(hits)))

The failure mode here is silent rather than loud: if one object uses chr1 and the other uses 1, findOverlaps returns zero hits and no error. Always check intersect(seqlevels(a), seqlevels(b)) before trusting an overlap count of zero, and standardize with seqlevelsStyle(). The second trap is coordinate convention. VCF and GRanges are 1-based inclusive, BED is 0-based half-open, and rtracklayer::import handles the conversion while a naive fread of a BED file does not.

If you want to go further and score whether a noncoding variant of yours sits in a functional element, two R-based approaches are worth knowing. monaLisa tests transcription factor motifs for enrichment across binned sequence sets, which is how you would ask whether the regions containing your variants are enriched for a particular factor’s binding motif 2. gkmSVM trains a gapped k-mer support vector machine on regulatory sequence and produces deltaSVM scores for the effect of a single-nucleotide change on predicted regulatory activity 3. Both give you a quantitative prior, not a conclusion.

4. Import RNA-seq across timepoints

Quantify outside R. We use salmon with a decoy-aware index built from the GENCODE transcriptome plus the primary assembly as decoys, and always with --gcBias and --seqBias, because a personal longitudinal series will otherwise pick up library-prep drift as apparent biology.

salmon quant -i gencode_v47_decoy_idx -l A \
  -1 t1_R1.fastq.gz -2 t1_R2.fastq.gz \
  --gcBias --seqBias --validateMappings \
  -p 16 -o quants/t1

Then import with tximport, which aggregates transcript-level estimates to genes and passes along the effective-length offsets DESeq2 needs.

library(tximport); library(DESeq2)

samples <- read.csv("meta/rnaseq_samples.csv")  # id, date, batch, fasted
files <- file.path("quants", samples$id, "quant.sf")
names(files) <- samples$id

tx2gene <- read.csv("ref/tx2gene_v47.csv")      # TXNAME, GENEID
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                ignoreTxVersion = TRUE)

dds <- DESeqDataSetFromTximport(txi, colData = samples,
                                design = ~ batch + fasted)
dds <- dds[rowSums(counts(dds) >= 10) >= 3, ]
dds <- DESeq(dds)

Here is the statistical reality of one subject. DESeq2’s negative binomial model estimates dispersion from replication, and with a single sample per condition there is no replication and the test is undefined. What makes a personal profile analyzable is repeated sampling of the same person: if you have six fasted draws and six fed draws, or eight draws before and eight after some change, the timepoints are your replicates and the design above is legitimate, with the caveat that observations from the same person on nearby days are correlated rather than independent. For that, limma-voom with duplicateCorrelation() is the better tool, because it estimates a within-subject correlation and adjusts the residual degrees of freedom accordingly.

For exploratory work on a series, skip hypothesis testing and use vst(dds, blind = FALSE) to get a homoscedastic matrix, then look at PCA and at the trajectory of specific genes over time. Note that PC1 in a personal series is frequently batch or blood cell composition, not anything about you. Deconvolution of cell-type proportions from whole blood, which is what a standard PAXgene RNA sample measures, explains a large share of variance and should be estimated and adjusted for before you interpret any single gene’s trend.

5. Find modules instead of testing genes one at a time

With a longitudinal series you are better positioned to ask which sets of genes move together over time than to ask whether one gene differs between two draws. Weighted gene correlation network analysis builds a signed co-expression network, picks a soft-thresholding power that makes the network approximately scale-free, and clusters genes into modules summarized by an eigengene, which is the first principal component of the module 4.

library(WGCNA)
enableWGCNAThreads(8)

expr <- t(assay(vst(dds, blind = TRUE)))
expr <- expr[, order(-apply(expr, 2, var))[1:8000]]

sft <- pickSoftThreshold(expr, powerVector = 1:20, networkType = "signed")
net <- blockwiseModules(expr, power = sft$powerEstimate,
                        networkType = "signed", TOMType = "signed",
                        minModuleSize = 30, mergeCutHeight = 0.25,
                        maxBlockSize = 10000, numericLabels = TRUE)
me <- moduleEigengenes(expr, net$colors)$eigengenes

WGCNA needs samples, not genes: below roughly fifteen observations the correlation estimates are too noisy and modules are unstable. If you have twenty or more timepoints, correlating module eigengenes against your biomarker panel or your glucose summaries is the most informative thing you can do with a personal RNA-seq series. This module-then-correlate pattern is the standard structure of integrative analyses in the literature, where differential expression is followed by network and pathway-level summarization rather than treated as the endpoint 5.

6. Handle proteomics missingness deliberately

Affinity proteomics arrives as a wide matrix already on a log2-ish scale (Olink NPX) or as linear relative units (SomaScan RFU, which you should log2 transform). Two things matter. First, values below the limit of detection are reported but unreliable, and dropping or imputing them changes your answer, so record which is which in a parallel logical matrix rather than silently filtering. Second, plate and batch effects in a series collected over a year are large enough to dominate biology.

library(limma)

npx <- as.matrix(read.csv("raw/olink_npx_wide.csv", row.names = 1))  # proteins x samples
below_lod <- read.csv("raw/olink_lod_flags.csv", row.names = 1) == TRUE
npx <- npx[rowMeans(below_lod) < 0.2, ]      # keep proteins detected in >80%

npx_adj <- removeBatchEffect(npx, batch = samples$plate)

design <- model.matrix(~ fasted, data = samples)
fit <- lmFit(npx_adj, design)
fit <- eBayes(fit, trend = TRUE)
topTable(fit, coef = "fasted", number = 20)

Use removeBatchEffect for visualization and include the batch term in the design matrix for testing. Doing both double-corrects and produces anticonservative p-values.

7. Join biomarkers and glucose on time

The last object to build is the one that makes the rest interpretable: a table indexed by time that carries lab values and glucose summaries next to the dates of your RNA and protein draws. The problem is almost entirely about time handling.

library(data.table); library(lubridate); library(iglu)

cgm <- fread("raw/dexcom_export.csv")
cgm[, time := ymd_hms(timestamp, tz = "America/New_York")]
setnames(cgm, "glucose_mg_dl", "gl")
cgm[, id := "self"]

daily <- rbind(
  as.data.table(cv_glu(cgm[, .(id, time, gl)])),
  as.data.table(in_range_percent(cgm[, .(id, time, gl)], target_ranges = list(c(70, 140)))),
  fill = TRUE
)

labs <- fread("raw/labs_long.csv")           # date, analyte, value, unit
labs[, date := as_date(date)]

joined <- merge(labs, cgm_daily_summary, by.x = "date", by.y = "date", all.x = TRUE)

Parse every timestamp with an explicit timezone. CGM exports are usually in device local time with no offset recorded, so a trip across timezones or a DST transition creates duplicated or missing hours that will show up as a spurious dip in your daily aggregates. Check cgm[, .N, by = as_date(time)] and confirm that a 5-minute-interval sensor gives you about 288 readings per day, flagging any day below roughly 70% completeness before you compute daily metrics on it. Glucose variability metrics are descriptive summaries of your own data. Whether any pattern in them is meaningful for you is a question for a physician, and abnormal values should go to one.

8. Write the pipeline down and make it machine-readable

Turn the steps above into a targets pipeline so that changing one input reruns only what depends on it, and export the final objects in formats an analysis agent can consume without your code: Parquet for tabular data, plain TSV for gene-level results with unversioned Ensembl IDs, and a sessioninfo::session_info() dump alongside every output. Large language models are now reasonable at writing and reviewing this kind of downstream analysis code and at summarizing results, and the quality of what they produce depends heavily on whether your outputs carry explicit units, identifiers, and provenance 6. They are unreliable narrators of biological significance, so treat generated interpretation as a hypothesis to check against primary literature.

Common problems

Chromosome naming and genome build mismatches cause more wasted hours than any other issue. Any overlap that returns zero hits should be suspected of chr1 versus 1 before it is believed, and any predictCoding run that warns about reference mismatches is telling you the build is wrong.

Gene identifier drift is next. Ensembl IDs carry versions (ENSG00000141510.17) and GENCODE releases change them, so ignoreTxVersion = TRUE in tximport plus a single pinned annotation release for the whole project is the only stable approach. Mapping to gene symbols at the end rather than the beginning avoids the many-to-one collisions that silently drop rows.

Trying to run DESeq2 or limma on two samples produces either an error about degrees of freedom or, worse, numbers that look like results. With one subject, testing requires repeated measures, and even then the within-subject correlation must be modeled.

Memory exhaustion on whole-genome VCFs is avoidable with ScanVcfParam(which = ...) and yieldSize chunking. If you find yourself needing the whole thing in memory, you probably want the variants as a filtered Parquet table queried with arrow or duckdb instead.

Finally, a mismatched R and Bioconductor pair will install older package versions without complaint. Run BiocManager::valid() after every install and keep renv.lock under version control.

Oak 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. Neha Kulkarni, Luca Alessandrì, Riccardo Panero, et al. Reproducible bioinformatics project: a community for reproducible bioinformatics analysis pipelines. BMC Bioinformatics, 2018. https://doi.org/10.1186/s12859-018-2296-x ↩

  2. Dania Machlab, Lukas Burger, Charlotte Soneson, et al. monaLisa: an R/Bioconductor package for identifying regulatory motifs. Bioinformatics, 2022. https://doi.org/10.1093/bioinformatics/btac102 ↩

  3. Mahmoud Ghandi, Morteza Mohammad-Noori, Narges Ghareghani, et al. gkmSVM: an R package for gapped-kmer SVM. Bioinformatics, 2016. https://doi.org/10.1093/bioinformatics/btw203 ↩

  4. Peter Langfelder, Steve Horvath. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics, 2008. https://doi.org/10.1186/1471-2105-9-559 ↩

  5. Siwei Liu, Xie Xie, Huajiang Lei, et al. Identification of Key circRNAs/lncRNAs/miRNAs/mRNAs and Pathways in Preeclampsia Using Bioinformatics Analysis. Medical Science Monitor, 2019. https://doi.org/10.12659/msm.912801 ↩

  6. Zhigang Meng, Zhikai Yang, Mingxun Zhu, et al. Large language models in bioinformatics: a comprehensive survey. Frontiers in Genetics, 2026. https://doi.org/10.3389/fgene.2026.1797863 ↩