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.2755549 0.2365942
## 2 p.(A107C) 107 A C missense 0.1134388 0.2255906
## 3 p.(A107D) 107 A D missense 0.3873864 0.2108217
## 4 p.(A107E) 107 A E missense -0.6415181 0.2534535
## 5 p.(A107F) 107 A F missense -0.1401362 0.2692807
## 6 p.(A107G) 107 A G missense -1.0571610 0.2417313
## lfsr lfsr.neg lfsr.pos test.neg test.pos label
## 1 1.220757e-01 1.220757e-01 0.87792430 FALSE FALSE Neutral
## 2 3.075340e-01 6.924660e-01 0.30753402 FALSE FALSE Neutral
## 3 3.306753e-02 9.669325e-01 0.03306753 FALSE FALSE Neutral
## 4 5.685139e-03 5.685139e-03 0.99431486 TRUE FALSE Neg
## 5 3.013891e-01 3.013891e-01 0.69861090 FALSE FALSE Neutral
## 6 6.119409e-06 6.119409e-06 0.99999388 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 67 359 102
## missense 2965 3682 3614
## synonymous 100 26 390In 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.(D5M) 5 D M missense 3.327057 0.3870681
## 2 p.(G508del) 508 G del deletion 3.204218 0.3741755
## 3 p.(D5W) 5 D W missense 2.724841 0.3856855
## 4 p.(A506L) 506 A L missense 2.636134 0.3641000
## 5 p.(N156del) 156 N del deletion 2.608736 0.2494461
## lfsr lfsr.neg lfsr.pos test.neg test.pos label
## 1 4.143899e-18 1 4.143899e-18 FALSE TRUE LOF
## 2 5.478980e-18 1 5.478980e-18 FALSE TRUE LOF
## 3 8.034936e-13 1 8.034936e-13 FALSE TRUE LOF
## 4 2.241158e-13 1 2.241158e-13 FALSE TRUE LOF
## 5 6.725221e-26 1 6.725221e-26 FALSE TRUE LOF## variants position wildtype mutation type mean sd
## 1 p.(A219V) 219 A V missense -3.617156 0.3442180
## 2 p.(L505C) 505 L C missense -3.040896 0.4894418
## 3 p.(H304M) 304 H M missense -2.904934 0.4207863
## 4 p.(C473S) 473 C S missense -2.777315 0.4695054
## 5 p.(A21E) 21 A E missense -2.719693 0.3956385
## lfsr lfsr.neg lfsr.pos test.neg test.pos label
## 1 3.954126e-26 3.954126e-26 1 TRUE FALSE GOF
## 2 2.599326e-10 2.599326e-10 1 TRUE FALSE GOF
## 3 2.535298e-12 2.535298e-12 1 TRUE FALSE GOF
## 4 1.655296e-09 1.655296e-09 1 TRUE FALSE GOF
## 5 3.117210e-12 3.117210e-12 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.
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 test.neg
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <lgl>
## 1 2 -0.736 0.261 1.18 0.484 0.266 0.734 0.266 FALSE
## 2 3 -0.773 0.205 0.615 0.256 0.104 0.896 0.104 FALSE
## 3 4 0.297 0.231 0.959 0.382 0.622 0.378 0.378 FALSE
## 4 5 1.41 0.195 0.702 0.274 0.978 0.0225 0.0225 FALSE
## 5 6 -0.602 0.208 0.737 0.314 0.207 0.793 0.207 FALSE
## 6 7 0.364 0.252 1.23 0.464 0.616 0.384 0.384 FALSE
## # ℹ 2 more variables: test.pos <lgl>, label <chr>