Why OncoSeq

A cancer genomics study rarely stops at a list of differentially expressed genes. The same tumours carry mutations, copy-number changes and clinical follow-up, and the interesting questions sit between those layers: which transcriptional programme accompanies a driver mutation, whether a copy-number gain actually raises expression, whether a signature predicts outcome once stage is accounted for.

OncoSeq keeps all of that in one object. OncoSeqExperiment extends SummarizedExperiment, so everything you already know about Bioconductor containers applies, and adds three layers plus a record of every analysis you run.

library(OncoSeq)
library(SummarizedExperiment)

Getting data in

From your own quantification

importCounts() takes a count matrix and a table of sample annotation, either as objects or as paths to delimited files, and checks that their sample identifiers agree.

ose <- importCounts("counts.tsv", "samples.tsv",
                    mutations = "cohort.maf")

importTximport() summarises salmon or kallisto output to the gene level, and asOncoSeq() wraps a SummarizedExperiment, DESeqDataSet or DGEList you already have.

From public repositories

importGDC() pulls a TCGA project from the Genomic Data Commons, assembling expression, clinical data and, on request, mutations and copy number into one object.

blca <- importGDC("TCGA-BLCA", sampleType = "Primary Tumor",
                  mutations = TRUE, cnv = TRUE)

importCBioPortal() does the same against cBioPortal’s curated studies, which is quicker and needs no download cache.

cbioStudies("bladder")[, c("studyId", "name")]
blca <- importCBioPortal("blca_tcga_pan_can_atlas_2018",
                         genes = c("TP53", "RB1", "FGFR3", "ERBB2"))

The example objects

Two datasets ship with the package. exampleOncoSeq() simulates a small count-scale cohort, which is what the differential-expression sections below use. blcaExample is a real 420-gene by 400-sample slice of the TCGA bladder cohort, used for the genomic and survival sections.

ose <- exampleOncoSeq()
ose
#> class: OncoSeqExperiment
#> dim: 300 features x 40 samples
#> assays(1): counts
#> colData(5): condition, age, stage, os_time, os_event
#> rowData(0):
#> mutations: 42 variants in 8 genes across 20 samples
#> cnv: 10 genes x 40 samples
#> analyses: none

Quality control

qcMetrics() writes per-sample depth and detection statistics into colData(), where they are available as covariates or filters.

ose <- qcMetrics(ose)
head(colData(ose)[, c("library_size", "detected_genes")])
#> DataFrame with 6 rows and 2 columns
#>     library_size detected_genes
#>        <numeric>      <numeric>
#> T01        95420            297
#> T02        93531            299
#> T03        98575            298
#> T04        92543            297
#> T05        94196            299
#> T06        91894            298

detectOutliers() flags samples that sit far from the bulk of the data, both by distance in principal-component space and by low correlation with the other samples. It reports; it does not remove.

out <- detectOutliers(ose)
out[out$outlier, ]
#> DataFrame with 1 row and 5 columns
#>     pc_distance median_correlation distance_outlier correlation_outlier
#>       <numeric>          <numeric>        <logical>           <logical>
#> N17     8.81888           0.869381            FALSE                TRUE
#>       outlier
#>     <logical>
#> N17      TRUE

filterLowExpression() drops features too weakly expressed to test. With a grouping variable it uses edgeR::filterByExpr().

ose <- filterLowExpression(ose, group = "condition")

Normalisation

ose <- normalizeCounts(ose, method = "vst")
assayNames(ose)
#> [1] "counts" "vst"

The raw counts are never overwritten; each method adds an assay. A principal-component plot is the quickest check that the groups separate.

plotSamplePCA(ose, colourBy = "condition")

Samples in the space of the first two principal components.

Differential expression

One interface covers three engines. The contrast is written the same way for all of them.

ose <- diffExpress(ose, ~ condition,
                   contrast = c("condition", "tumour", "normal"),
                   engine = "DESeq2")
