Papillary thyroid carcinoma (PTC) is predominantly driven by either a
BRAFV600E mutation or a RAS-family mutation.
Agrawal et al. (2014) (Cell, doi: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:
brs_fit() — given a cohort with known
driver mutations for at least some samples, computes the two reference
centroids.predict() — scores any sample
against those centroids, including samples that were never labeled.thyroidBRS is installed from Bioconductor with
BiocManager:
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("thyroidBRS")Then load it:
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 BRAFV600E-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.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.
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.
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.
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))This is the point of the package: extending classification beyond the exome-sequenced subset.
preds_new <- predict(fit, expr[, unlabeled])
head(preds_new)
#> sample brs_score brs_scaled brs_class
#> 1 sample_37 -0.4191807 -0.9281559 Braf-like
#> 2 sample_38 0.2960093 0.5685397 Ras-like
#> 3 sample_39 -0.4516275 -1.0000000 Braf-like
#> 4 sample_40 0.3211279 0.6167844 Ras-like
#> 5 sample_41 -0.3681664 -0.8151993 Braf-like
#> 6 sample_42 0.3671792 0.7052344 Ras-like
table(preds_new$brs_class)
#>
#> Braf-like Ras-like
#> 12 12preds_all <- predict(fit, expr)
validate_brs(preds_all, labels, fit = fit)
#> Warning: Every compared sample also defined the centroids: this is a
#> resubstitution estimate, not validation. See ?validate_brs.
#> $n
#> [1] 36
#>
#> $n_concordant
#> [1] 36
#>
#> $pct_concordant
#> [1] 100
#>
#> $n_resubstituted
#> [1] 36
#>
#> $table
#> predicted
#> published Braf-like Ras-like
#> Braf-like 24 0
#> Ras-like 0 12Note 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:
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)
#> [1] 1brs_score() wraps brs_fit() +
predict() for the common case of scoring the same cohort
you fit on:
all_preds <- brs_score(expr, labels)
head(all_preds)
#> sample brs_score brs_scaled brs_class
#> 1 sample_1 -0.5023254 -0.9399192 Braf-like
#> 2 sample_2 -0.4793919 -0.8970075 Braf-like
#> 3 sample_3 -0.5214914 -0.9757814 Braf-like
#> 4 sample_4 -0.3989581 -0.7465049 Braf-like
#> 5 sample_5 -0.4654293 -0.8708817 Braf-like
#> 6 sample_6 -0.4003326 -0.7490769 Braf-likebrs_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:
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)
#> [1] TRUE
head(predict(fit_se, se, assay = "logexpr"))
#> sample brs_score brs_scaled brs_class
#> 1 sample_1 -0.5023254 -0.9399192 Braf-like
#> 2 sample_2 -0.4793919 -0.8970075 Braf-like
#> 3 sample_3 -0.5214914 -0.9757814 Braf-like
#> 4 sample_4 -0.3989581 -0.7465049 Braf-like
#> 5 sample_5 -0.4654293 -0.8708817 Braf-like
#> 6 sample_6 -0.4003326 -0.7490769 Braf-likeOn 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().
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:
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.
data-raw/validate_tcga.R runs the full reproduction. It
lives in the package
repository, 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:
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.
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] SummarizedExperiment_1.43.0 Biobase_2.73.2
#> [3] GenomicRanges_1.65.4 Seqinfo_1.3.2
#> [5] IRanges_2.47.5 S4Vectors_0.51.10
#> [7] BiocGenerics_0.59.12 generics_0.1.4
#> [9] MatrixGenerics_1.25.0 matrixStats_1.5.0
#> [11] thyroidBRS_0.99.7 BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] Matrix_1.7-6 jsonlite_2.0.0 compiler_4.6.1
#> [4] BiocManager_1.30.27 jquerylib_0.1.4 yaml_2.3.12
#> [7] fastmap_1.2.0 lattice_0.23-1 XVector_0.53.0
#> [10] R6_2.6.1 S4Arrays_1.13.2 knitr_1.52
#> [13] DelayedArray_0.39.8 maketools_1.3.2 bslib_0.12.0
#> [16] rlang_1.3.0 cachem_1.1.0 xfun_0.61
#> [19] sass_0.4.10 sys_3.4.3 otel_0.2.0
#> [22] SparseArray_1.13.4 cli_3.6.6 digest_0.6.39
#> [25] grid_4.6.1 lifecycle_1.0.5 evaluate_1.0.5
#> [28] buildtools_1.0.0 abind_1.4-8 rmarkdown_2.32
#> [31] tools_4.6.1 htmltools_0.5.9Agrawal 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