Interpret Rosace Results
To have a more comprehensive understanding and clearer view of the score distribution, we’ve made available a precomputed object based on the complete OCT1 dataset. This may help you draw parallel to your own data. libra
# obtain the data
data("oct1_rosace_scored")
# extract scores
scores.data <- OutputScore(oct1_rosace_scored, name = "1SM73_ROSACE", sig.test = 0.05)
head(scores.data)## variants position wildtype mutation type mean sd
## 1 p.(A107A) 107 A A synonymous -0.2766760 0.2337295
## 2 p.(A107C) 107 A C missense 0.1195266 0.2269136
## 3 p.(A107D) 107 A D missense 0.3834164 0.2090132
## 4 p.(A107E) 107 A E missense -0.6452607 0.2549238
## 5 p.(A107F) 107 A F missense -0.1401713 0.2644013
## 6 p.(A107G) 107 A G missense -1.0386441 0.2606083
## lfsr lfsr.neg lfsr.pos test.neg test.pos label
## 1 1.182571e-01 1.182571e-01 0.88174285 FALSE FALSE Neutral
## 2 2.991837e-01 7.008163e-01 0.29918374 FALSE FALSE Neutral
## 3 3.329643e-02 9.667036e-01 0.03329643 FALSE FALSE Neutral
## 4 5.683809e-03 5.683809e-03 0.99431619 TRUE FALSE Neg
## 5 2.980054e-01 2.980054e-01 0.70199463 FALSE FALSE Neutral
## 6 3.367477e-05 3.367477e-05 0.99996633 TRUE FALSE NegData visualization
First, we can visualize the distribution of scores across mutation
types with the scoreDensity function.
scoreDensity(scores.data,
hist = FALSE,
savedir = "../tests/results/stan/assayset_full/plot/",
name = "DensityPlot_1SM73")
For other visualizations including visualizing position-wise distributions, refer to the Visualizing functional score results vignette.
Hypothesis testing
A common goal in DMS experiments is to determine which mutations
cause a GOF or a LOF. One natural approach is to conduct hypothesis
testing to check whether a score is significantly different from 0.
lfsr stands for local false sign rate, and is defined as
the minimum of the density on the left of 0 (lfsr.pos) and
on the right of 0 (lfsr.neg). The OutputScore
function automatically does a two-sided test based on the
lfsr with threshold 0.05: if lfsr.neg <
0.05/2, we are confident that the mass of the distribution is on the
left side of 0 and label the variant “Neg”, and vice versa.
For the OCT1 data, positive scores represent loss-of-function and negative scores represent gain-of-function. You can update the label to “LOF”, “Neutral”, and “GOF” (instead of “Neg”, “Neutral”, and “Pos”) by the following function.
scores.data <- scores.data %>%
mutate(label = case_when(
label == "Neg" ~ "GOF",
label == "Pos" ~ "LOF",
label == "Neutral" ~ "Neutral"
))Keep in mind that setting a threshold is a somewhat counter-Bayesian approach. The threshold should be adaptable, depending on the nature of the experiment and distribution of the data. You may be more or less conservative with thresholding based on the context and prior knowledge.
##
## GOF LOF Neutral
## deletion 53 368 107
## missense 2841 3758 3662
## synonymous 103 26 387In the growth screen, hypothesis testing for negative scores is much harder than that for the positive scores because random dropout of a variant is common, even for the synonymous mutations, especially when the mass of LOF variant score is shifting towards positive (faster growth rate than synonymous). As a result, it is expected to see more false discovery in the “GOF” bin in the OCT1 data. One might adjust the threshold to be more stringent when discovering “Neg” variants, and looser when discovering “Pos” variants.
Rank the variants
Sometimes, we want to select the most “LOF” or “GOF” variants for
validation. Our most confident set of variants are those that have a
high absolute value of the score and small lfsr close to
0.
For example, the 5 most “LOF” and “GOF” variants in the OCT1 data can be found by the function below.
## variants position wildtype mutation type mean sd
## 1 p.(N156del) 156 N del deletion 2.624989 0.2537912
## 2 p.(L211R) 211 L R missense 2.586086 0.2678294
## 3 p.(Y127Q) 127 Y Q missense 2.581790 0.2009293
## 4 p.(V466D) 466 V D missense 2.485453 0.2608994
## 5 p.(L25K) 25 L K missense 2.441372 0.1932449
## lfsr lfsr.neg lfsr.pos test.neg test.pos label
## 1 2.248601e-25 1 2.248601e-25 FALSE TRUE LOF
## 2 2.324369e-22 1 2.324369e-22 FALSE TRUE LOF
## 3 4.342803e-38 1 4.342803e-38 FALSE TRUE LOF
## 4 8.135244e-22 1 8.135244e-22 FALSE TRUE LOF
## 5 6.894090e-37 1 6.894090e-37 FALSE TRUE LOF## variants position wildtype mutation type mean sd
## 1 p.(F460H) 460 F H missense -2.169316 0.2889278
## 2 p.(P79G) 79 P G missense -1.914219 0.2858500
## 3 p.(P120A) 120 P A missense -1.803614 0.2948415
## 4 p.(Y279G) 279 Y G missense -1.795650 0.2917418
## 5 p.(V103H) 103 V H missense -1.791982 0.2895312
## lfsr lfsr.neg lfsr.pos test.neg test.pos label
## 1 2.998224e-14 2.998224e-14 1 TRUE FALSE GOF
## 2 1.066723e-11 1.066723e-11 1 TRUE FALSE GOF
## 3 4.760717e-10 4.760717e-10 1 TRUE FALSE GOF
## 4 3.755610e-10 3.755610e-10 1 TRUE FALSE GOF
## 5 3.022497e-10 3.022497e-10 1 TRUE FALSE GOFPosition-level estiamte
When running with the position mode (default), Rosace infers position-specific mean estimate (“phi”) and position-specific variance estimate (“sigma2”). The hypothesis testing columns can be interpreted the same way as the ones in the variant-specific section.
If the user includes synonymous and/or missense mutations when running Rosace, it would exclude those variants from the position-level estimates.
scores.data.list <- OutputScore(oct1_rosace_scored, pos.info = TRUE, name = "1SM73_ROSACE", sig.test = 0.05)
scores.var <- scores.data.list$df_variant
scores.pos <- scores.data.list$df_position
head(scores.pos)## # A tibble: 6 × 11
## pos phi_mean phi_sd sigma2_mean sigma2_sd lfsr.neg lfsr.pos lfsr
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2 -0.273 0.173 0.488 0.196 0.0577 9.42e- 1 5.77e- 2
## 2 3 -0.209 0.153 0.364 0.147 0.0860 9.14e- 1 8.60e- 2
## 3 4 0.469 0.176 0.586 0.225 0.996 3.77e- 3 3.77e- 3
## 4 5 1.25 0.147 0.349 0.138 1 8.55e-18 8.55e-18
## 5 6 -0.145 0.155 0.397 0.154 0.174 8.26e- 1 1.74e- 1
## 6 7 0.469 0.218 0.911 0.341 0.985 1.55e- 2 1.55e- 2
## # ℹ 3 more variables: test.neg <lgl>, test.pos <lgl>, label <chr>