head(getAnalysis(ose, "de"))
#>   feature   baseMean    log2FC     lfcSE       stat       pvalue         padj
#> 1 GENE215  582.67470 -3.077663 0.1937232 -15.886909 7.808401e-57 2.217586e-54
#> 2 GENE031  938.02804 -3.005172 0.1928173 -15.585591 9.121583e-55 1.295265e-52
#> 3 GENE237 1189.00215  2.278510 0.2272104  10.028194 1.145929e-23 1.084813e-21
#> 4 GENE160   76.36814 -2.048323 0.2144237  -9.552688 1.263744e-21 8.972580e-20
#> 5 GENE020   55.99033 -2.495868 0.2681566  -9.307501 1.308756e-20 7.433736e-19
#> 6 GENE104  314.90242 -1.438728 0.1649423  -8.722611 2.718588e-18 1.286798e-16
#>   engine
#> 1 DESeq2
#> 2 DESeq2
#> 3 DESeq2
#> 4 DESeq2
#> 5 DESeq2
#> 6 DESeq2

Covariates go in the design and are adjusted for:

ose <- diffExpress(ose, ~ stage + condition,
                   contrast = c("condition", "tumour", "normal"))

Because the three engines make different distributional assumptions, agreement between them is a useful sanity check on a signature.

cmp <- compareEngines(ose, ~ condition,
                      contrast = c("condition", "tumour", "normal"))
cmp$overlap
#>    DESeq2     edgeR     limma consensus 
#>        26        25        24        23

plotVolcano(ose, nLabel = 8)

Volcano plot of the tumour versus normal contrast.

plotExpressionHeatmap(ose, nTopGenes = 30, annotationColumn = "condition")

Standardised expression of the most variable features.

Gene set enrichment

runGSEA() ranks every tested feature and calls fgsea; runORA() tests the significant list against each set with the hypergeometric distribution. Both store a table in the same format.

sets <- readGMT(system.file("extdata", "example_sets.gmt",
                            package = "OncoSeq"))
ose <- runGSEA(ose, sets, minSize = 3)
head(getAnalysis(ose, "gsea")[, c("set", "NES", "padj")])
#>                            set        NES      padj
#> 1    UROTHELIAL_BASAL_SQUAMOUS  1.5940620 0.2177594
#> 2        CHROMATIN_REMODELLING  1.4188844 0.3412017
#> 3          RTK_PI3K_SIGNALLING  1.1014779 0.6268657
#> 4           CELL_CYCLE_CONTROL  1.0272019 0.6268657
#> 5 UROTHELIAL_LUMINAL_PAPILLARY -0.9497487 0.6268657
#> 6          DNA_DAMAGE_RESPONSE -0.8797426 0.6268657

plotEnrichmentBar(ose, "gsea")

Normalised enrichment scores.

For real work, load a collection with readGMT() or, if you have msigdbr installed, msigdbrSets("Homo sapiens", "H").

The genomic layers

From here on we use the real bladder cohort.

data(blcaExample)
blcaExample
#> class: OncoSeqExperiment
#> dim: 420 features x 400 samples
#> assays(1): expression
#> colData(18): sampleId, patientId, AGE, SEX, RACE, AJCC_PATHOLOGIC_TUMOR_STAGE ...
#> rowData(0):
#> mutations: 14337 variants in 420 genes across 400 samples
#> cnv: 420 genes x 400 samples
#> analyses: none

Mutations

head(mutationFrequency(blcaExample, topN = 10))
#>     gene n_mutated n_samples frequency  ci_lower  ci_upper
#> 1   TP53       199       400    0.4975 0.4474271 0.5476101
#> 2    TTN       188       400    0.4700 0.4202191 0.5202279
#> 3  MUC16       118       400    0.2950 0.2507200 0.3423439
#> 4  KMT2D       113       400    0.2825 0.2388814 0.3293707
#> 5 ARID1A       105       400    0.2625 0.2200275 0.3085265
#> 6  KDM6A       105       400    0.2625 0.2200275 0.3085265

plotMutationFrequency(blcaExample, topN = 15)

Mutation frequency with exact binomial confidence intervals.

plotOncoprint(blcaExample, topN = 10)

Mutation landscape of the ten most frequently altered genes.

Tumour mutation burden is written back into colData():

blcaExample <- tumourMutationBurden(blcaExample, captureSize = 38)
summary(colData(blcaExample)$tmb)
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#>  0.05263  0.39474  0.71053  0.94322  1.18421 12.15789

