gimap performs analysis of dual-targeting CRISPR screening data, with the goal of aiding the identification of genetic interactions (e.g. cooperativity, synthetic lethality) in models of disease and other biological contexts. gimap analyzes functional genomic data generated by the pgPEN (paired guide RNAs for genetic interaction mapping) approach, quantifying growth effects of single and paired gene knockouts upon application of a CRISPR library. A multitude of CRISPR screen types can be used for this analysis, with helpful descriptions found in this review (https://www.nature.com/articles/s43586-021-00093-4). Use of pgPEN and GI-mapping in a paired gRNA format can be found here (https://pubmed.ncbi.nlm.nih.gov/34469736/).
library(gimap)
#> Warning: replacing previous import 'dplyr::group_rows' by
#> 'kableExtra::group_rows' when loading 'gimap'
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
output_dir <- "output_treatment"
dir.create(output_dir, showWarnings = FALSE)In this example we are going to examine a dataset with drug treatments and controls.
Let’s examine this example pgPEN counts table. It’s divided into columns containing:
id: an ID corresponding to the names of paired
guidesseq_1: gRNA sequence 1, targeting “paralog A”seq_2: gRNA sequence 2, targeting “paralog B”pretreatment: The pretreated controldmsoA: Vehicle only for Replicate AdmsoB: Vehicle only for Replicate Bdrug1A: Drug treatment for Replicate Adrug2B: Drug treatment for Replicate BFor the purposes of this tutorial, we can grab example data from the package.
example_data <- get_example_data("count_treatment")
#> Rows: 33170 Columns: 7
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: "\t"
#> chr (2): id, gene
#> dbl (5): pretreatment, dmsoA, dmsoB, drug1A, drug1B
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.The metadata you have may vary slightly from this but you’ll want to make sure you have the essential variables and information regarding how you collected your data.
colnames(example_data)
#> [1] "id" "gene" "pretreatment" "dmsoA" "dmsoB"
#> [6] "drug1A" "drug1B"We’re going to set up three datasets that we will provide to the
set_up() function to create a gimap dataset
object.
counts - the counts generated from pgPENpg_ids - the IDs that correspond to the rows of the
counts and specify the constructsample_metadata - metadata that describes the columns
of the counts including their timepointsFor this example we are using a treatment dataset where there is a pretreatment sample, two vehicle replicates and two drug treated replicates.
counts <- example_data %>%
select(c("pretreatment", "dmsoA", "dmsoB", "drug1A", "drug1B")) %>%
as.matrix()The next datasets are metadata that describe the dimensions of the count data.
pg_id are just the unique IDs listed in the same
order/sorted the same way as the count data and can be used for mapping
between the count data and the metadata. This is required and very
important because it is necessary to know the IDs and be able to map
them to pgRNA constructs and counts data.
Sample metadata is the information that describes the samples and is sorted the same order as the columns in the count data.
You need to have two columns in the metadata you provide. You’ll need
to specify the names of these columns in the
gimap_annotate() function.
col_names - Must match the colnames of the counts data
being submittedtreatments - A factor variable that describes the
treatments for these data. The label for the control specified in this
column will need to be supplied to the gimap_normalize()
function.
sample_metadata <- data.frame(
col_names = c("pretreatment", "dmsoA", "dmsoB", "drug1A", "drug1B"),
drug_treatment = as.factor(c("pretreatment", "dmso", "dmso", "drug", "drug"))
)We’ll need to provide counts, pg_ids and
sample_metadata to setup_data().
Now let’s setup our data using setup_data(). Optionally
we can provide the metadata in this function as well so that it is
stored with the data.
gimap_dataset <- setup_data(
counts = counts,
pg_ids = pg_ids,
sample_metadata = sample_metadata
)You’ll notice that this set up gives us a list of formatted data.
This contains the original counts we gave setup_data()
function but also normalized counts, and the total counts per
sample.
raw_counts: The original counts data that illustrates
the number of cells alive in the sample. This data has samples as the
columns and the paired guide constructs as rows.counts_per_sample: Add up all the counts for each
sample over all of the paired guide designs.count_norm - For each sample, the data is normalized
-log10(( counts +1) / total counts for the sample over all the pg designs ))
cpm - For each sample this is calculated by taking the
counts / total counts for the sample over all the pg designs)*1e6
log2cpm: log-2 transformed counts per million this is
calculated by log2(cpms + 1)
sample_metadata: Metadata that describes the samples.
This likely includes the time point information, replicates, sample IDs,
and any other additional information that is needed regarding the
experimental setup.
str(gimap_dataset)
#> List of 10
#> $ raw_counts : num [1:33170, 1:5] 240 244 24 175 50 128 24 65 311 66 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : NULL
#> .. ..$ : chr [1:5] "pretreatment" "dmsoA" "dmsoB" "drug1A" ...
#> $ counts_per_sample: Named num [1:5] 6214019 6662533 6640876 7115331 6894220
#> ..- attr(*, "names")= chr [1:5] "pretreatment" "dmsoA" "dmsoB" "drug1A" ...
#> $ transformed_data :List of 3
#> ..$ count_norm: num [1:33170, 1:5] 4.41 4.4 5.4 4.55 5.09 ...
#> .. ..- attr(*, "dimnames")=List of 2
#> .. .. ..$ : NULL
#> .. .. ..$ : chr [1:5] "pretreatment" "dmsoA" "dmsoB" "drug1A" ...
#> ..$ cpm : num [1:33170, 1:5] 38.62 39.27 3.86 28.16 8.05 ...
#> .. ..- attr(*, "dimnames")=List of 2
#> .. .. ..$ : NULL
#> .. .. ..$ : chr [1:5] "pretreatment" "dmsoA" "dmsoB" "drug1A" ...
#> ..$ log2_cpm : num [1:33170, 1:5] 5.31 5.33 2.28 4.87 3.18 ...
#> .. ..- attr(*, "dimnames")=List of 2
#> .. .. ..$ : NULL
#> .. .. ..$ : chr [1:5] "pretreatment" "dmsoA" "dmsoB" "drug1A" ...
#> $ metadata :List of 2
#> ..$ pg_ids : tibble [33,170 × 1] (S3: tbl_df/tbl/data.frame)
#> .. ..$ id: chr [1:33170] "AADAC_AADACL2_pg1" "AADAC_AADACL2_pg10" "AADAC_AADACL2_pg11" "AADAC_AADACL2_pg12" ...
#> ..$ sample_metadata:'data.frame': 5 obs. of 2 variables:
#> .. ..$ col_names : chr [1:5] "pretreatment" "dmsoA" "dmsoB" "drug1A" ...
#> .. ..$ drug_treatment: Factor w/ 3 levels "dmso","drug",..: 3 1 1 2 2
#> $ filtered_data :List of 5
#> ..$ filter_step_run : logi FALSE
#> ..$ metadata_pg_ids : NULL
#> ..$ transformed_log2_cpm : NULL
#> ..$ removed_pg_ids : NULL
#> ..$ all_reps_zerocount_ids: NULL
#> $ comparisons : NULL
#> $ annotation : NULL
#> $ normalized_log_fc: NULL
#> $ crispr_score : NULL
#> $ results : NULL
#> - attr(*, "class")= chr [1:2] "list" "gimap_dataset"Let’s see how many rows of data we have.
nrow(gimap_dataset$transformed_data$log2_cpm)
#> [1] 33170The first step is running some quality checks on our data. The
run_qc() function will create a report we can look at to
assess this.
The report includes several visualizations of raw/unfiltered data:
log2cpm values for each individual sample:
this helps users identify samples that do not have a normal distribution
of reads and inform the upcoming filtering steps.This report also includes several visualizations after filters are applied:
There is a filter that flags pgRNA constructs where any of the time points have a count of zero. - We include a bar plot that shows the number of pgRNA constructs which have counts of zero in either 0, 1, 2, or 3 replicates. - We include a table that specifies how many pgRNAs would be filtered out by applying this filter.
There is a filter that flags pgRNA constructs that have low log2 CPM counts for the day 0 or plasmid time point. - The histogram of the log2 CPM values of pgRNA constructs at the plasmid time point mentioned earlier does have a dashed line specifying the lower outlier (or a user defined cutoff) and pgRNA constructs with a plasmid log2 CPM lower than that value can be filtered out. - We include a table that specifies how many pgRNAs would be filtered out by applying this filter.
There is a filter that flags pgRNA constructs that have low log2 CPM counts for the day 0 or plasmid time point. - The histogram of the log2 CPM values of pgRNA constructs at the plasmid time point mentioned earlier does have a dashed line specifying the lower outlier (or a user defined cutoff) and pgRNA constructs with a plasmid log2 CPM lower than that value can be filtered out. - We include a table that specifies how many pgRNAs would be filtered out by applying this filter.
run_qc(gimap_dataset,
output_file = "example_qc_report.Rmd",
overwrite = TRUE,
quiet = TRUE)You can take a look at an example QC report here.
After considering the QC report and which filters are
appropriate/desired for your data, you can apply filters to the data
using the gimap_filter function.
gimap_filtered <- gimap_dataset %>%
gimap_filter()How many rows of data do we have now after filtering?
nrow(gimap_filtered$filtered_data$transformed_log2_cpm)
#> [1] 32402
str(gimap_dataset$filtered_data)
#> List of 5
#> $ filter_step_run : logi FALSE
#> $ metadata_pg_ids : NULL
#> $ transformed_log2_cpm : NULL
#> $ removed_pg_ids : NULL
#> $ all_reps_zerocount_ids: NULLLet’s take a look at how many rows of data we have left:
nrow(gimap_dataset$filtered_data$transformed_log2_cpm)
#> NULLAs you can see from the output above, there are fewer pgRNA constructs in the filtered dataset following completion of filtering.
The filtering step also stores two tables of information that you may want to use or report.
$filtered_data$removed_pg_ids is a table that has the
pgRNA construct IDs that are removed following completion of filtering
in the id column and the relevant filter(s) that led to
removal as a comma separated list in the relevantFilters
column$filtered_data$all_reps_zerocount_ids is a table that
lists the IDs of pgRNA constructs which had a count of 0 for all final
timepoint replicates. These pgRNA constructs are NOT necessarily
filtered outNow that you’ve performed QC and filtering, the rest of the pipeline can be run
cell_line your data uses so that the correct
corresponding DepMap data is used for annotation. This function is
gimap_annotate()
gimap_normalize()
function. timepoints needs to be specified pointing to the
correct column names from the sample_data passed to the
setup_data() function earlier.calc_crispr()
function.calc_gi() function.
gimap_dataset <- gimap_filtered %>%
gimap_annotate(cell_line = "PC9") %>%
# Whatever is specified for "control_name" is what will be used to normalize other data points
gimap_normalize(
treatments = "drug_treatment",
control_name = "pretreatment"
) %>%
calc_crispr() %>%
calc_gi()
#> Annotating Data
#> Rows: 1884 Columns: 3
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (3): gene, gene_symbol, entrez_id
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> Normalizing Log Fold Change
#>
#> The following 2 IDs were not found in the annotation data:
#> c("DHFR2_nt1_g1", "DHFR2_nt1_g2")
#>
#> The input data for the IDs which were not found in the annotation data has been filtered out and will not be included in the analysis output.
#>
#> Calculating CRISPR score
#>
#> Calculating Genetic Interaction scoresHere’s what’s included in the GI Scores table:
head(gimap_dataset$gi_scores)
#> # A tibble: 6 × 8
#> pgRNA_target_double rep double_target_gi_score single_target_gi_score_1
#> <chr> <chr> <dbl> <dbl>
#> 1 AADAC_AADACL2 dmsoA_dmso 1.24 0.428
#> 2 AADAC_AADACL2 dmsoA_dmso 1.24 0.428
#> 3 AADAC_AADACL2 dmsoA_dmso 1.24 -0.109
#> 4 AADAC_AADACL2 dmsoA_dmso 1.24 -0.109
#> 5 AADAC_AADACL2 dmsoB_dmso 0.638 0.112
#> 6 AADAC_AADACL2 dmsoB_dmso 0.638 0.112
#> # ℹ 4 more variables: single_target_gi_score_2 <dbl>,
#> # expected_crispr_single_1 <dbl>, expected_crispr_single_2 <dbl>,
#> # expected_crispr_double <dbl>Take a look at a preview of the CRISPR scores:
head(gimap_dataset$crispr_score)
#> # A tibble: 6 × 13
#> pg_ids rep double_crispr_score single_crispr_score_1 single_crispr_score_2
#> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 AADAC_A… dmso… 1.63 0.399 1.19
#> 2 AADAC_A… dmso… 1.63 0.399 0.0519
#> 3 AADAC_A… dmso… 1.63 0.00407 1.19
#> 4 AADAC_A… dmso… 1.63 0.00407 0.0519
#> 5 AADAC_A… dmso… 0.716 -0.223 1.15
#> 6 AADAC_A… dmso… 0.716 -0.223 -0.353
#> # ℹ 8 more variables: pgRNA_target_double <chr>, gene_symbol_1 <chr>,
#> # gene_symbol_2 <chr>, mean_single_target_crispr_1 <dbl>,
#> # mean_single_target_crispr_2 <dbl>, pgRNA1_seq <chr>, pgRNA2_seq <chr>,
#> # mean_double_control_crispr <dbl>Take a look at the overall stats:
gimap_dataset$results$overall
#> # A tibble: 8 × 6
#> # Groups: rep [4]
#> rep term estimate std.error statistic p.value
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 dmsoA_dmso (Intercept) 0.0257 0.00160 16.0 6.07e- 58
#> 2 dmsoA_dmso mean_expected_single_crispr 0.441 0.00228 194. 0
#> 3 dmsoB_dmso (Intercept) -0.0435 0.00176 -24.7 5.59e-135
#> 4 dmsoB_dmso mean_expected_single_crispr 0.451 0.00232 195. 0
#> 5 drug1A_drug (Intercept) 0.122 0.00300 40.7 0
#> 6 drug1A_drug mean_expected_single_crispr 0.424 0.00221 192. 0
#> 7 drug1B_drug (Intercept) -0.0165 0.00258 -6.40 1.56e- 10
#> 8 drug1B_drug mean_expected_single_crispr 0.429 0.00227 189. 0 Lastly, let’s take a look at the by target results. We can arrange by
results using dplyr::arrange()
gimap_dataset$results$by_target %>%
dplyr::arrange(fdr_vals_wil_dmsoA_dmso)
#> # A tibble: 1,030 × 17
#> targets p_val_ttest_dmsoA_dmso p_val_ttest_dmsoB_dmso p_val_ttest_drug1A_d…¹
#> <chr> <dbl> <dbl> <dbl>
#> 1 RBBP7_R… 0.00000547 0.00233 0.235
#> 2 AP2A1_A… 0.000174 0.000298 0.0515
#> 3 ATP6V1E… 0.0000251 0.0384 0.319
#> 4 KCTD2_K… 0.00000696 0.0121 0.0103
#> 5 GRIP2_G… 0.00000720 0.187 0.122
#> 6 MPPED2_… 0.0000250 0.328 0.426
#> 7 CSTF2T_… 0.000132 0.195 0.880
#> 8 LHFPL6_… 0.0000130 0.985 0.356
#> 9 HRASLS2… 0.0000154 0.966 0.120
#> 10 PSMB6_P… 0.0000310 0.621 0.0127
#> # ℹ 1,020 more rows
#> # ℹ abbreviated name: ¹p_val_ttest_drug1A_drug
#> # ℹ 13 more variables: p_val_ttest_drug1B_drug <dbl>,
#> # p_val_wil_dmsoA_dmso <dbl>, p_val_wil_dmsoB_dmso <dbl>,
#> # p_val_wil_drug1A_drug <dbl>, p_val_wil_drug1B_drug <dbl>,
#> # fdr_vals_ttest_dmsoA_dmso <dbl>, fdr_vals_ttest_dmsoB_dmso <dbl>,
#> # fdr_vals_ttest_drug1A_drug <dbl>, fdr_vals_ttest_drug1B_drug <dbl>, …We can save the genetic interactions scores like this:
Similarly let’s save the CRISPR scores.
We can extract and save the stats results to the TSV
readr::write_tsv(
gimap_dataset$results$overall,
file.path(output_dir, "overall_stats.tsv")
)
readr::write_tsv(
gimap_dataset$results$by_target,
file.path(output_dir, "by_target_stats.tsv")
)We can save all these data as an RDS.
saveRDS(gimap_dataset, "gimap_dataset_final_treatment.RDS")This is just for provenance purposes.
sessionInfo()
#> R version 4.4.0 (2024-04-24)
#> Platform: x86_64-apple-darwin20
#> Running under: macOS 15.0.1
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> time zone: America/New_York
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] dplyr_1.1.4 gimap_0.1.0
#>
#> loaded via a namespace (and not attached):
#> [1] tidyr_1.3.1 sass_0.4.9 utf8_1.2.4 generics_0.1.3
#> [5] xml2_1.3.6 stringi_1.8.4 hms_1.1.3 digest_0.6.37
#> [9] magrittr_2.0.3 RColorBrewer_1.1-3 evaluate_1.0.1 grid_4.4.0
#> [13] timechange_0.3.0 fastmap_1.2.0 jsonlite_1.8.9 backports_1.5.0
#> [17] purrr_1.0.2 fansi_1.0.6 viridisLite_0.4.2 scales_1.3.0
#> [21] textshaping_0.4.0 jquerylib_0.1.4 cli_3.6.3 crayon_1.5.3
#> [25] rlang_1.1.4 bit64_4.5.2 munsell_0.5.1 withr_3.0.1
#> [29] cachem_1.1.0 yaml_2.3.10 parallel_4.4.0 tools_4.4.0
#> [33] tzdb_0.4.0 colorspace_2.1-1 ggplot2_3.5.1 kableExtra_1.4.0
#> [37] broom_1.0.7 curl_5.2.3 vctrs_0.6.5 R6_2.5.1
#> [41] lifecycle_1.0.4 lubridate_1.9.3 snakecase_0.11.1 stringr_1.5.1
#> [45] bit_4.5.0 fs_1.6.4 htmlwidgets_1.6.4 vroom_1.6.5
#> [49] ragg_1.3.3 janitor_2.2.0 pkgconfig_2.0.3 desc_1.4.3
#> [53] pkgdown_2.1.1 pillar_1.9.0 bslib_0.8.0.9000 gtable_0.3.6
#> [57] glue_1.8.0 systemfonts_1.1.0 xfun_0.48 tibble_3.2.1
#> [61] tidyselect_1.2.1 rstudioapi_0.17.1 knitr_1.48 htmltools_0.5.8.1
#> [65] rmarkdown_2.28 svglite_2.1.3 readr_2.1.5 pheatmap_1.0.12
#> [69] compiler_4.4.0