A guide to the GEMINI R package

Introduction

GEMINI follows a basic workflow:

• Create Input (gemini_create_input)
• Initialize Model (gemini_initialize)
• Perform Coordinate-Ascent Variational Inference (CAVI) (gemini_inference)
• Score Interactions (gemini_score)

Using counts data derived from a combination CRISPR screen, and annotations for both samples/replicates and guide/gene IDs, GEMINI can identify genetic interactions such as synthetic lethality and recovery.

Model

GEMINI uses the following model:

\(D_{gihjl} = N(x_{gi} * y_{gl} + x_{hj}*x_{hl} + x_{gihj}*s_{ghl}, \tau_{gihjl})\)

• \(D\) is the observation, the log-fold change of a guide pair in a cell line.
• \(x\) is the reagent-level effect, reflecting the consistency of a single guide or guide pair across samples
• \(y\) is the gene-level effect, reflecting the effect of gene knockout in a single sample.
• \(s\) is the combination effect, reflecting the effect of the interaction between a gene pair

• \(\tau\) is the precision of the data

• \(g\) is a gene targeted by guide \(i\)
• \(h\) is a gene targeted by guide \(j\)

• \(l\) is the sample

Installation

The stable version used in the publication Zamanighomi et al., 2019 can be found on Bioconductor:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("gemini")

For the latest development version, you can download from GitHub using the devtools package:

if (!requireNamespace("devtools", quietly = TRUE))
    install.packages("devtools")

devtools::install_github("foo/bar", build_opts = c("--no-resave-data", "--no-manual"), build_vignettes = TRUE)

Import Big Papi data

This is the data published in Najm et al., 2018 (doi: 10.1038/nbt.4048). All data was obtained from Supplementary Tables 2 and 3.

library("gemini")
data("counts", "guide.annotation", "sample.replicate.annotation", package = "gemini")

Input

GEMINI takes a counts matrix as follows:

knitr::kable(head(counts[,1:5]), caption = "Counts matrix", align = 'l')

GEMINI also requires sample/replicate annotation and guide/gene annotation:

knitr::kable(head(sample.replicate.annotation), caption = "Sample/replicate annotations")

knitr::kable(head(guide.annotation[,1:3]), caption = "Guide/gene annotation")

These can be used to create a gemini.input object using the gemini_create_input function.

Input <- gemini_create_input(counts.matrix = counts,
                    sample.replicate.annotation = sample.replicate.annotation,
                    guide.annotation = guide.annotation,
                    ETP.column = 'pDNA', 
                    gene.column.names = c("U6.gene", "H1.gene"),
                    sample.column.name = "samplename",
                    verbose = TRUE)

## Merging sample annotations with colnames of counts.matrix...

## Merging guide annotations with rownames()...

## Created gemini input object.

# Note: ETP column can also be specified by column index 
# (e.g. ETP.column = c(1))

Pre-processing

GEMINI requires log-fold changes as an input, which are calculated using the gemini_calculate_lfc function. A pseudo-count (CONSTANT) of \(32\) is used by default.

Input %<>% gemini_calculate_lfc(normalize = TRUE, 
                                CONSTANT = 32)

Initialization and Inference

To initialize the CAVI approach, a gemini.model object is created using the gemini_initialize function. For large libraries, more cores can be specified to speed up the initialization (see ?gemini_parallelization).

To note, it is highly recommended that at least one negative control gene should be specified here, and all other genes should be paired with at least one negative control. We use CD81 in this example.

In the absence of a negative control gene, the median LFC of all guide pairs targeting each gene is used to initialize the CAVI approach. However, this is only reasonable in all-by-all format screens.

Also to note, the pattern_split argument must describe a separator used in the rownames of the counts.matrix. For example, a semi-colon (“;”) is used in the Big Papi data and therefore specified here. On the other hand, pattern_join is specified by the user to join the gene pairs. For example, if the U6.gene is BRCA2, and the H1.gene is PARP1, the output will be “BRCA2{pattern_join}PARP1”.

Model <- gemini_initialize(Input = Input, 
                  nc_gene = "CD81", 
                  pattern_join = ';',
                  pattern_split = ';', 
                  cores = 1,
                  verbose = TRUE)

## Initializing x

## Initializing y

## Initializing s

## Initializing tau

## Updating mae

## Model initialized.

Inference is performed with the gemini_inference function. For large libraries, more cores can be specified to speed up the inference (see ?gemini_parallelization).

Model %<>% gemini_inference(cores = 1,
                            verbose = FALSE)

Convergence