Note that the burden computed here reflects only the genes in this trimmed panel, not the whole exome; on a complete MAF it is the conventional non-synonymous-variants-per-megabase figure.

Copy number

head(cnvSummary(blcaExample)$perGene)
#>     gene n_gain n_loss n_samples freq_gain freq_loss
#> 1 CDKN2A      1    129       400    0.0025    0.3225
#> 2   E2F3     62      1       400    0.1550    0.0025
#> 3 VPS13B     54      0       400    0.1350    0.0000
#> 4  PPARG     53      0       400    0.1325    0.0000
#> 5   UBR5     51      0       400    0.1275    0.0000
#> 6  CCND1     43      0       400    0.1075    0.0000

plotCNVFrequency(blcaExample, topN = 15)

Gain and loss frequency across the cohort.

Whether a copy-number change actually moves the transcript is a separate question:

head(cnvExpressionCorrelation(blcaExample), 5)
#>     gene   n  estimate       pvalue         padj
#> 1 CDKN2A 400 0.7470393 1.423398e-72 5.978273e-70
#> 2   E2F3 400 0.6557226 1.610271e-50 3.381570e-48
#> 3   ICE1 400 0.6303721 1.073342e-45 1.502679e-43
#> 4  PSIP1 400 0.6248941 1.035275e-44 1.087039e-42
#> 5  VIRMA 400 0.6209017 5.250229e-44 4.410192e-42

Survival

survivalFit() stratifies by whatever layer the feature lives in. It searches expression, mutation, copy number and clinical data in that order unless you name the layer.

fit <- survivalFit(blcaExample, "TP53", by = "mutation")
fit
#> Kaplan-Meier analysis of TP53 (mutation) 
#> groups: wild-type (n=201), mutant (n=199) 
#> log-rank p: 0.703 
#> hazard ratio: 0.944
plotKM(fit, timeLabel = "Months")

Overall survival by TP53 mutation status.

For a continuous feature, choose how to split. The "optimal" rule searches for the most separating threshold, which inflates significance; the result reports both the raw and the search-corrected p-value, and only the latter should be quoted.

survivalFit(blcaExample, "MKI67", by = "expression", split = "optimal")

A multivariable model mixes clinical covariates with expression (expr:) and mutation (mut:) terms:

coxModel(blcaExample, c("AGE", "expr:TP53", "mut:RB1"), name = NULL)
#>            term       coef       HR  ci_lower ci_upper          se         z
#> 1           AGE 0.03102412 1.031510 1.0156721 1.047596 0.007894821 3.9296793
#> 2     expr.TP53 0.01649459 1.016631 0.8794485 1.175213 0.073957917 0.2230267
#> 3 mut.RB1mutant 0.14819206 1.159736 0.8066775 1.667317 0.185219425 0.8000892
#>         pvalue
#> 1 8.505923e-05
#> 2 8.235147e-01
#> 3 4.236591e-01

survivalScreen() fits one Cox model per gene and corrects across genes.

screen <- survivalScreen(blcaExample, name = NULL)
head(screen, 5)
#>       gene       HR ci_lower ci_upper        z       pvalue   n      padj
#> 1 ADAMTS12 1.371737 1.172387 1.604984 3.944961 7.981285e-05 400 0.0159883
#> 2     TCHH 1.342803 1.158546 1.556364 3.914255 9.068394e-05 400 0.0159883
#> 3     LRP1 1.349185 1.158220 1.571636 3.846317 1.199065e-04 400 0.0159883
#> 4    AHNAK 1.378670 1.167490 1.628049 3.785463 1.534228e-04 400 0.0159883
#> 5   AHNAK2 1.360811 1.157500 1.599833 3.731508 1.903369e-04 400 0.0159883

Integrating layers

The question that motivates the whole container: what does a driver mutation do to the transcriptome?

colData(blcaExample)$TP53_status <- mutationStatus(blcaExample, "TP53")
res <- diffExpress(blcaExample, ~ TP53_status,
                   contrast = c("TP53_status", "mutant", "wild-type"),
                   engine = "limma", assay = "expression", voom = FALSE,
                   name = "de_TP53")
