No result object created

[sofiedemeyer]

I tried running SCEPTRE on my dataset with the following code:

`cells = intersect(colnames(expression_matrix), colnames(perturbation_matrix)) expression_matrix = expression_matrix[,cells] perturbation_matrix = perturbation_matrix[,cells]

lg_mRNA_lib_size <- log(colSums(expression_matrix)) covariate_matrix <- data.frame(lg_mRNA_lib_size = lg_mRNA_lib_size)

gene_gRNA_pairs <- expand.grid(gene_id = row.names(expression_matrix), gRNA_id = row.names(perturbation_matrix)) storage_dir <- tempdir() side <- "left" pod_sizes <- c(gene = nrow(expression_matrix), gRNA = nrow(perturbation_matrix), pair = nrow(gene_gRNA_pairs)) result <- run_sceptre_in_memory(storage_dir, expression_matrix, perturbation_matrix, covariate_matrix, gene_gRNA_pairs, side, pod_sizes)

result$p_value_adj <- p.adjust(p = result$p_value, method = "BH")`

However, when running the command run_sceptre_in_memory it runs like everything is fine, but I do not get a result object returned. CSV-files with the expression and the perturbation matrix can be found here: https://we.tl/t-tc9qIpdnGx

What could be the problem?

[timothy-barry]

Sorry -- closed issue by accident. Reopening it.

[timothy-barry]

Thanks for opening the issue.

I think the issue likely is that you are passing a numeric perturbation matrix rather than a binary perturbation matrix to the run_sceptre_in_memory function. The function assumes that the perturbation matrix is binary, but to make the function easier to use, I updated it so that it automatically converts a numeric perturbation matrix into a binary perturbation matrix (via a simple thresholding operation) when the user passes a numeric perturbation matrix. I also updated run_sceptre_in_memory so that it sets the gene, gRNA and pair "pod sizes" to reasonable default values. These "pod sizes" have no biological or statistical significance and only affect the amount of parallelization that is used.

Please download the sceptre package again and try the following code. I stored your matrices "expression_matrix.csv" and "perturbation_matrix.csv" on my desktop, so you will have to update the file paths on lines 9 and 15 to point to the corresponding locations on your machine. Could you please try this and let me know if it works?

library(sceptre)
library(dplyr)
library(readr)

# set seed for reproducibility
set.seed(4)

# gene expression matrix
gene_exp_df <- read_csv("~/Desktop/expression_matrix.csv", col_names = TRUE) %>% rename("gene_id" = "...1")
gene_ids <- gene_exp_df$gene_id
gene_exp_matrix <- gene_exp_df %>% select(-gene_id) %>% as.matrix()
row.names(gene_exp_matrix) <- gene_ids

# perturbation matrix
pert_df <- read_csv("~/Desktop/perturbation_matrix.csv", col_names = TRUE) %>% rename("pert_id" = "...1")
pert_ids <- pert_df$pert_id
pert_matrix <- pert_df %>% select(-pert_id) %>% as.matrix()
row.names(pert_matrix) <- pert_ids

# Carry out the analysis
# 1. ensure gene expression matrix and perturbation matrix have the same cells in the same order
cell_barcodes <- intersect(colnames(pert_matrix), colnames(gene_exp_matrix))
gene_exp_matrix <- gene_exp_matrix[, cell_barcodes]
pert_matrix <- pert_matrix[, cell_barcodes]

# 2. compute the covariate matrix; include log gene and gRNA lib sizes.
covariate_matrix <- data.frame(lg_mRNA_lib_size = colSums(gene_exp_matrix) %>% log(),
                               lg_gRNA_lib_size = colSums(pert_matrix) %>% log)

# 3. create the data frame of gene-gRNA pairs; analyze only 50 pairs for simplicity
gRNA_gene_pairs <- expand.grid(gRNA_id =  row.names(pert_matrix),
                               gene_id = row.names(gene_exp_matrix)) %>% slice_sample(n = 50)

# 4. Run method
temp_dir <- tempdir()
result <- run_sceptre_in_memory(storage_dir = temp_dir,
                                expression_matrix = gene_exp_matrix,
                                perturbation_matrix = pert_matrix,
                                covariate_matrix = covariate_matrix,
                                gene_gRNA_pairs = gRNA_gene_pairs,
                                side = "left",
                                pod_sizes = c(gene = 5, gRNA = 5, pair = 5))

# Examine the result and correct the p-values
result$p_value_adj <- p.adjust(p = result$p_value, method = "BH")
discoveries <- result %>% filter(p_value_adj < 0.1) %>% select(gene_id, gRNA_id, p_value_adj)

[sofiedemeyer]

Thanks! I tried it, and is works fine now..

[timothy-barry]

Good to hear. Would it be OK for me to include the above script (that uses your data) as an example on the website? I'd just link to this script as an example.

[sofiedemeyer]

Yes, you can use the script with my data as an example. For the data loading part, I changed that to only two commands:

gene_exp_matrix = as.matrix(read.csv("expression_matrix.csv", row.names = 1, header = T)) pert_matrix = as.matrix(read.csv("perturbation_matrix.csv", header = T, row.names = 1))

Might be simpler..

From: Tim @.> Reply to: Katsevich-Lab/sceptre @.> Date: Friday, 7 January 2022 at 16:38 To: Katsevich-Lab/sceptre @.> Cc: sofiedemeyer @.>, Author @.> Subject: Re: [Katsevich-Lab/sceptre] No result object created (Issue #9) Resent from: @.>

Good to hear. Would it be OK for me to include the above script (that uses your data) as an example on the website? I'd just link to this script as an example.

— Reply to this email directly, view it on GitHubhttps://github.com/Katsevich-Lab/sceptre/issues/9#issuecomment-1007504031, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AM2U3GBRWGWKD7QCAUCT5PLUU4CGFANCNFSM5LMRD7RQ. Triage notifications on the go with GitHub Mobile for iOShttps://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Androidhttps://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub. You are receiving this because you authored the thread.Message ID: @.> [ { @.": "http://schema.org", @.": "EmailMessage", "potentialAction": { @.": "ViewAction", "target": "https://github.com/Katsevich-Lab/sceptre/issues/9#issuecomment-1007504031", "url": "https://github.com/Katsevich-Lab/sceptre/issues/9#issuecomment-1007504031", "name": "View Issue" }, "description": "View this Issue on GitHub", "publisher": { @.***": "Organization", "name": "GitHub", "url": "https://github.com" } } ]

[timothy-barry]

Ah, that definitely is simpler. Thanks and best of luck!