knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
For more background on gimap and the calculations done here, read here
Besides installing the gimap package, you will also need to install wget if you do not already have it installed. This will allow you to download the annotation files needed to run gimap
.
To install the gimap
package you will need to run:
install.packages("gimap")
Or you can install the development version from GitHub:
install.packages("remotes") remotes::install_github("FredHutch/gimap")
library(gimap)
library(dplyr)
First we can create a folder we will save files to.
output_dir <- "output_timepoints" dir.create(output_dir, showWarnings = FALSE)
example_data <- get_example_data("count")
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 timepointscounts <- example_data %>% select(c("Day00_RepA", "Day22_RepA", "Day22_RepB", "Day22_RepC")) %>% as.matrix()
pg_id
are just the unique IDs listed in the same order/sorted the same way as the count data.
pg_ids <- example_data %>% dplyr::select("id")
Sample metadata is the information that describes the samples and is sorted the same order as the columns in the count data.
sample_metadata <- data.frame( col_names = c("Day00_RepA", "Day22_RepA", "Day22_RepB", "Day22_RepC"), day = as.numeric(c("0", "22", "22", "22")), rep = as.factor(c("RepA", "RepA", "RepB", "RepC")) )
We'll need to provide example_counts
, pg_ids
and sample_metadata
to setup_data()
.
gimap_dataset <- setup_data( counts = counts, pg_ids = pg_ids, sample_metadata = sample_metadata )
It's ideal to run quality checks first. The run_qc()
function will create a report we can look at to assess this.
run_qc(gimap_dataset, output_file = file.path(output_dir, "example_qc_report.Rmd"), overwrite = TRUE, plots_dir = "plots", quiet = TRUE )
You can take a look at an example QC report here.
gimap_dataset <- gimap_dataset %>% gimap_filter() %>% gimap_annotate(cell_line = "HELA") %>% gimap_normalize( timepoints = "day" ) %>% calc_gi()
Genetic interaction is calculated by:
pgRNA_target
- what gene(s) were targeted by this the original pgRNAs for these datamean_expected_cs
- the average expected genetic interaction scoremean_gi_score
- the average observer genetic interaction scoretarget_type
- describes whether the CRISPR design is targeting two genes ("gene_gene"), or a gene and a non targeting control ("gene_ctrl") or a targeting control and a gene ("ctrl_gene").p_val
- p values from the testing whether a double knockout construct is significantly different in its genetic interaction score from single targets. fdr
- False discovery rate corrected p valuesgimap_dataset$gi_scores %>% dplyr::arrange(fdr) %>% head() %>% knitr::kable(format = "html")
You can remove any samples from these plots by altering the reps_to_drop
argument.
plot_exp_v_obs_scatter(gimap_dataset) # Save it to a file ggsave(file.path(output_dir, "exp_v_obs_scatter.png"))
plot_rank_scatter(gimap_dataset) # Save it to a file ggsave(file.path(output_dir, "plot_rank_scatter.png"))
plot_volcano(gimap_dataset) # Save it to a file ggsave(file.path(output_dir, "volcano_plot.png"))
We can pick out a specific pair to plot.
# "CNOT8_CNOT7" is top result so let's plot that plot_targets(gimap_dataset, target1 = "CNOT8", target2 = "CNOT7") # Save it to a file ggsave(file.path(output_dir, "CNOT8_CNOT7.png"))
We can save all these data as an RDS or the genetic interaction scores themselves to a tsv file.
saveRDS(gimap_dataset, "gimap_dataset_final.RDS")
readr::write_tsv(gimap_dataset$gi_scores, "gi_scores.tsv")
This is just for provenance purposes.
sessionInfo()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.