head(getAnalysis(res, "de_TP53"), 5)
#>   feature baseMean     log2FC lfcSE      stat       pvalue         padj engine
#> 1  CDKN2A 1229.132  1228.6262    NA  6.844353 2.903085e-11 1.219296e-08  limma
#> 2   CENPF 1927.917   724.2118    NA  6.283604 8.659861e-10 1.268080e-07  limma
#> 3  NUP188 1615.998   462.9302    NA  6.275930 9.057716e-10 1.268080e-07  limma
#> 4    MDM2 3914.314 -4516.8990    NA -6.065656 3.050089e-09 3.202594e-07  limma
#> 5  CDKN1A 5476.273 -2181.0172    NA -5.717367 2.120827e-08 1.781495e-06  limma

On a count-scale object, integrateMutationExpression() does the same in one call and handles the group-size checks for you.

Running everything at once

runOncoSeq() chains quality control, filtering, normalisation, differential expression, enrichment, mutation burden and a survival screen, skipping any step whose inputs are missing.

full <- runOncoSeq(exampleOncoSeq(), ~ condition,
                   contrast = c("condition", "tumour", "normal"),
                   geneSets = sets, engine = "edgeR",
                   normalisation = "vst", verbose = FALSE)
names(analyses(full))
#> [1] "de"              "gsea"            "ora"             "survival_screen"

Every step is recorded, so an object can account for itself:

vapply(provenance(full)$steps, `[[`, character(1), "step")
#> [1] "qcMetrics"            "filterLowExpression"  "normalizeCounts"     
#> [4] "diffExpress"          "runGSEA"              "runORA"              
#> [7] "tumourMutationBurden" "survivalScreen"       "runOncoSeq"

Session information

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] OncoSeq_0.99.1              SummarizedExperiment_1.43.0
#>  [3] Biobase_2.73.2              GenomicRanges_1.65.4       
#>  [5] Seqinfo_1.3.2               IRanges_2.47.5             
#>  [7] S4Vectors_0.51.10           BiocGenerics_0.59.12       
#>  [9] generics_0.1.4              MatrixGenerics_1.25.0      
#> [11] matrixStats_1.5.0           BiocStyle_2.41.0           
#> 
#> loaded via a namespace (and not attached):
#>  [1] fastmatch_1.1-8     gtable_0.3.6        xfun_0.62          
#>  [4] bslib_0.12.0        ggplot2_4.0.3       httr2_1.3.0        
#>  [7] ggrepel_0.9.8       lattice_0.23-1      vctrs_0.7.3        
#> [10] tools_4.6.1         parallel_4.6.1      tibble_3.3.1       
#> [13] pkgconfig_2.0.3     Matrix_1.7-6        data.table_1.18.6.1
#> [16] RColorBrewer_1.1-3  S7_0.2.2            lifecycle_1.0.5    
#> [19] compiler_4.6.1      farver_2.1.2        statmod_1.5.2      
#> [22] fgsea_1.39.2        DESeq2_1.53.6       codetools_0.2-20   
#> [25] htmltools_0.5.9     sys_3.4.3           buildtools_1.0.0   
#> [28] sass_0.4.10         yaml_2.3.12         pillar_1.11.1      
#> [31] jquerylib_0.1.4     BiocParallel_1.47.0 limma_3.99.2       
#> [34] DelayedArray_0.39.8 cachem_1.1.0        abind_1.4-8        
#> [37] locfit_1.5-9.12     tidyselect_1.2.1    digest_0.6.39      
#> [40] dplyr_1.2.1         labeling_0.4.3      splines_4.6.1      
#> [43] maketools_1.3.2     cowplot_1.2.0       fastmap_1.2.0      
#> [46] grid_4.6.1          cli_3.6.6           SparseArray_1.13.4 
#> [49] magrittr_2.0.5      S4Arrays_1.13.2     survival_3.8-12    
#> [52] withr_3.0.3         edgeR_4.99.10       scales_1.4.0       
#> [55] rmarkdown_2.32      XVector_0.53.0      otel_0.2.0         
#> [58] evaluate_1.0.5      knitr_1.52          rlang_1.3.0        
#> [61] Rcpp_1.1.2          glue_1.8.1          BiocManager_1.30.27
#> [64] jsonlite_2.0.0      R6_2.6.1

Analysing cancer transcriptomes and genomes with OncoSeq