After running gemini_inference, the resulting convergence rate can be visualized. If the model is divergent, alternative parameters for the priors should be selected until convergence is achieved. This can be done through cross-validation. Details on this will be added soon.

gemini_plot_mae(Model)

Scoring and Visualization

To score genetic interactions, use the gemini_score function.

To note, at least one positive control gene (pc_gene) should be specified to remove interactions involving individually lethal genes. If no positive control is explicitly specified, lethality is estimated as described in Zamanighomi et al. In short, the most lethal 1st percentile of individual gene effects is treated as the positive control value.

Additionally, non-interacting gene pairs (nc_pairs) should be specified for the calculation of p-values and false discovery rates (FDRs). If not specified, only scores are calculated, which reflects the relative strengths of interactions in each sample.

# Use non-interacting gene pairs as nc_pairs. 
# A caveat here is that this set is constructed only using negative controls! 
# This probably leads to biased sampling of the null distribution, 
# thereby overestimating the number of significant hits, but still is useful in this case.
nc_pairs <- grep("6T|HPRT", rownames(Model$s), value = TRUE)

# An example of some nc_pairs...
head(nc_pairs, n = 5)

## [1] "6T;HPRT intron"       "HPRT intron;UBC"      "HPRT intron;HSP90AA1"
## [4] "CHEK2;HPRT intron"    "AKT1;HPRT intron"

Score <- gemini_score(Model = Model,
             pc_gene = "EEF2",
             nc_pairs = nc_pairs)

GEMINI’s scores for interactions can be seen in the Score object. Here, we show the top 10 interacting gene pairs in A549, with their scores across all cell lines.

knitr::kable(Score$strong[order(Score$strong[,"A549"], decreasing = TRUE)[1:10],], caption = "Strong scores for top 10 interactions from A549 in all samples")

Significant interactions can be identified through the FDR and p-value slots in the Score object. Again, we show the top 10 interacting gene pairs in A549, but now present FDRs across all cell lines.

knitr::kable(Score$fdr_strong[order(Score$fdr_strong[,"A549"], decreasing = FALSE)[1:10],], caption = "FDRs for top 10 interactions in A549")

To visualize these interactions, we can use the gemini_boxplot function. For example, in BRCA2-PARP1:

gemini_boxplot(Model = Model, 
               g = "BRCA2",
               h = "PARP1",
               nc_gene = "CD81",
               sample = "A549",
               show_inference = TRUE, 
               identify_guides = TRUE
               )

We can see that GEMINI makes adjustments to the individual gene effects of both BRCA2 and PARP1, and adjusts the sample-independent values of each to account for variation in the screen.

This boxplot can also be re-colored by the adjustment made to each guide or guide pair through the sample-independent effects. Guides/guide pairs with the least adjustment to x (grey) are considered to have the least variation within the screen.

gemini_boxplot(Model = Model, 
               g = "BRCA2",
               h = "PARP1",
               nc_gene = "CD81",
               sample = "A549",
               show_inference = TRUE, 
               color_x = TRUE
               )

Summary

GEMINI can be run on any counts matrix from a pairwise screen. GEMINI computes log fold changes and infers sample-dependent and sample-independent effects using CAVI. Then, GEMINI calculates interaction strength and significance.

The manuscript and more visualization tools will be made available soon.

Session Info

sessionInfo()

## R version 3.6.1 (2019-07-05)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows Server 2012 R2 x64 (build 9600)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=C                          
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] gemini_1.0.0
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.2        highr_0.8         pillar_1.4.2     
##  [4] compiler_3.6.1    mixtools_1.1.0    tools_3.6.1      
##  [7] digest_0.6.22     evaluate_0.14     tibble_2.1.3     
## [10] gtable_0.3.0      lattice_0.20-38   pkgconfig_2.0.3  
## [13] rlang_0.4.1       Matrix_1.2-17     parallel_3.6.1   
## [16] yaml_2.2.0        xfun_0.10         withr_2.1.2      
## [19] dplyr_0.8.3       stringr_1.4.0     knitr_1.25       
## [22] segmented_1.0-0   grid_3.6.1        tidyselect_0.2.5 
## [25] glue_1.3.1        R6_2.4.0          pbmcapply_1.5.0  
## [28] survival_2.44-1.1 rmarkdown_1.16    ggplot2_3.2.1    
## [31] purrr_0.3.3       magrittr_1.5      scales_1.0.0     
## [34] htmltools_0.4.0   MASS_7.3-51.4     splines_3.6.1    
## [37] assertthat_0.2.1  colorspace_1.4-1  labeling_0.3     
## [40] stringi_1.4.3     lazyeval_0.2.2    munsell_0.5.0    
## [43] crayon_1.3.4