At this point, I think it is fair to say that the two cell types are
sufficiently different that they do not really belong together in a
single analysis.
Our main question of
interest
The data structure hs_macr contains our primary macrophages, which
are, as shown above, the data we can really sink our teeth into.
Note, we expect some errors when running the combine_de_tables()
because not all methods I use are comfortable using the ratio or ratios
contrasts we added in the ‘extras’ argument. As a result, when we
combine them into the larger output tables, those peculiar contrasts
fail. This does not stop it from writing the rest of the results,
however.
## test = deseq_pairwise(normalize_expt(hs_macr, filter=TRUE), model_batch = "svaseq", filter = TRUE, extra_contrasts = tmrc2_human_extra)
hs_macr_de <- all_pairwise(
hs_macr, model_batch = "svaseq", parallel = FALSE,
filter = TRUE,
extra_contrasts = tmrc2_human_extra)
##
## inf_z22 inf_z23 infsb_z22 infsb_z23 uninf_none uninfsb_none
## 11 12 12 11 4 4
## Removing 0 low-count genes (11756 remaining).
## Setting 2374 low elements to zero.
## transform_counts: Found 2374 values equal to 0, adding 1 to the matrix.
## Starting basic pairwise comparison.
## Basic step 0/3: Normalizing data.
## Basic step 0/3: Converting data.
## Basic step 0/3: Transforming data.
## Basic step 1/3: Creating mean and variance tables.
## Basic step 2/3: Performing 21 comparisons.
## Basic step 3/3: Creating faux DE Tables.
## Basic: Returning tables.
## Starting DESeq2 pairwise comparisons.
## The data should be suitable for EdgeR/DESeq/EBSeq.
## If they freak out, check the state of the count table
## and ensure that it is in integer counts.
## DESeq2 step 1/5: Including a matrix of batch estimates in the deseq model.
## converting counts to integer mode
## DESeq2 step 2/5: Estimate size factors.
## DESeq2 step 3/5: Estimate dispersions.
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## Using a parametric fitting seems to have worked.
## DESeq2 step 4/5: nbinomWaldTest.
## The contrast z23drugnodrug is not in the results.
## If this is not an extra contrast, then this is an error.
## The contrast z23z22drug is not in the results.
## If this is not an extra contrast, then this is an error.
## Starting edgeR pairwise comparisons.
## The data should be suitable for EdgeR/DESeq/EBSeq.
## If they freak out, check the state of the count table
## and ensure that it is in integer counts.
## EdgeR step 1/9: Importing and normalizing data.
## EdgeR step 2/9: Estimating the common dispersion.
## EdgeR step 3/9: Estimating dispersion across genes.
## EdgeR step 4/9: Estimating GLM Common dispersion.
## EdgeR step 5/9: Estimating GLM Trended dispersion.
## EdgeR step 6/9: Estimating GLM Tagged dispersion.
## EdgeR step 7/9: Running glmFit, switch to glmQLFit by changing the argument 'edger_test'.
## EdgeR step 8/9: Making pairwise contrasts.
## Length Class Mode
## title 1 -none- character
## notes 1 -none- character
## initial_metadata 71 data.frame list
## expressionset 1 ExpressionSet S4
## design 71 data.frame list
## conditions 54 -none- character
## batches 54 -none- character
## samplenames 54 -none- character
## colors 54 -none- character
## state 5 -none- list
## libsize 54 -none- numeric
## original_expressionset 1 ExpressionSet S4
## normalized 6 -none- list
## best_libsize 54 -none- numeric
## norm_result 6 -none- list
## Starting limma pairwise comparison.
## libsize was not specified, this parameter has profound effects on limma's result.
## Using the libsize from expt$best_libsize.
## Limma step 1/6: choosing model.
## Limma step 2/6: running limma::voom(), switch with the argument 'which_voom'.
## Using normalize.method = quantile for voom.
## Limma step 3/6: running lmFit with method: ls.
## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)
## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.
## Limma step 6/6: Writing limma outputs.
## Limma step 6/6: 1/17: Creating table: infsbz23_vs_infsbz22. Adjust = BH
## Limma step 6/6: 2/17: Creating table: infz22_vs_infsbz22. Adjust = BH
## Limma step 6/6: 3/17: Creating table: infz23_vs_infsbz22. Adjust = BH
## Limma step 6/6: 4/17: Creating table: uninfnone_vs_infsbz22. Adjust = BH
## Limma step 6/6: 5/17: Creating table: uninfsbnone_vs_infsbz22. Adjust = BH
## Limma step 6/6: 6/17: Creating table: infz22_vs_infsbz23. Adjust = BH
## Limma step 6/6: 7/17: Creating table: infz23_vs_infsbz23. Adjust = BH
## Limma step 6/6: 8/17: Creating table: uninfnone_vs_infsbz23. Adjust = BH
## Limma step 6/6: 9/17: Creating table: uninfsbnone_vs_infsbz23. Adjust = BH
## Limma step 6/6: 10/17: Creating table: infz23_vs_infz22. Adjust = BH
## Limma step 6/6: 11/17: Creating table: uninfnone_vs_infz22. Adjust = BH
## Limma step 6/6: 12/17: Creating table: uninfsbnone_vs_infz22. Adjust = BH
## Limma step 6/6: 13/17: Creating table: uninfnone_vs_infz23. Adjust = BH
## Limma step 6/6: 14/17: Creating table: uninfsbnone_vs_infz23. Adjust = BH
## Limma step 6/6: 15/17: Creating table: uninfsbnone_vs_uninfnone. Adjust = BH
## Limma step 6/6: 16/17: Creating table: z23drugnodrug_vs_z22drugnodrug. Adjust = BH
## Limma step 6/6: 17/17: Creating table: z23z22drug_vs_z23z22nodrug. Adjust = BH
## Limma step 6/6: 1/6: Creating table: infsbz22. Adjust = BH
## Limma step 6/6: 2/6: Creating table: infsbz23. Adjust = BH
## Limma step 6/6: 3/6: Creating table: infz22. Adjust = BH
## Limma step 6/6: 4/6: Creating table: infz23. Adjust = BH
## Limma step 6/6: 5/6: Creating table: uninfnone. Adjust = BH
## Limma step 6/6: 6/6: Creating table: uninfsbnone. Adjust = BH
## Length Class Mode
## title 1 -none- character
## notes 1 -none- character
## initial_metadata 71 data.frame list
## expressionset 1 ExpressionSet S4
## design 71 data.frame list
## conditions 54 -none- character
## batches 54 -none- character
## samplenames 54 -none- character
## colors 54 -none- character
## state 5 -none- list
## libsize 54 -none- numeric
## original_expressionset 1 ExpressionSet S4
## normalized 6 -none- list
## best_libsize 54 -none- numeric
## norm_result 6 -none- list
tmp_keepers <- tmrc2_human_keepers[13]
hs_macr_table <- combine_de_tables(
hs_macr_de,
keepers = tmrc2_human_keepers,
excel = glue("analyses/macrophage_de/hs_macr_drug_zymo_table_testing_macr_only-v{ver}.xlsx"))
## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.
## Did not find z22drugnodrug or z23drugnodrug.
## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.
## Did not find z22drugnodrug or z23drugnodrug.
## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.
## Did not find z23z22nodrug or z23z22drug.
## Warning in combine_extracted_plots(entry_name, combined, wanted_denominator, : I think this
## is an extra contrast table, the plots may be weird.
## Did not find z23z22nodrug or z23z22drug.
## Adding venn plots for z23nosb_vs_uninf.
## Adding venn plots for z22nosb_vs_uninf.
## Adding venn plots for z23nosb_vs_z22nosb.
## Adding venn plots for z23sb_vs_z22sb.
## Adding venn plots for z23sb_vs_z23nosb.
## Adding venn plots for z22sb_vs_z22nosb.
## Adding venn plots for z23sb_vs_sb.
## Adding venn plots for z22sb_vs_sb.
## Adding venn plots for z23sb_vs_uninf.
## Adding venn plots for z22sb_vs_uninf.
## Adding venn plots for sb_vs_uninf.
## Adding venn plots for extra_z2322.
## Adding venn plots for extra_drugnodrug.
combined_to_tsv(hs_macr_table, "macrophage")
hs_macr_sig <- extract_significant_genes(
hs_macr_table,
excel = glue("analyses/macrophage_de/hs_macr_drug_zymo_sig-v{ver}.xlsx"))
## There is no deseq_logfc column in the table.
## The columns are: ensembltranscriptid, ensemblgeneid, version, transcriptversion, hgncsymbol, description, genebiotype, cdslength, chromosomename, strand, startposition, endposition, transcript, edger_logfc, edger_adjp, limma_logfc, limma_adjp, edger_logcpm, edger_lr, edger_p, limma_ave, limma_t, limma_b, limma_p, limma_adjp_ihw, edger_adjp_ihw
## There is no deseq_logfc column in the table.
## The columns are: ensembltranscriptid, ensemblgeneid, version, transcriptversion, hgncsymbol, description, genebiotype, cdslength, chromosomename, strand, startposition, endposition, transcript, edger_logfc, edger_adjp, limma_logfc, limma_adjp, edger_logcpm, edger_lr, edger_p, limma_ave, limma_t, limma_b, limma_p, limma_adjp_ihw, edger_adjp_ihw
## There is no basic_logfc column in the table.
## The columns are: ensembltranscriptid, ensemblgeneid, version, transcriptversion, hgncsymbol, description, genebiotype, cdslength, chromosomename, strand, startposition, endposition, transcript, edger_logfc, edger_adjp, limma_logfc, limma_adjp, edger_logcpm, edger_lr, edger_p, limma_ave, limma_t, limma_b, limma_p, limma_adjp_ihw, edger_adjp_ihw
## There is no basic_logfc column in the table.
## The columns are: ensembltranscriptid, ensemblgeneid, version, transcriptversion, hgncsymbol, description, genebiotype, cdslength, chromosomename, strand, startposition, endposition, transcript, edger_logfc, edger_adjp, limma_logfc, limma_adjp, edger_logcpm, edger_lr, edger_p, limma_ave, limma_t, limma_b, limma_p, limma_adjp_ihw, edger_adjp_ihw
hs_macr_highsig <- extract_significant_genes(
hs_macr_table, min_mean_exprs = high_expression, exprs_column = high_expression_column,
excel = glue("analyses/macrophage_de/hs_macr_drug_zymo_highsig-v{ver}.xlsx"))
## Warning in get_sig_genes(this_table, lfc = lfc, p = p, z = z, n = n, column = this_fc_column,
## : The column deseq_basemean does not appears to be in the table, cannot filter by expression.
## Warning in get_sig_genes(this_table, lfc = lfc, p = p, z = z, n = n, column = this_fc_column,
## : The column deseq_basemean does not appears to be in the table, cannot filter by expression.
## Warning in get_sig_genes(this_table, lfc = lfc, p = p, z = z, n = n, column = this_fc_column,
## : The column deseq_basemean does not appears to be in the table, cannot filter by expression.
## Warning in get_sig_genes(this_table, lfc = lfc, p = p, z = z, n = n, column = this_fc_column,
## : The column deseq_basemean does not appears to be in the table, cannot filter by expression.
## There is no deseq_logfc column in the table.
## The columns are: ensembltranscriptid, ensemblgeneid, version, transcriptversion, hgncsymbol, description, genebiotype, cdslength, chromosomename, strand, startposition, endposition, transcript, edger_logfc, edger_adjp, limma_logfc, limma_adjp, edger_logcpm, edger_lr, edger_p, limma_ave, limma_t, limma_b, limma_p, limma_adjp_ihw, edger_adjp_ihw
## There is no deseq_logfc column in the table.
## The columns are: ensembltranscriptid, ensemblgeneid, version, transcriptversion, hgncsymbol, description, genebiotype, cdslength, chromosomename, strand, startposition, endposition, transcript, edger_logfc, edger_adjp, limma_logfc, limma_adjp, edger_logcpm, edger_lr, edger_p, limma_ave, limma_t, limma_b, limma_p, limma_adjp_ihw, edger_adjp_ihw
## There is no basic_logfc column in the table.
## The columns are: ensembltranscriptid, ensemblgeneid, version, transcriptversion, hgncsymbol, description, genebiotype, cdslength, chromosomename, strand, startposition, endposition, transcript, edger_logfc, edger_adjp, limma_logfc, limma_adjp, edger_logcpm, edger_lr, edger_p, limma_ave, limma_t, limma_b, limma_p, limma_adjp_ihw, edger_adjp_ihw
## There is no basic_logfc column in the table.
## The columns are: ensembltranscriptid, ensemblgeneid, version, transcriptversion, hgncsymbol, description, genebiotype, cdslength, chromosomename, strand, startposition, endposition, transcript, edger_logfc, edger_adjp, limma_logfc, limma_adjp, edger_logcpm, edger_lr, edger_p, limma_ave, limma_t, limma_b, limma_p, limma_adjp_ihw, edger_adjp_ihw
Our main questions
in U937
Let us do the same comparisons in the U937 samples, though I will not
do the extra contrasts, primarily because I think the dataset is less
likely to support them.
u937_de <- all_pairwise(u937_expt, model_batch = "svaseq", filter = TRUE)
##
## inf_z22 inf_z23 infsb_z22 infsb_z23 uninf_none uninfsb_none
## 3 3 3 3 1 1
## Removing 0 low-count genes (10751 remaining).
## Setting 5 low elements to zero.
## transform_counts: Found 5 values equal to 0, adding 1 to the matrix.
u937_table <- combine_de_tables(
u937_de,
keepers = u937_keepers,
excel = glue("analyses/macrophage_de/u937_drug_zymo_table-v{ver}.xlsx"))
## Error in names(x) <- value: 'names' attribute [8] must be the same length as the vector [6]
combined_to_tsv(u937_table, celltype = "u937")
## Error in combined_to_tsv(u937_table, celltype = "u937"): object 'u937_table' not found
u937_sig <- extract_significant_genes(
u937_table,
excel = glue("analyses/macrophage_de/u937_drug_zymo_sig-v{ver}.xlsx"))
## Error in extract_significant_genes(u937_table, excel = glue("analyses/macrophage_de/u937_drug_zymo_sig-v{ver}.xlsx")): object 'u937_table' not found
u937_highsig <- extract_significant_genes(
u937_table, min_mean_exprs = high_expression, exprs_column = high_expression_column,
excel = glue("analyses/macrophage_de/u937_drug_zymo_highsig-v{ver}.xlsx"))
## Error in extract_significant_genes(u937_table, min_mean_exprs = high_expression, : object 'u937_table' not found
Compare (no)Sb
z2.3/z2.2 treatments among macrophages
upset_plots_hs_macr <- upsetr_sig(
hs_macr_sig, both = TRUE,
contrasts = c("z23sb_vs_z22sb", "z23nosb_vs_z22nosb"))
upset_plots_hs_macr[["both"]]
groups <- upset_plots_hs_macr[["both_groups"]]
shared_genes <- attr(groups, "elements")[groups[[2]]] %>%
gsub(pattern = "^gene:", replacement = "")
length(shared_genes)
## [1] 387
shared_gp <- simple_gprofiler(shared_genes)
## No results to show
## Please make sure that the organism is correct or set significant = FALSE
## No results to show
## Please make sure that the organism is correct or set significant = FALSE
shared_gp[["pvalue_plots"]][["MF"]]
shared_gp[["pvalue_plots"]][["BP"]]
shared_gp[["pvalue_plots"]][["REAC"]]
drug_genes <- attr(groups, "elements")[groups[["z23sb_vs_z22sb"]]] %>%
gsub(pattern = "^gene:", replacement = "")
drugonly_gp <- simple_gprofiler(drug_genes)
## No results to show
## Please make sure that the organism is correct or set significant = FALSE
## No results to show
## Please make sure that the organism is correct or set significant = FALSE
## No results to show
## Please make sure that the organism is correct or set significant = FALSE
## No results to show
## Please make sure that the organism is correct or set significant = FALSE
## No results to show
## Please make sure that the organism is correct or set significant = FALSE
## No results to show
## Please make sure that the organism is correct or set significant = FALSE
drugonly_gp[["pvalue_plots"]][["BP"]]
I want to try something, directly include the u937 data in this…
both_sig <- hs_macr_sig
names(both_sig[["deseq"]][["ups"]]) <- paste0("macr_", names(both_sig[["deseq"]][["ups"]]))
names(both_sig[["deseq"]][["downs"]]) <- paste0("macr_", names(both_sig[["deseq"]][["downs"]]))
u937_deseq <- u937_sig[["deseq"]]
## Error in eval(expr, envir, enclos): object 'u937_sig' not found
names(u937_deseq[["ups"]]) <- paste0("u937_", names(u937_deseq[["ups"]]))
## Error in paste0("u937_", names(u937_deseq[["ups"]])): object 'u937_deseq' not found
names(u937_deseq[["downs"]]) <- paste0("u937_", names(u937_deseq[["downs"]]))
## Error in paste0("u937_", names(u937_deseq[["downs"]])): object 'u937_deseq' not found
both_sig[["deseq"]][["ups"]] <- c(both_sig[["deseq"]][["ups"]], u937_deseq[["ups"]])
## Error in eval(expr, envir, enclos): object 'u937_deseq' not found
both_sig[["deseq"]][["downs"]] <- c(both_sig[["deseq"]][["ups"]], u937_deseq[["downs"]])
## Error in eval(expr, envir, enclos): object 'u937_deseq' not found
summary(both_sig[["deseq"]][["ups"]])
## Length Class Mode
## macr_z23nosb_vs_uninf 49 data.frame list
## macr_z22nosb_vs_uninf 49 data.frame list
## macr_z23nosb_vs_z22nosb 49 data.frame list
## macr_z23sb_vs_z22sb 49 data.frame list
## macr_z23sb_vs_z23nosb 49 data.frame list
## macr_z22sb_vs_z22nosb 49 data.frame list
## macr_z23sb_vs_sb 49 data.frame list
## macr_z22sb_vs_sb 49 data.frame list
## macr_z23sb_vs_uninf 49 data.frame list
## macr_z22sb_vs_uninf 49 data.frame list
## macr_sb_vs_uninf 49 data.frame list
## macr_extra_z2322 0 data.frame list
## macr_extra_drugnodrug 0 data.frame list
upset_plots_both <- upsetr_sig(
both_sig, both = TRUE,
contrasts = c("macr_z23sb_vs_z22sb", "macr_z23nosb_vs_z22nosb",
"u937_z23sb_vs_z22sb", "u937_z23nosb_vs_z22nosb"))
upset_plots_both$both
Compare DE
results from macrophages and U937 samples
Looking a bit more closely at these, I think the u937 data is too
sparse to effectively compare.
macr_u937_comparison <- compare_de_results(hs_macr_table, u937_table)
## Error: object 'u937_table' not found
macr_u937_comparison$lfc_heat
## Error in eval(expr, envir, enclos): object 'macr_u937_comparison' not found
macr_u937_venns <- compare_significant_contrasts(hs_macr_sig, second_sig_tables = u937_sig,
contrasts = "z23sb_vs_z23nosb")
## Error: object 'u937_sig' not found
## Error in eval(expr, envir, enclos): object 'macr_u937_venns' not found
macr_u937_venns$down_plot
## Error in eval(expr, envir, enclos): object 'macr_u937_venns' not found
macr_u937_venns_v2 <- compare_significant_contrasts(hs_macr_sig, second_sig_tables = u937_sig,
contrasts = "z22sb_vs_z22nosb")
## Error in compare_significant_contrasts(hs_macr_sig, second_sig_tables = u937_sig, : object 'u937_sig' not found
macr_u937_venns_v2$up_plot
## Error in eval(expr, envir, enclos): object 'macr_u937_venns_v2' not found
macr_u937_venns_v2$down_plot
## Error in eval(expr, envir, enclos): object 'macr_u937_venns_v2' not found
macr_u937_venns_v3 <- compare_significant_contrasts(hs_macr_sig, second_sig_tables = u937_sig,
contrasts = "sb_vs_uninf")
## Error in compare_significant_contrasts(hs_macr_sig, second_sig_tables = u937_sig, : object 'u937_sig' not found
macr_u937_venns_v3$up_plot
## Error in eval(expr, envir, enclos): object 'macr_u937_venns_v3' not found
macr_u937_venns_v3$down_plot
## Error in eval(expr, envir, enclos): object 'macr_u937_venns_v3' not found
Compare
macrophage/u937 with respect to z2.3/z2.2
comparison_df <- merge(hs_macr_table[["data"]][["z23sb_vs_z22sb"]],
u937_table[["data"]][["z23sb_vs_z22sb"]],
by = "row.names")
## Error in h(simpleError(msg, call)): error in evaluating the argument 'y' in selecting a method for function 'merge': object 'u937_table' not found
macru937_z23z22_plot <- plot_linear_scatter(comparison_df[, c("deseq_logfc.x", "deseq_logfc.y")])
## Error in is.data.frame(x): object 'comparison_df' not found
macru937_z23z22_plot$scatter
## Error in eval(expr, envir, enclos): object 'macru937_z23z22_plot' not found
comparison_df <- merge(hs_macr_table[["data"]][["z23nosb_vs_z22nosb"]],
u937_table[["data"]][["z23nosb_vs_z22nosb"]],
by = "row.names")
## Error in h(simpleError(msg, call)): error in evaluating the argument 'y' in selecting a method for function 'merge': object 'u937_table' not found
macru937_z23z22_plot <- plot_linear_scatter(comparison_df[, c("deseq_logfc.x", "deseq_logfc.y")])
## Error in is.data.frame(x): object 'comparison_df' not found
macru937_z23z22_plot$scatter
## Error in eval(expr, envir, enclos): object 'macru937_z23z22_plot' not found
Add donor to the
contrasts, no sva
no_power_fact <- paste0(pData(hs_macr)[["donor"]], "_",
pData(hs_macr)[["condition"]])
table(pData(hs_macr)[["donor"]])
##
## d01 d02 d09 d81
## 13 14 13 14
## no_power_fact
## d01_inf_z22 d01_inf_z23 d01_infsb_z22 d01_infsb_z23 d01_uninf_none
## 2 3 3 3 1
## d01_uninfsb_none d02_inf_z22 d02_inf_z23 d02_infsb_z22 d02_infsb_z23
## 1 3 3 3 3
## d02_uninf_none d02_uninfsb_none d09_inf_z22 d09_inf_z23 d09_infsb_z22
## 1 1 3 3 3
## d09_infsb_z23 d09_uninf_none d09_uninfsb_none d81_inf_z22 d81_inf_z23
## 2 1 1 3 3
## d81_infsb_z22 d81_infsb_z23 d81_uninf_none d81_uninfsb_none
## 3 3 1 1
hs_nopower <- set_expt_conditions(hs_macr, fact = no_power_fact)
##
## d01_inf_z22 d01_inf_z23 d01_infsb_z22 d01_infsb_z23 d01_uninf_none
## 2 3 3 3 1
## d01_uninfsb_none d02_inf_z22 d02_inf_z23 d02_infsb_z22 d02_infsb_z23
## 1 3 3 3 3
## d02_uninf_none d02_uninfsb_none d09_inf_z22 d09_inf_z23 d09_infsb_z22
## 1 1 3 3 3
## d09_infsb_z23 d09_uninf_none d09_uninfsb_none d81_inf_z22 d81_inf_z23
## 2 1 1 3 3
## d81_infsb_z22 d81_infsb_z23 d81_uninf_none d81_uninfsb_none
## 3 3 1 1
hs_nopower <- subset_expt(hs_nopower, subset="macrophagezymodeme!='none'")
## subset_expt(): There were 54, now there are 46 samples.
hs_nopower_nosva_de <- all_pairwise(hs_nopower, model_batch = FALSE, filter = TRUE)
##
## d01_inf_z22 d01_inf_z23 d01_infsb_z22 d01_infsb_z23 d02_inf_z22 d02_inf_z23
## 2 3 3 3 3 3
## d02_infsb_z22 d02_infsb_z23 d09_inf_z22 d09_inf_z23 d09_infsb_z22 d09_infsb_z23
## 3 3 3 3 3 2
## d81_inf_z22 d81_inf_z23 d81_infsb_z22 d81_infsb_z23
## 3 3 3 3
nopower_keepers <- list(
"d01_zymo" = c("d01infz23", "d01infz22"),
"d01_sbzymo" = c("d01infsbz23", "d01infsbz22"),
"d02_zymo" = c("d02infz23", "d02infz22"),
"d02_sbzymo" = c("d02infsbz23", "d02infsbz22"),
"d09_zymo" = c("d09infz23", "d09infz22"),
"d09_sbzymo" = c("d09infsbz23", "d09infsbz22"),
"d81_zymo" = c("d81infz23", "d81infz22"),
"d81_sbzymo" = c("d81infsbz23", "d81infsbz22"))
hs_nopower_nosva_table <- combine_de_tables(
hs_nopower_nosva_de, keepers = nopower_keepers,
excel = glue("analyses/macrophage_de/hs_nopower_table-v{ver}.xlsx"))
## Error in names(x) <- value: 'names' attribute [8] must be the same length as the vector [6]
## extra_contrasts = extra)
hs_nopower_nosva_sig <- extract_significant_genes(
hs_nopower_nosva_table,
excel = glue("analyses/macrophage_de/hs_nopower_nosva_sig-v{ver}.xlsx"))
## Error in extract_significant_genes(hs_nopower_nosva_table, excel = glue("analyses/macrophage_de/hs_nopower_nosva_sig-v{ver}.xlsx")): object 'hs_nopower_nosva_table' not found
d01d02_zymo_nosva_comp <- merge(hs_nopower_nosva_table[["data"]][["d01_zymo"]],
hs_nopower_nosva_table[["data"]][["d02_zymo"]],
by="row.names")
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'merge': object 'hs_nopower_nosva_table' not found
d0102_zymo_nosva_plot <- plot_linear_scatter(d01d02_zymo_nosva_comp[, c("deseq_logfc.x", "deseq_logfc.y")])
## Error in is.data.frame(x): object 'd01d02_zymo_nosva_comp' not found
d0102_zymo_nosva_plot$scatter
## Error in eval(expr, envir, enclos): object 'd0102_zymo_nosva_plot' not found
d0102_zymo_nosva_plot$correlation
## Error in eval(expr, envir, enclos): object 'd0102_zymo_nosva_plot' not found
d0102_zymo_nosva_plot$lm_rsq
## Error in eval(expr, envir, enclos): object 'd0102_zymo_nosva_plot' not found
d09d81_zymo_nosva_comp <- merge(hs_nopower_nosva_table[["data"]][["d09_zymo"]],
hs_nopower_nosva_table[["data"]][["d81_zymo"]],
by="row.names")
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'merge': object 'hs_nopower_nosva_table' not found
d0981_zymo_nosva_plot <- plot_linear_scatter(d09d81_zymo_nosva_comp[, c("deseq_logfc.x", "deseq_logfc.y")])
## Error in is.data.frame(x): object 'd09d81_zymo_nosva_comp' not found
d0981_zymo_nosva_plot$scatter
## Error in eval(expr, envir, enclos): object 'd0981_zymo_nosva_plot' not found
d0981_zymo_nosva_plot$correlation
## Error in eval(expr, envir, enclos): object 'd0981_zymo_nosva_plot' not found
d0981_zymo_nosva_plot$lm_rsq
## Error in eval(expr, envir, enclos): object 'd0981_zymo_nosva_plot' not found
d01d81_zymo_nosva_comp <- merge(hs_nopower_nosva_table[["data"]][["d01_zymo"]],
hs_nopower_nosva_table[["data"]][["d81_zymo"]],
by="row.names")
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'merge': object 'hs_nopower_nosva_table' not found
d0181_zymo_nosva_plot <- plot_linear_scatter(d01d81_zymo_nosva_comp[, c("deseq_logfc.x", "deseq_logfc.y")])
## Error in is.data.frame(x): object 'd01d81_zymo_nosva_comp' not found
d0181_zymo_nosva_plot$scatter
## Error in eval(expr, envir, enclos): object 'd0181_zymo_nosva_plot' not found
d0181_zymo_nosva_plot$correlation
## Error in eval(expr, envir, enclos): object 'd0181_zymo_nosva_plot' not found
d0181_zymo_nosva_plot$lm_rsq
## Error in eval(expr, envir, enclos): object 'd0181_zymo_nosva_plot' not found
upset_plots_nosva <- upsetr_sig(hs_nopower_nosva_sig, both=TRUE,
contrasts=c("d01_zymo", "d02_zymo", "d09_zymo", "d81_zymo"))
## Error in upsetr_sig(hs_nopower_nosva_sig, both = TRUE, contrasts = c("d01_zymo", : object 'hs_nopower_nosva_sig' not found
## Error in eval(expr, envir, enclos): object 'upset_plots_nosva' not found
## Error in eval(expr, envir, enclos): object 'upset_plots_nosva' not found
## Error in eval(expr, envir, enclos): object 'upset_plots_nosva' not found
## The 7th element in the both groups list is the set shared among all donors.
## I don't feel like writing out x:y:z:a
groups <- upset_plots_nosva[["both_groups"]]
## Error in eval(expr, envir, enclos): object 'upset_plots_nosva' not found
shared_genes <- attr(groups, "elements")[groups[[7]]] %>%
gsub(pattern = "^gene:", replacement = "")
## Error in (function (cond) : error in evaluating the argument 'x' in selecting a method for function 'gsub': subscript out of bounds
shared_gp <- simple_gprofiler(shared_genes)
## No results to show
## Please make sure that the organism is correct or set significant = FALSE
## No results to show
## Please make sure that the organism is correct or set significant = FALSE
shared_gp$pvalue_plots$MF
shared_gp$pvalue_plots$BP
shared_gp$pvalue_plots$REAC
shared_gp$pvalue_plots$WP
Add donor to the
contrasts, sva
hs_nopower_sva_de <- all_pairwise(hs_nopower, model_batch = "svaseq", filter = TRUE)
##
## d01_inf_z22 d01_inf_z23 d01_infsb_z22 d01_infsb_z23 d02_inf_z22 d02_inf_z23
## 2 3 3 3 3 3
## d02_infsb_z22 d02_infsb_z23 d09_inf_z22 d09_inf_z23 d09_infsb_z22 d09_infsb_z23
## 3 3 3 3 3 2
## d81_inf_z22 d81_inf_z23 d81_infsb_z22 d81_infsb_z23
## 3 3 3 3
## Removing 0 low-count genes (11720 remaining).
## Setting 2174 low elements to zero.
## transform_counts: Found 2174 values equal to 0, adding 1 to the matrix.
nopower_keepers <- list(
"d01_zymo" = c("d01infz23", "d01infz22"),
"d01_sbzymo" = c("d01infsbz23", "d01infsbz22"),
"d02_zymo" = c("d02infz23", "d02infz22"),
"d02_sbzymo" = c("d02infsbz23", "d02infsbz22"),
"d09_zymo" = c("d09infz23", "d09infz22"),
"d09_sbzymo" = c("d09infsbz23", "d09infsbz22"),
"d81_zymo" = c("d81infz23", "d81infz22"),
"d81_sbzymo" = c("d81infsbz23", "d81infsbz22"))
hs_nopower_sva_table <- combine_de_tables(
hs_nopower_sva_de, keepers = nopower_keepers,
excel = glue("analyses/macrophage_de/hs_nopower_table-v{ver}.xlsx"))
## Error in names(x) <- value: 'names' attribute [8] must be the same length as the vector [6]
## extra_contrasts = extra)
hs_nopower_sva_sig <- extract_significant_genes(
hs_nopower_sva_table,
excel = glue("analyses/macrophage_de/hs_nopower_sva_sig-v{ver}.xlsx"))
## Error in extract_significant_genes(hs_nopower_sva_table, excel = glue("analyses/macrophage_de/hs_nopower_sva_sig-v{ver}.xlsx")): object 'hs_nopower_sva_table' not found
d01d02_zymo_sva_comp <- merge(hs_nopower_sva_table[["data"]][["d01_zymo"]],
hs_nopower_sva_table[["data"]][["d02_zymo"]],
by="row.names")
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'merge': object 'hs_nopower_sva_table' not found
d0102_zymo_sva_plot <- plot_linear_scatter(d01d02_zymo_sva_comp[, c("deseq_logfc.x", "deseq_logfc.y")])
## Error in is.data.frame(x): object 'd01d02_zymo_sva_comp' not found
d0102_zymo_sva_plot$scatter
## Error in eval(expr, envir, enclos): object 'd0102_zymo_sva_plot' not found
d0102_zymo_sva_plot$correlation
## Error in eval(expr, envir, enclos): object 'd0102_zymo_sva_plot' not found
d0102_zymo_sva_plot$lm_rsq
## Error in eval(expr, envir, enclos): object 'd0102_zymo_sva_plot' not found
d09d81_zymo_sva_comp <- merge(hs_nopower_sva_table[["data"]][["d09_zymo"]],
hs_nopower_sva_table[["data"]][["d81_zymo"]],
by="row.names")
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'merge': object 'hs_nopower_sva_table' not found
d0981_zymo_sva_plot <- plot_linear_scatter(d09d81_zymo_sva_comp[, c("deseq_logfc.x", "deseq_logfc.y")])
## Error in is.data.frame(x): object 'd09d81_zymo_sva_comp' not found
d0981_zymo_sva_plot$scatter
## Error in eval(expr, envir, enclos): object 'd0981_zymo_sva_plot' not found
d0981_zymo_sva_plot$correlation
## Error in eval(expr, envir, enclos): object 'd0981_zymo_sva_plot' not found
d0981_zymo_sva_plot$lm_rsq
## Error in eval(expr, envir, enclos): object 'd0981_zymo_sva_plot' not found
d01d81_zymo_sva_comp <- merge(hs_nopower_sva_table[["data"]][["d01_zymo"]],
hs_nopower_sva_table[["data"]][["d81_zymo"]],
by="row.names")
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'merge': object 'hs_nopower_sva_table' not found
d0181_zymo_sva_plot <- plot_linear_scatter(d01d81_zymo_sva_comp[, c("deseq_logfc.x", "deseq_logfc.y")])
## Error in is.data.frame(x): object 'd01d81_zymo_sva_comp' not found
d0181_zymo_sva_plot$scatter
## Error in eval(expr, envir, enclos): object 'd0181_zymo_sva_plot' not found
d0181_zymo_sva_plot$correlation
## Error in eval(expr, envir, enclos): object 'd0181_zymo_sva_plot' not found
d0181_zymo_sva_plot$lm_rsq
## Error in eval(expr, envir, enclos): object 'd0181_zymo_sva_plot' not found
upset_plots_sva <- upsetr_sig(hs_nopower_sva_sig, both=TRUE,
contrasts=c("d01_zymo", "d02_zymo", "d09_zymo", "d81_zymo"))
## Error in upsetr_sig(hs_nopower_sva_sig, both = TRUE, contrasts = c("d01_zymo", : object 'hs_nopower_sva_sig' not found
## Error in eval(expr, envir, enclos): object 'upset_plots_sva' not found
## Error in eval(expr, envir, enclos): object 'upset_plots_sva' not found
## Error in eval(expr, envir, enclos): object 'upset_plots_sva' not found
## The 7th element in the both groups list is the set shared among all donors.
## I don't feel like writing out x:y:z:a
groups <- upset_plots_sva[["both_groups"]]
## Error in eval(expr, envir, enclos): object 'upset_plots_sva' not found
shared_genes <- attr(groups, "elements")[groups[[7]]] %>%
gsub(pattern = "^gene:", replacement = "")
## Error in (function (cond) : error in evaluating the argument 'x' in selecting a method for function 'gsub': subscript out of bounds
shared_gp <- simple_gprofiler(shared_genes)
## No results to show
## Please make sure that the organism is correct or set significant = FALSE
## No results to show
## Please make sure that the organism is correct or set significant = FALSE
shared_gp$pvalue_plots$MF
shared_gp$pvalue_plots$BP
shared_gp$pvalue_plots$REAC
shared_gp$pvalue_plots$WP
Donor
comparison
hs_donors <- set_expt_conditions(hs_macr, fact = "donor")
##
## d01 d02 d09 d81
## 13 14 13 14
donor_de <- all_pairwise(hs_donors, model_batch="svaseq", filter=TRUE)
##
## d01 d02 d09 d81
## 13 14 13 14
## Removing 0 low-count genes (11756 remaining).
## Setting 1225 low elements to zero.
## transform_counts: Found 1225 values equal to 0, adding 1 to the matrix.
donor_table <- combine_de_tables(
donor_de,
excel=glue("analyses/macrophage_de/donor_tables-v{ver}.xlsx"))
## Error in names(x) <- value: 'names' attribute [8] must be the same length as the vector [6]
donor_sig <- extract_significant_genes(
donor_table,
excel = glue("analyses/macrophage_de/donor_sig-v{ver}.xlsx"))
## Error in extract_significant_genes(donor_table, excel = glue("analyses/macrophage_de/donor_sig-v{ver}.xlsx")): object 'donor_table' not found
Primary query
contrasts
The final contrast in this list is interesting because it depends on
the extra contrasts applied to the all_pairwise() above. In my way of
thinking, the primary comparisons to consider are either cross-drug or
cross-strain, but not both. However I think in at least a few instances
Olga is interested in strain+drug / uninfected+nodrug.
Write contrast
results
Now let us write out the xlsx file containing the above contrasts.
The file with the suffix _table-version will therefore contain all genes
and the file with the suffix _sig-version will contain only those deemed
significant via our default criteria of DESeq2 |logFC| >= 1.0 and
adjusted p-value <= 0.05.