--- title: "Introduction to thyroidBRS" author: "Camila Barrios-Cervantes, Daniel Pérez-Calixto, Hugo Tovar" package: thyroidBRS output: BiocStyle::html_document: toc_float: true vignette: > %\VignetteIndexEntry{Introduction to thyroidBRS} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup-knitr, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") ``` ## Introduction Papillary thyroid carcinoma (PTC) is predominantly driven by either a `BRAF`^V600E^ mutation or a RAS-family mutation. Agrawal et al. (2014) (*Cell*, [doi:10.1016/j.cell.2014.09.050](https://doi.org/10.1016/j.cell.2014.09.050)) showed that the two driver groups have distinguishable transcriptional profiles, and defined the **BRAF-RAS Score (BRS)** from a 71-gene signature (their Figure S7A). The study reported a BRS for the 391 tumors that had **both** exome and RNA sequencing — exome data was needed to say which tumors belonged to the BRAF-mutant reference group and which to the RAS-mutant one, but **not** to score a tumor once those groups and their centroids are fixed. Any additional sample with expression data can be scored against the same centroids. `thyroidBRS` implements exactly that two-step idea: 1. **`brs_fit()`** — given a cohort with known driver mutations for at least some samples, computes the two reference centroids. 2. **`predict()`** — scores *any* sample against those centroids, including samples that were never labeled. ## Installation `thyroidBRS` is installed from Bioconductor with `BiocManager`: ```{r install, eval = FALSE} if (!requireNamespace("BiocManager", quietly = TRUE)) { install.packages("BiocManager") } BiocManager::install("thyroidBRS") ``` Then load it: ```{r setup} library(thyroidBRS) ``` ## How the score is defined Extended Experimental Procedures 14.1 defines the score of a tumor `t` as the difference between its distance to the BRAF centroid `c(B)` and to the RAS centroid `c(R)`: $$\mathrm{BRS}(t) = \lVert v(t) - c(B) \rVert_2 - \lVert v(t) - c(R) \rVert_2$$ where `v(t)` is the sample's expression across the signature, z-scored gene-by-gene against the reference cohort, and each norm is divided by the square root of the number of signature genes, so the score does not depend on how many of them survived. Negative scores are BRAF^V600E^-like, positive ones RAS-like; `brs_class` spells the two `"Braf-like"` and `"Ras-like"`. `predict()` returns two flavours of the same quantity: - `brs_score` — the value above. Stable per sample: it does not change when you score a different set of samples against the same centroids. - `brs_scaled` — the same value rescaled so the scored set spans $[-1, +1]$, which is the axis the published figures use. This is a property of the *set*, so it moves if the cohort changes and is undefined for a single sample. ## What this reproduces This package is a re-implementation from the paper's description; the original code was never released. It reproduces the published *classification*, not the published numbers: | Compared against | Result | | --- | --- | | Driver-mutation status, 10-fold CV | 99.7% | | Published class labels, samples held out of the fit | 96.4% | | Published continuous score, Spearman | 0.941 | None of these is a clean generalization estimate. The cross-validation re-fits the centroids in each fold, but the signature does not move: those 71 genes were selected on these same TCGA tumors, so gene selection sits outside the loop and 99.7% is an upper bound. The 96.4% is the cleanest figure, since those tumors carry other drivers and entered neither the signature derivation nor the fit. The 96.4% still measures agreement with the original classifier, not with a biological truth. Expression pipelines have moved since 2014 (RSEM on hg19 then, STAR on GENCODE v36 now) and that difference is not reconstructable, so the scores this returns are close to but not the same as the ones in the 2014 supplementary table. Classify and rank with them; do not expect them to match value for value. ## Which samples define the reference groups This is the easiest thing to get wrong, so it is worth being explicit. In the paper, set B is the **BRAF-V600E-mutant** tumors and set R is the **RAS-mutant** tumors. The groups are defined by genotype. They are *not* the BRAF-like/RAS-like classes the BRS produces. In the published TCGA-THCA tables those are the column `BRAF_RAF_class`, which is exactly `sign(BRAF_RAF_score)` — the classifier's own output. Fitting centroids on it is circular, and it pulls tumors with other drivers into both reference groups: 43 other-driver tumors and 1 RAS-mutant tumor land in the "BRAF" group, and 68 other-driver tumors in the "RAS" group. On TCGA-THCA the two fits happen to end up close — the axes correlate at 0.998 and only 2 of 505 tumors change class — so the practical damage is small there. The problem is what it does to the *evaluation*: if you fit on the published classes and then score the same samples, agreeing with them measures nothing. `brs_fit()` accepts either vocabulary, so nothing breaks if you have `"Braf-like"`/`"Ras-like"` labels, but prefer the mutation columns. ## A worked example `thyroidBRS` ships the signature but no expression data. This example uses a small simulated cohort so it runs quickly; the same calls work unchanged on real data. ```{r simulate-data} set.seed(1) genes <- na.omit(brs_genes$current_symbol) samples <- paste0("sample_", seq_len(60)) # Two groups with a real difference across the signature, so the example # shows the classifier doing something rather than fitting noise. The 24 # samples without a label are a mix of both kinds, as a real cohort is. group <- rep(c("BRAF_V600E", "RAS", "unknown"), c(24, 12, 24)) up_in_ras <- brs_genes$up_in[!is.na(brs_genes$current_symbol)] == "RAS" shift <- ifelse(up_in_ras, 1, -1) side <- ifelse(group == "BRAF_V600E", -0.5, 0.5) side[group == "unknown"] <- rep(c(-0.5, 0.5), length.out = 24) expr <- matrix(rnorm(length(genes) * length(samples)), nrow = length(genes), dimnames = list(genes, samples)) expr <- expr + outer(shift, side) # Mutation labels for the 36 samples that have exome data... labels <- setNames(group[seq_len(36)], samples[seq_len(36)]) # ...and 24 with expression only. unlabeled <- setdiff(samples, names(labels)) ``` ### Fit the reference centroids ```{r fit} fit <- brs_fit(expr, labels) fit ``` ### Score the previously-unlabeled samples This is the point of the package: extending classification beyond the exome-sequenced subset. ```{r classify-new} preds_new <- predict(fit, expr[, unlabeled]) head(preds_new) table(preds_new$brs_class) ``` ### Check against the labels you do have ```{r validate} preds_all <- predict(fit, expr) validate_brs(preds_all, labels, fit = fit) ``` Note the `n_resubstituted` field and the warning. Every one of those 36 samples also defined the centroids, so this number is a resubstitution estimate, not a validation. It is a useful smoke test — if it is *not* high, something is wrong — but it is not evidence that the classifier generalizes. How optimistic is it? At the dimensions of the real TCGA analysis (70 genes, 391 samples split 272/119), resubstitution on data with no signal whatsoever still lands around 68%. For a real estimate, hold samples out: ```{r crossval} set.seed(2) n_folds <- 4 fold <- sample(rep(seq_len(n_folds), length.out = length(labels))) # Fit on all folds but one, score the held-out fold, and name each call by # its sample. cv_fold <- function(k) { f <- brs_fit(expr, labels[fold != k]) held <- names(labels)[fold == k] setNames(predict(f, expr[, held, drop = FALSE])$brs_class, held) } cv <- unlist(lapply(seq_len(n_folds), cv_fold))[names(labels)] truth <- ifelse(labels == "BRAF_V600E", "Braf-like", "Ras-like") mean(cv == truth) ``` ### One-step alternative `brs_score()` wraps `brs_fit()` + `predict()` for the common case of scoring the same cohort you fit on: ```{r one-shot} all_preds <- brs_score(expr, labels) head(all_preds) ``` ## Expression containers `brs_fit()`, `predict()` and `brs_score()` also take a `SummarizedExperiment` or an `ExpressionSet`, so data from `TCGAbiolinks` or `GEOquery` goes in without unwrapping. With the simulated cohort from above: ```{r containers, message = FALSE} library(SummarizedExperiment) # colData() has to carry driver-mutation status. Whatever you name the # column, it must not be the published BRAF-like/RAS-like class: that is # the output of this same classifier, and fitting on it is circular. # Samples without exome data get NA and are simply not used for the fit. se <- SummarizedExperiment( assays = list(logexpr = expr), colData = DataFrame(driver = replace(group, group == "unknown", NA), row.names = samples) ) fit_se <- brs_fit(se, "driver", assay = "logexpr") identical(fit_se, fit) head(predict(fit_se, se, assay = "logexpr")) ``` On a GDC object the calls are the same, with the assay named in both and the transform requested: `brs_fit(se, "driver", assay = "tpm_unstranded", log2_transform = TRUE)`, then `predict(fit, se, assay = "tpm_unstranded")`. Two things to watch. Pick the assay explicitly — the first assay of a GDC object is raw counts, which is not what you want, and `brs_fit()` warns if you leave the choice to it. And `labels` may be given as a single string, in which case it is read as a column of `colData()` / `pData()`. ## Microarray data, and what does not transfer The signature behaves well on microarrays. On GSE33630 (Affymetrix HG-U133 Plus 2.0, 49 papillary carcinomas) all 70 usable genes are present, the first principal component of the signature accounts for 46% of the variance, and 68 of those 70 genes load on it with the sign their Figure S7A block predicts. The biology is intact. What does not survive a platform change is the **zero threshold**. A `brs_fit` stores the per-gene means and SDs of the cohort it was fitted on; feed it data from a different platform and those recentre the samples somewhere else entirely. Scoring GSE33630 against TCGA RNA-seq centroids ranks the tumors essentially perfectly — Spearman 0.97 against the array's own signature axis — and then calls all 49 of them BRAF-like. So: - **Same platform and pipeline as the fit** — everything works, and the frozen standardization is what makes the score robust. Scoring a cohort that is 100% RAS-like still gets 97.5% of them right, because nothing about the score depends on the composition of the set being scored. - **A different platform** — fit new centroids on that platform if you have driver-mutation labels for it. If you do not, use `brs_score` as a *ranking* and do not trust the class call. `predict()` warns when every sample lands on the same side of zero, which is the usual symptom. Re-standardizing against the new cohort instead is tempting and does fix the cross-platform threshold, but it buys that by making every sample's score depend on the company it keeps: applied to a cohort that happens to be all RAS-like, it gets 8.4% of them right. That is why it is not the default. ## Using real TCGA data `data-raw/validate_tcga.R` runs the full reproduction. It lives in the [package repository](https://github.com/camilasauria/thyroidBRS), not in the installed package. The shape of it is below; the chunk is not evaluated because it downloads the whole TCGA-THCA expression set from GDC: ```{r tcga, eval = FALSE} library(TCGAbiolinks) # 1. A log2(TPM + 1) matrix of primary tumors, gene symbols x patient barcode. q <- GDCquery( project = "TCGA-THCA", data.category = "Transcriptome Profiling", data.type = "Gene Expression Quantification", workflow.type = "STAR - Counts" ) GDCdownload(q) se <- GDCprepare(q) # keep assay "tpm_unstranded", primary tumors (-01), collapse duplicated # gene symbols, then log2(x + 1) -> expr # 2. Reference labels: driver mutation, NOT the published class. ref <- read.csv(system.file("extdata", "thca_reference.csv", package = "thyroidBRS")) mut <- setNames(ref$driver_group, ref$patient) mut <- mut[mut %in% c("BRAF_V600E", "RAS")] # 3. Fit, then score everyone. fit <- brs_fit(expr, mut) preds <- predict(fit, expr) # 4. Validate against something that did not enter the fit. published <- setNames(ref$published_class, ref$patient) validate_brs(preds, published[!is.na(published)], fit = fit) ``` `thca_reference.csv` has one row per primary tumor of the study (496): `sample` (aliquot barcode), `patient`, `driver_group` (`BRAF_V600E`, `RAS` or `OTHER`), `has_exome`, `published_brs` and `published_class`. `OTHER` covers both tumors with another driver and the 94 that were never exome sequenced; `has_exome` tells them apart. For the fit this does not matter, since `OTHER` is simply "unlabeled", but it does if you use the table for anything else. Doing this today scores 505 primary tumors — 114 more than the published 391 — and agrees with the published score at Spearman 0.941, with 96.4% class concordance on the 111 tumors that were not part of the reference groups — tumors that are neither BRAF-V600E- nor RAS-mutant, and that were all exome sequenced. ## Session info ```{r session-info} sessionInfo() ``` ## Reference Agrawal N, Akbani R, Aksoy BA, et al. Integrated Genomic Characterization of Papillary Thyroid Carcinoma. *Cell*. 2014;159(3):676-690. [doi:10.1016/j.cell.2014.09.050](https://doi.org/10.1016/j.cell.2014.09.050)