ipRGC bulk RNAseq – All Samples – This time with feeling (transcript-level abundances)!

1 Introduction

This document is intended to repeat with minimal modification the work performed by Theresa when analyzing the RNASeq IPRGC data. This is being performed in an apptainer (ugh I hate that name) container and so I will likely remove some of the libraries invoked in it.

Not all of these are available in the container, and I do not think I have the version of Theresa’s helper_functions.R which was used for this. I did find a version of it that appears to work without any problems, though.

At the end of this document, I am going to try to re-implement it with some minor changes to see if I can more easily reconcile the results to the tables provided by Colenso.

library("hpgltools")
library("ggplot2")
library("patchwork")
library("plotly")
library("RColorBrewer")
library("DT")
library("dplyr")
library("openxlsx")
library("ComplexHeatmap")
library("circlize")
library("gprofiler2")
source("/data/helper_functions.R")
library("gplots")
library("ggrepel")

#setwd("~/scratch_02/reesyxan/mmusculus_iprgc_2019")

2 M. musculus

2.1 Annotations

I am using mm38_100.

mm_annot <- load_biomart_annotations(species = "mmusculus")

## The biomart annotations file already exists, loading from it.

mm_annot <- mm_annot[["annotation"]]
mm_annot[["txid"]] <- paste0(mm_annot[["ensembl_transcript_id"]], ".", mm_annot[["version"]])
rownames(mm_annot) <- make.names(mm_annot[["ensembl_gene_id"]], unique=TRUE)

mm_tx_annot <- mm_annot
## yeah, the ID tag in my gff file has transcript: as a prefix, so the featureCounts outputs
## includes it...  I will fix that inpost.
rownames(mm_tx_annot) <- make.names(mm_tx_annot[["ensembl_transcript_id"]], unique = TRUE)

tx_gene_map <- mm_annot[, c("txid", "ensembl_gene_id")]

So, I now have a table of mouse annotations.

2.2 Hisat2 expressionset

mm38_hisat <- create_se("sample_sheets/samplesheet_bulkrnaseq_Opn4.xlsx",
                        gene_info = mm_tx_annot, file_column = "salmon_count_table")

## Reading the sample metadata.

## Checking the state of the condition column.

## Checking the state of the batch column.

## Checking the condition factor.

## The sample definitions comprises: 57 rows(samples) and 27 columns(metadata fields).

## There are 184622 rows with trailing numbers, removing them.

## Matched 184622 annotations and counts out of 184622 count rows and 916710 annotations.

## The final summarized experiment has 184622 rows and 27 columns.

3 Color Set Up

color_condition <- mm38_hisat$initial_metadata$condition
color_condition <- factor(color_condition, levels = c("Retina_het", "Retina_ko",
                                                      "SCN_het", "SCN_ko", "dLGN_het", "dLGN_ko"))

color_values <- c( "#00BFC4",  "#9AD2D3", "#F8766D", "#FBADA7", "#7CAE00", "#D3EE8F")
col_fun = colorRamp2(c(0, 3, 6), c("blue", "pink", "red"))

4 Mapped Counts Per Sample Plot

library(forcats)
mm38_hisat <- set_conditions(mm38_hisat, fact = "batch")

## The numbers of samples by condition are:

## 
##  1  2  3 
## 18 18 21

## Setting colors with no pre-defined colors, using the Dark2 palette.

libsize_df <- as.data.frame(plot_libsize(mm38_hisat)$table)
rownames(libsize_df) <- mm38_hisat$design$sampleid
libsize_df$samplenum <- paste0(mm38_hisat$design[rownames(libsize_df), "samplenumb"], "_",
                               mm38_hisat$design[rownames(libsize_df), "batch"])
libsize_df$sample_id <-  rownames(mm38_hisat$design)
libsize_df <- libsize_df[order(libsize_df$sum),]

ggplot(libsize_df, aes(x = fct_rev(fct_infreq(samplenum, sum)), y = sum, fill = condition)) +
  geom_bar(stat = "identity", color = "black") +
  theme_classic() +
  #theme(axis.text.x= element_text(fct_rev(fct_infreq(libsize_df$samplenum, libsize_df$sum))), #remove x axis labels
  #    axis.ticks.x=element_blank()) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  scale_fill_manual(values = c("#60449b", "#009f6f", "orange")) +
  xlab("Sample") +
  ylab("Sum of Mapped Gene Counts") +
  labs(fill = "")

plot_libsize(mm38_hisat)

## Library sizes of 57 samples, 
## ranging from 3,037,266 to 19,545,726.

5 Normalize Expression

First we will filter out the genes which have low counts across all samples. Then I will normalize both without and with SVA correction.

mm38_filt <- normalize(mm38_hisat, filter = TRUE)

## Removing 128079 low-count genes (56543 remaining).

des <- colData(mm38_filt)
new_factor <- paste0(des[["location"]], "_", des[["genotype"]])
colData(mm38_hisat)[["location_genotype"]] <- new_factor

mm38_norm <- normalize(mm38_hisat, filter = TRUE, convert = "cpm",
                       norm = "quant", transform="log2")

## Removing 128079 low-count genes (56543 remaining).

## transform_counts: Found 1386506 values equal to 0, adding 1 to the matrix.

mm38_norm <- set_conditions(mm38_norm, fact = "location_genotype")

## The numbers of samples by condition are:

## 
## Retina_het  Retina_ko  Retina_wt    SCN_het     SCN_ko     SCN_wt   dLGN_het 
##         11          4          4         11          4          4         11 
##    dLGN_ko    dLGN_wt 
##          4          4

## Setting colors with no pre-defined colors, using the Dark2 palette.

6 Melanopsin Expression

Here is a plot of the Opn4 expression across all retina samples split by batch. The 4 quadrants in the plot are

het:1 – retina het samples in batch 1
het:2 – retina het samples in batch 2
ko:1 – retina ko samples in batch 1
ko:2 – retina ko samples in batch 2

This is a Tx comparison; I need to either drop this or hunt down the set of transcript IDs (which I did)

library(RColorBrewer)
colors <- brewer.pal(name = "Spectral", n = 2)
mm38_cpm <- normalize(mm38_filt, convert = "cpm")
df <- data.frame(exprs = assay(mm38_cpm)["ENSMUSG00000021799",],
                 genotype = as.factor(des[["genotype"]]),
                 location = as.factor(des[["location"]]),
                 sampleid = des[["sampleid"]],
                 samplenum = des[["samplenumb"]],
                 batch = des[["batch"]])

## Set postscript arguments
setEPS()
postscript("bulkRNAseq_Opn4Expression_Retina.eps")
cpm_df <- df |>
  filter(location == "Retina" | location == "retina")

## This is where I most miss ggstatsplot
## I am going to sit and write my own version of it.
ggplot(cpm_df, aes(x = samplenum, y = exprs, fill = genotype, group = genotype)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.8)) +
  xlab("Sample") +
  ylab("Opn4 Counts Per Million (CPMs)") +
  #scale_fill_manual(values = c("#FFC107", "#004D40", "orange"),
  #                  aesthetics = "fill") +
  ggtitle("Retina Melanopsin Expression") +
  facet_wrap(~genotype:batch, scales = "free") +
  ylim(c(0, 25)) +
  theme_classic()
dev.off()

The WT samples in batch 3 have about double the Opn4 expression as the het samples in batch 3 (and batch 2).

7 PCA

pca_norm <- plot_pca(mm38_norm, plot_labels = FALSE)
pca_norm[["plot"]]

## rownames(pca_norm$result$v)
labels <- des[["sampleid"]]
## labels <- paste0(mm38_norm$design$sampleid, "_", mm38_norm$design$samplenumb)
all_df <-  cbind(as.data.frame(pca_norm$result$v), labels)

shape_genotype <- mm38_norm$design$genotype
color_condition <- mm38_norm$condition

## filter(sample_id %in% c("S3_Het_Ret_S1", "S4_Het_Ret_S4", "S5_Het_Ret_S16", "S1_Het_Ret_S1", "S2_Het_Ret_S4", "S3_Het_Ret_S7", "S4_Het_Ret_S10", "S5_Het_Ret_S13")) |>

pca_norm_plot <-  ggplot(all_df, aes(x = PC1, y = PC2, label = labels)) +
  geom_point(aes(color = color_condition, shape = shape_genotype),
             size = 5) +
  #scale_color_manual(values = color_values,
  #                  aesthetics = "color") +
  theme_classic() +
  labs(color="Cell Compartment/Genotype")   +
  ggtitle("PCA with No Batch Correction")  +
  labs(x = paste0("PC1: ",  pca_norm$prop_var[1], "% variance"),
       y =  paste0("PC2: ",  pca_norm$prop_var[2], "% variance"))

#setEPS()                                             # Set postscript arguments
#postscript("bulkRNAseq_PCA_allsamples_labeledS5.eps")                           # Start graphics device driver
all_df[["labels"]] <- rownames(all_df)
pca_norm_plot + ggtitle("PCA Normalized Counts") +
geom_label_repel(data=filter(all_df,
                             labels %in% labels), #c('s5_het_dlgn', 's5_het_ret', 's5_het_scn')),
    nudge_x = 0.15,
    nudge_y = 0.15)

#dev.off()

7.1 Look at the first batch of sequencing

#batch 1 plot
all_df[rownames(all_df) %in% tolower(rownames(des[des[["batch"]] == "1",])),] |>
  ggplot(aes(x = PC1, y = PC2)) +
  geom_point(aes(color = color_condition[which(rownames(all_df) %in% tolower(rownames(des[des$batch == "1",])))],
                 shape = shape_genotype[which(rownames(all_df) %in% tolower(rownames(des[des$batch == "1",])))]),
             size = 5) +
  scale_color_manual(values = color_values,
                     aesthetics = "color") +
  theme_classic() +
  labs(color="Cell Compartment/Genotype", shape = "Genotype")   +
  ggtitle("PCA - Batch 1 RNAseq")  +
  labs(x = paste0("PC1: ",  pca_norm$prop_var[1], "% variance"),
       y =  paste0("PC2: ",  pca_norm$prop_var[2], "% variance"))

I think there is an easier way to do the above

b1 <- subset_se(mm38_norm, subset = "batch == '1'")
plot_pca(b1)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by Retina_het, Retina_ko, SCN_het, SCN_ko, dLGN_het, dLGN_ko
## Shapes are defined by 1.

b1_b2 <- subset_se(mm38_norm, subset = "batch == '1'|batch == '2'")
plot_pca(b1_b2)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by Retina_het, Retina_ko, SCN_het, SCN_ko, dLGN_het, dLGN_ko
## Shapes are defined by 1, 2.

#batch 1&2 plot
pca_batch1and2 <- all_df |>
  ggplot(aes(x = PC1, y = PC2)) +
  geom_point(aes(color = color_condition,
                 shape = shape_genotype),
             size = 5) +
  scale_color_manual(values = color_values,
                     aesthetics = "color") +
  theme_classic() +
  labs(color="Cell Compartment/Genotype", shape = "Genotype")   +
  ggtitle("PCA - Batch 1&2 RNAseq")  +
  labs(x = paste0("PC1: ",  pca_norm$prop_var[1], "% variance"),
       y =  paste0("PC2: ",  pca_norm$prop_var[2], "% variance"))


#setEPS()
#postscript(file = "rnaseq_pca-batch1and2.eps")
#pca_batch1and2
#dev.off()

In this pca plot, we see the retina samples cluster on the right in the PC1 direction and the brain samples cluster on the left in the PC 1 direction.

As we saw before, there is disparity between some of the hets. As it turns out, all of the new Het samples cluster with the KO samples as our one weirdo did before. So now it seems that it is quite easy to get some genetic drift that leads to these mice’s expression to behave in weird ways.

As we did before, let’s look at the Opn4 expression across all the retina samples to see if the Het samples that behave like KOs have lower Opn4 expression.

8 A closer look at the similar KOs and Hets

Let’s exclude the 2 Hets that are the most dissimilar from the bunch and look at a PCA plot of the rest of the samples.

## The following subset is using an older sample ID which I am unable to find anywhere.
## As a result, I rewrote that portion of Theresa's pca plot and labeled everything
## according to the extant naming scheme, I will now try and figure out which
## samples were removed.
## mm38_subset <- subset_se(
##   mm38_filt,
## subset = "sampleid != 'S3_Het_dLGN'&sampleid != 'S3_Het_Ret'&sampleid != 'S3_Het_SCN'&sampleid != 'S4_Het_dLGN'&sampleid != 'S4_Het_Ret'&sampleid != 'S4_Het_SCN'")
## S3_Het_dLGN_S2 is an outlier marked by Theresa
## S3_Het_Ret_S1 is an outlier
## S3_Het_SCN_S3 is an outlier
## S4_Het_dLGN_S11 is somewhat an outlier, I might be inclined to leave it.
## S4_Het_Ret_S4 is a huge outlier
## S4_Het_SCN_S6 is an outlier

## If I were choosing likely samples, I think I would go for:
##S4_Het_Ret_S4, S3_Het_Ret_S1, S3_Het_SCN_s3, S3_Het_dLGN_S2, S4_Het_dLGN_S5, S4_Het_SCN_S6

## Either way, I am adding a couple columns to the sample sheet with these calls for included samples.

mm38_subset <- subset_se(mm38_filt, subset = "included=='yes'")

mm38_subset <- normalize(mm38_subset, filter=TRUE, convert="cpm",
                         norm = "quant", transform="log2")

## Removing 2435 low-count genes (54108 remaining).

## transform_counts: Found 1107260 values equal to 0, adding 1 to the matrix.

mm38_subset <- set_conditions(mm38_subset, fact = "location_time")

## The numbers of samples by condition are:

## 
## Retina_het  Retina_ko  Retina_wt    SCN_het     SCN_ko     SCN_wt   dLGN_het 
##          9          4          4          9          4          4          9 
##    dLGN_ko    dLGN_wt 
##          4          4

## Setting colors with no pre-defined colors, using the Dark2 palette.

pca_norm_subset <- plot_pca(mm38_subset, plot_labels = FALSE)
pca_norm_subset

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by Retina_het, Retina_ko, Retina_wt, SCN_het, SCN_ko, SCN_wt, dLGN_het, dLGN_ko, dLGN_wt
## Shapes are defined by 1, 2, 3.

labels <- rownames(pca_norm_subset$result$v)

all_df <-  cbind(as.data.frame(pca_norm_subset$result$v), labels)
shape_genotype <- mm38_subset$design$genotype
color_condition <- pca_norm_subset$table$condition

pca_norm_plot <-  ggplot(pca_norm_subset$table, aes(x = PC1, y = PC2)) +
  geom_point(aes(color = condition),
             size = 5) +
  scale_color_manual(values = color_values,
                     aesthetics = "color") +
  theme_classic() +
  labs(color="Cell Compartment/Genotype")   +
  ggtitle("PCA with No Batch Correction")  +
  labs(x = paste0("PC1: ",  pca_norm_subset$prop_var[1], "% variance"),
       y =  paste0("PC2: ",  pca_norm_subset$prop_var[2], "% variance"))

pca_norm_plot <-  ggplot(all_df, aes(x = PC1, y = PC2, label = labels)) +
  geom_point(aes(color = color_condition, shape = shape_genotype),
             size = 5) +
  #scale_color_manual(values = color_values,
  #                  aesthetics = "color") +
  theme_classic() +
  labs(color="Cell Compartment/Genotype")   +
  ggtitle("PCA with No Batch Correction")  +
  labs(x = paste0("PC1: ",  pca_norm_subset$prop_var[1], "% variance"),
       y =  paste0("PC2: ",  pca_norm_subset$prop_var[2], "% variance"))

#setEPS()                                             # Set postscript arguments
#postscript("bulkRNAseq_PCA_allsamples_labeledS5.eps")                           # Start graphics device driver
pca_norm_plot + ggtitle("PCA Normalized Counts -- Subset") +
  geom_label_repel(data=filter(all_df,
                               labels %in% c('s4_het_scn_s12', 's4_het_ret', 's4_het_scn', 's5_het_dlgn')),
                   nudge_x = 0.15,
                   nudge_y = 0.1)

pca_norm_plot + ggtitle("PCA Normalized Counts")# +

pca_norm_plot

# geom_label_repel(data=filter(all_df,
#                               labels %in% labels), #c('s4_het_scn', 's5_het_ret', 's5_het_scn', 's5_het_dlgn')),
#                   nudge_x = 0.15,
#                   nudge_y = 0.15)
#dev.off()

There really doesn’t seem to be a discernible difference globally between the hets and kos within each location. We can do a DE with these samples between het/kos to quantify this.

9 DE of KOs versus Low expressor Hets

mm38_filt_sub <- subset_se(mm38_filt, subset = "sampleid != 'S3_Het_dLGN'&sampleid != 'S3_Het_Ret'&sampleid != 'S3_Het_SCN'&sampleid != 'S4_Het_dLGN' & sampleid != 'S4_Het_Ret' & sampleid != 'S4_Het_SCN'")

mm38_dlgn_subset <- subset_se(mm38_filt_sub, subset = "location_time == 'dLGN_het' | location_time == 'dLGN_ko'")
mm38_ret_subset <- subset_se(mm38_filt_sub, subset = "location_time == 'Retina_het' | location_time == 'Retina_ko'")
mm38_scn_subset <- subset_se(mm38_filt_sub, subset = "location_time == 'SCN_het' | location_time == 'SCN_ko'")

mm38_dlgn_subset <- normalize(mm38_dlgn_subset, filter = TRUE)

## Removing 16951 low-count genes (39592 remaining).

mm38_scn_subset <- normalize(mm38_scn_subset, filter = TRUE)

## Removing 16859 low-count genes (39684 remaining).

mm38_ret_subset <- normalize(mm38_ret_subset, filter = TRUE)

## Removing 18328 low-count genes (38215 remaining).

mm38_dlgn_subset <- set_batches(mm38_dlgn_subset, fact = "condition")

## The number of samples by batch are:

## 
## 1 2 3 
## 6 6 3

mm38_scn_subset <- set_batches(mm38_scn_subset, fact = "condition")

## The number of samples by batch are:

## 
## 1 2 3 
## 6 6 3

mm38_ret_subset <- set_batches(mm38_ret_subset, fact = "condition")

## The number of samples by batch are:

## 
## 1 2 3 
## 6 6 3

mm38_dlgn_subset <- set_conditions(mm38_dlgn_subset, fact = "location_time")

## The numbers of samples by condition are:

## 
## dLGN_het  dLGN_ko 
##       11        4

## Setting colors with no pre-defined colors, using the Dark2 palette.

mm38_scn_subset <- set_conditions(mm38_scn_subset, fact = "location_time")

## The numbers of samples by condition are:

## 
## SCN_het  SCN_ko 
##      11       4

## Setting colors with no pre-defined colors, using the Dark2 palette.

mm38_ret_subset <- set_conditions(mm38_ret_subset, fact = "location_time")

## The numbers of samples by condition are:

## 
## Retina_het  Retina_ko 
##         11          4

## Setting colors with no pre-defined colors, using the Dark2 palette.

test <- set_conditions(mm38_filt, fact = "location_time") |>
normalize(transform = "log2", convert = "cpm", batch = "sva")

## The numbers of samples by condition are:

## 
## Retina_het  Retina_ko  Retina_wt    SCN_het     SCN_ko     SCN_wt   dLGN_het 
##         11          4          4         11          4          4         11 
##    dLGN_ko    dLGN_wt 
##          4          4

## Setting colors with no pre-defined colors, using the Dark2 palette.

## transform_counts: Found 460060 values less than 0.

## transform_counts: Found 460060 values equal to 0, adding 1 to the matrix.

plot_pca(test)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by Retina_het, Retina_ko, Retina_wt, SCN_het, SCN_ko, SCN_wt, dLGN_het, dLGN_ko, dLGN_wt
## Shapes are defined by 1, 2, 3.

I do not think one should really do the DE with the data separate like this unless there is some compelling reason. I am going to repeat without this constraint and see how similar/different the results are. In addition, the current revision of all_pairwise() is much more comfortable with more complex models, even with sva; so that might make this more fun.

mm_de_dlgn <- all_pairwise(mm38_dlgn_subset, model_fstring = "~ 0 + condition",
                           model_svs = "svaseq", filter = TRUE)

## dLGN_het  dLGN_ko 
##       11        4

## Removing 0 low-count genes (39592 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## Basic step 0/3: Transforming data.

## Setting 219467 entries to zero.

## This received a matrix of SVs.

## Error in DESeqDataSet(se, design = design, ignoreRank) : 
##   some values in assay are not integers

## conditions
## dLGN_het  dLGN_ko 
##       11        4

## conditions
## dLGN_het  dLGN_ko 
##       11        4

## conditions
## dLGN_het  dLGN_ko 
##       11        4

## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The deseq results are not of class hpgltools::deseq_pairwise.

mm_de_dlgn

## A pairwise differential expression with results from: basic, deseq, ebseq, edger, limma, noiseq.

## The factors used are: dLGN_het, dLGN_ko.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 1 comparisons.

## The logFC agreement among the methods follows:

##                 dLGN___LGN
## basic_vs_dream   0.5048447
## basic_vs_ebseq   0.5902495
## basic_vs_edger   0.3869976
## basic_vs_limma   0.5131146
## basic_vs_noiseq  0.7710287
## dream_vs_ebseq   0.5027929
## dream_vs_edger   0.7223475
## dream_vs_limma   0.9804184
## dream_vs_noiseq  0.3396771
## ebseq_vs_edger   0.5814362
## ebseq_vs_limma   0.4727963
## ebseq_vs_noiseq  0.7726008
## edger_vs_limma   0.7157456
## edger_vs_noiseq  0.3700917
## limma_vs_noiseq  0.3317524

mm_table_dlgn <- combine_de_tables(mm_de_dlgn,
                                   excel = "39bulk_tables_tx/mm_de_dlgn_table.xlsx")

## Looking for subscript invalid names, end of extract_keepers.

## Adding the coefficient plot to dLGN_ko_vs_dLGN_het.
## Adding the coefficient plot to dLGN_ko_vs_dLGN_het.
## Adding the coefficient plot to dLGN_ko_vs_dLGN_het.
## Adding the coefficient plot to dLGN_ko_vs_dLGN_het.
## Adding the coefficient plot to dLGN_ko_vs_dLGN_het.
## Adding the coefficient plot to dLGN_ko_vs_dLGN_het.

mm_sig_dlgn <- extract_significant_genes(mm_table_dlgn,
                                         excel = "39bulk_tables_tx/mm_de_dlgn_sig.xlsx")

## Unable to discern the methods to test significance, according_to is NULL.

mm_de_ret <- all_pairwise(mm38_ret_subset, model_fstring = "~ 0 + condition",
                          model_svs = "svaseq", filter = TRUE)

## Retina_het  Retina_ko 
##         11          4

## Removing 0 low-count genes (38215 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## Basic step 0/3: Transforming data.

## Setting 165632 entries to zero.

## This received a matrix of SVs.

## Error in DESeqDataSet(se, design = design, ignoreRank) : 
##   some values in assay are not integers

## conditions
## Retina_het  Retina_ko 
##         11          4

## conditions
## Retina_het  Retina_ko 
##         11          4

## conditions
## Retina_het  Retina_ko 
##         11          4

## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The deseq results are not of class hpgltools::deseq_pairwise.

mm_de_ret

## A pairwise differential expression with results from: basic, deseq, ebseq, edger, limma, noiseq.

## The factors used are: Retina_het, Retina_ko.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 1 comparisons.

## The logFC agreement among the methods follows:

##                 Rtn_k_v_R_
## basic_vs_dream   0.4567153
## basic_vs_ebseq   0.5316853
## basic_vs_edger   0.3292782
## basic_vs_limma   0.4681497
## basic_vs_noiseq  0.5870014
## dream_vs_ebseq   0.4335788
## dream_vs_edger   0.6679701
## dream_vs_limma   0.9889980
## dream_vs_noiseq  0.4626451
## ebseq_vs_edger   0.6107728
## ebseq_vs_limma   0.4217901
## ebseq_vs_noiseq  0.7164748
## edger_vs_limma   0.6669577
## edger_vs_noiseq  0.5256508
## limma_vs_noiseq  0.4617358

mm_table_ret <- combine_de_tables(mm_de_ret,
                                  excel = "39bulk_tables_tx/mm_de_ret_table.xlsx")

## Looking for subscript invalid names, end of extract_keepers.

## Adding the coefficient plot to Retina_ko_vs_Retina_het.
## Adding the coefficient plot to Retina_ko_vs_Retina_het.
## Adding the coefficient plot to Retina_ko_vs_Retina_het.
## Adding the coefficient plot to Retina_ko_vs_Retina_het.
## Adding the coefficient plot to Retina_ko_vs_Retina_het.
## Adding the coefficient plot to Retina_ko_vs_Retina_het.

mm_sig_ret <- extract_significant_genes(mm_table_ret,
                                         excel = "39bulk_tables_tx/mm_de_ret_sig.xlsx")

## Unable to discern the methods to test significance, according_to is NULL.

mm_de_scn <- all_pairwise(mm38_scn_subset, model_fstring = "~ 0 + condition",
                          model_svs = "svaseq", filter = TRUE)

## SCN_het  SCN_ko 
##      11       4

## Removing 0 low-count genes (39684 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## Basic step 0/3: Transforming data.

## Setting 201034 entries to zero.

## This received a matrix of SVs.

## Error in DESeqDataSet(se, design = design, ignoreRank) : 
##   some values in assay are not integers

## conditions
## SCN_het  SCN_ko 
##      11       4

## conditions
## SCN_het  SCN_ko 
##      11       4

## conditions
## SCN_het  SCN_ko 
##      11       4

## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The deseq results are not of class hpgltools::deseq_pairwise.

mm_de_scn

## A pairwise differential expression with results from: basic, deseq, ebseq, edger, limma, noiseq.

## The factors used are: SCN_het, SCN_ko.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 1 comparisons.

## The logFC agreement among the methods follows:

##                 SCN___SCN_
## basic_vs_dream   0.5541243
## basic_vs_ebseq   0.5750793
## basic_vs_edger   0.4334210
## basic_vs_limma   0.5651717
## basic_vs_noiseq  0.7471217
## dream_vs_ebseq   0.5613151
## dream_vs_edger   0.7242821
## dream_vs_limma   0.9846958
## dream_vs_noiseq  0.3919830
## ebseq_vs_edger   0.6739467
## ebseq_vs_limma   0.5196392
## ebseq_vs_noiseq  0.7754524
## edger_vs_limma   0.7127799
## edger_vs_noiseq  0.4396494
## limma_vs_noiseq  0.3742456

mm_table_scn <- combine_de_tables(mm_de_scn, excel = "39bulk_tables_tx/mm_de_scn_table.xlsx")

## Looking for subscript invalid names, end of extract_keepers.

## Adding the coefficient plot to SCN_ko_vs_SCN_het.
## Adding the coefficient plot to SCN_ko_vs_SCN_het.
## Adding the coefficient plot to SCN_ko_vs_SCN_het.
## Adding the coefficient plot to SCN_ko_vs_SCN_het.
## Adding the coefficient plot to SCN_ko_vs_SCN_het.
## Adding the coefficient plot to SCN_ko_vs_SCN_het.

mm_sig_scn <- extract_significant_genes(mm_table_scn, according_to = "deseq",
                                        excel = "39bulk_tables_tx/mm_de_scn_sig.xlsx")

## Unable to discern the methods to test significance, according_to is NULL.

9.1 SCN

keepers <- names(mm_de_scn$deseq$all_tables)

## Error in `mm_de_scn$deseq$all_tables`:
## ! $ operator is invalid for atomic vectors

res_tbls_scn <- c()
for (tbl in keepers) {
  deseq_tbl <- get_de_table(mm_de_scn, table = tbl)
  res_tbls_scn[[tbl]] <- merge(deseq_tbl, hisat_annot, by = "row.names")
  fctr <- c(paste(gsub("\\_.*","",tbl), "\nEnriched"),
            paste(sub('.*\\_', '', tbl), "\nEnriched"))
  res_tbls_scn[[tbl]] <- set_sig_limma(res_tbls_scn[[tbl]], factors = fctr)
  res_tbls_scn[[paste0(tbl, "volc_plotly")]] <- volc_plot(res_tbls_scn[[tbl]],
                                                          title = paste("Volcano Plot:", tbl),
                                                          type = "plotly",
                                                          tbl = "limma",
                                                          gene_name_col = "uniprot_gn_symbol")
  res_tbls_scn[[paste0(tbl, "volc_ggplot")]] <- volc_plot(res_tbls_scn[[tbl]],
                                                          title = paste("Volcano Plot:", tbl),
                                                          tbl = "limma")
}

## Error:
## ! object 'keepers' not found

res_tbls_scn$SCNko_vs_SCNhetvolc_ggplot

## NULL

res_tbls_scn$SCNko_vs_SCNhetvolc_plotly

## NULL

res_tbls_scn$SCNko_vs_SCNhet |>
  filter(adj.P.Val <= 0.05) |>
  filter(abs(logFC) >= 1) |>
  select("ensembl_gene_id", "uniprot_gn_symbol", "Significance", "description", "baseMean", "logFC", "adj.P.Val") |>
  arrange(logFC) |>
  DT::datatable(rownames = FALSE)

## Error in `UseMethod()`:
## ! no applicable method for 'filter' applied to an object of class "NULL"

9.2 dLGN

keepers <- names(mm_de_dlgn$deseq$all_tables)

## Error in `mm_de_dlgn$deseq$all_tables`:
## ! $ operator is invalid for atomic vectors

res_tbls_dlgn <- c()
for (tbl in keepers) {
  deseq_tbl <- get_de_table(mm_de_dlgn, table = tbl)
  res_tbls_dlgn[[tbl]] <- merge(deseq_tbl, hisat_annot, by = "row.names")
  res_tbls_dlgn[[tbl]] <- set_sig_limma(res_tbls_dlgn[[tbl]], factors = c(paste(gsub("\\_.*","",tbl), "\nEnriched"),
                                                                          paste(sub('.*\\_', '', tbl), "\nEnriched")))
  res_tbls_dlgn[[paste0(tbl, "volc_plotly")]] <- volc_plot(res_tbls_dlgn[[tbl]],
                                                           title = paste("Volcano Plot:", tbl),
                                                           type = "plotly",
                                                           tbl = "limma",
                                                           gene_name_col = "uniprot_gn_symbol")
  res_tbls_dlgn[[paste0(tbl, "volc_ggplot")]] <- volc_plot(res_tbls_dlgn[[tbl]],
                                                           title = paste("Volcano Plot:", tbl),
                                                           tbl = "limma")
}

## Error:
## ! object 'keepers' not found

res_tbls_dlgn$dLGNko_vs_dLGNhetvolc_ggplot

## NULL

res_tbls_dlgn$dLGNko_vs_dLGNhetvolc_plotly

## NULL

res_tbls_dlgn$dLGNko_vs_dLGNhet |>
  filter(adj.P.Val <= 0.05) |>
  filter(abs(logFC) >= 1) |>
  select("ensembl_gene_id", "uniprot_gn_symbol", "Significance", "description", "baseMean", "logFC", "adj.P.Val") |>
  arrange(logFC) |>
  DT::datatable(rownames = FALSE)

## Error in `UseMethod()`:
## ! no applicable method for 'filter' applied to an object of class "NULL"

9.3 Retina

keepers <- names(mm_de_ret$deseq$all_tables)

## Error in `mm_de_ret$deseq$all_tables`:
## ! $ operator is invalid for atomic vectors

res_tbls_ret <- c()
for (tbl in keepers) {
  deseq_tbl <- get_de_table(mm_de_ret, table = tbl)
  res_tbls_ret[[tbl]] <- merge(deseq_tbl, hisat_annot, by = "row.names")
  res_tbls_ret[[tbl]] <- set_sig_limma(res_tbls_ret[[tbl]], factors = c(paste(gsub("\\_.*","",tbl), "\nEnriched"),
                                                                        paste(sub('.*\\_', '', tbl), "\nEnriched")))
  res_tbls_ret[[paste0(tbl, "volc_plotly")]] <- volc_plot(res_tbls_ret[[tbl]],
                                                          title = paste("Volcano Plot:", tbl),
                                                          type = "plotly",
                                                          tbl = "limma",
                                                          gene_name_col = "uniprot_gn_symbol")
  res_tbls_ret[[paste0(tbl, "volc_ggplot")]] <- volc_plot(res_tbls_ret[[tbl]],
                                                          title = paste("Volcano Plot:", tbl),
                                                          tbl = "limma")
}

## Error:
## ! object 'keepers' not found

res_tbls_ret$Retinako_vs_Retinahetvolc_ggplot

## NULL

res_tbls_ret$Retinako_vs_Retinahetvolc_plotly

## NULL

res_tbls_ret$Retinako_vs_Retinahet |>
  filter(adj.P.Val <= 0.05) |>
  filter(abs(logFC) >= 1) |>
  select("ensembl_gene_id", "uniprot_gn_symbol", "Significance", "description", "baseMean", "logFC", "adj.P.Val") |>
  arrange(logFC) |>
  DT::datatable(rownames = FALSE)

## Error in `UseMethod()`:
## ! no applicable method for 'filter' applied to an object of class "NULL"

10 DE of all samples with SVA batch correction for high v low expressors

mm38_dlgn_subset <- subset_se(
  mm38_filt, subset = "location_time == 'dLGN_het'| location_time == 'dLGN_ko'")
mm38_ret_subset <- subset_se(
  mm38_filt, subset = "location_time == 'Retina_het'| location_time == 'Retina_ko'")
mm38_scn_subset <- subset_se(
  mm38_filt, subset = "location_time == 'SCN_het'| location_time == 'SCN_ko'")

mm38_dlgn_subset <- normalize(mm38_dlgn_subset, filter = TRUE)

## Removing 16951 low-count genes (39592 remaining).

mm38_scn_subset <- normalize(mm38_scn_subset, filter = TRUE)

## Removing 16859 low-count genes (39684 remaining).

mm38_ret_subset <- normalize(mm38_ret_subset, filter = TRUE)

## Removing 18328 low-count genes (38215 remaining).

mm38_dlgn_subset <- set_batches(mm38_dlgn_subset, fact = "condition")

## The number of samples by batch are:

## 
## 1 2 3 
## 6 6 3

mm38_scn_subset <- set_batches(mm38_scn_subset, fact = "condition")

## The number of samples by batch are:

## 
## 1 2 3 
## 6 6 3

mm38_ret_subset <- set_batches(mm38_ret_subset, fact = "condition")

## The number of samples by batch are:

## 
## 1 2 3 
## 6 6 3

mm38_dlgn_subset <- set_conditions(mm38_dlgn_subset, fact = "location_time")

## The numbers of samples by condition are:

## 
## dLGN_het  dLGN_ko 
##       11        4

## Setting colors with no pre-defined colors, using the Dark2 palette.

mm38_scn_subset <- set_conditions(mm38_scn_subset, fact = "location_time")

## The numbers of samples by condition are:

## 
## SCN_het  SCN_ko 
##      11       4

## Setting colors with no pre-defined colors, using the Dark2 palette.

mm38_ret_subset <- set_conditions(mm38_ret_subset, fact = "location_time")

## The numbers of samples by condition are:

## 
## Retina_het  Retina_ko 
##         11          4

## Setting colors with no pre-defined colors, using the Dark2 palette.

mm_de_dlgn_sva <- all_pairwise(mm38_dlgn_subset, model_fstring = "~ 0 + condition",
                               model_svs = "svaseq", filter = TRUE)

## dLGN_het  dLGN_ko 
##       11        4

## Removing 0 low-count genes (39592 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## Basic step 0/3: Transforming data.

## Setting 219467 entries to zero.

## This received a matrix of SVs.

## Error in DESeqDataSet(se, design = design, ignoreRank) : 
##   some values in assay are not integers

## conditions
## dLGN_het  dLGN_ko 
##       11        4

## conditions
## dLGN_het  dLGN_ko 
##       11        4

## conditions
## dLGN_het  dLGN_ko 
##       11        4

## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The deseq results are not of class hpgltools::deseq_pairwise.

mm_de_ret_sva <- all_pairwise(mm38_ret_subset, model_fstring = "~ 0 + condition",
                               model_svs = "svaseq", filter = TRUE)

## Retina_het  Retina_ko 
##         11          4

## Removing 0 low-count genes (38215 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## Basic step 0/3: Transforming data.

## Setting 165632 entries to zero.

## This received a matrix of SVs.

## Error in DESeqDataSet(se, design = design, ignoreRank) : 
##   some values in assay are not integers

## conditions
## Retina_het  Retina_ko 
##         11          4

## conditions
## Retina_het  Retina_ko 
##         11          4

## conditions
## Retina_het  Retina_ko 
##         11          4

## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The deseq results are not of class hpgltools::deseq_pairwise.

mm_de_scn_sva <- all_pairwise(mm38_scn_subset, model_fstring = "~ 0 + condition",
                               model_svs = "svaseq", filter = TRUE)

## SCN_het  SCN_ko 
##      11       4

## Removing 0 low-count genes (39684 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## Basic step 0/3: Transforming data.

## Setting 201034 entries to zero.

## This received a matrix of SVs.

## Error in DESeqDataSet(se, design = design, ignoreRank) : 
##   some values in assay are not integers

## conditions
## SCN_het  SCN_ko 
##      11       4

## conditions
## SCN_het  SCN_ko 
##      11       4

## conditions
## SCN_het  SCN_ko 
##      11       4

## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The deseq results are not of class hpgltools::deseq_pairwise.

mm_de_scn_sva

## A pairwise differential expression with results from: basic, deseq, ebseq, edger, limma, noiseq.

## The factors used are: SCN_het, SCN_ko.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 1 comparisons.

## The logFC agreement among the methods follows:

##                 SCN___SCN_
## basic_vs_dream   0.5541243
## basic_vs_ebseq   0.5750793
## basic_vs_edger   0.4334210
## basic_vs_limma   0.5651717
## basic_vs_noiseq  0.7471217
## dream_vs_ebseq   0.5613151
## dream_vs_edger   0.7242821
## dream_vs_limma   0.9846958
## dream_vs_noiseq  0.3919830
## ebseq_vs_edger   0.6739467
## ebseq_vs_limma   0.5196392
## ebseq_vs_noiseq  0.7754524
## edger_vs_limma   0.7127799
## edger_vs_noiseq  0.4396494
## limma_vs_noiseq  0.3742456

mm_scn_sva_table <- combine_de_tables(
  mm_de_scn_sva, excel = "excel/mm_subset_scn_sva_table.xlsx")

## Looking for subscript invalid names, end of extract_keepers.

## Adding the coefficient plot to SCN_ko_vs_SCN_het.
## Adding the coefficient plot to SCN_ko_vs_SCN_het.
## Adding the coefficient plot to SCN_ko_vs_SCN_het.
## Adding the coefficient plot to SCN_ko_vs_SCN_het.
## Adding the coefficient plot to SCN_ko_vs_SCN_het.
## Adding the coefficient plot to SCN_ko_vs_SCN_het.

mm_scn_sva_table

## A set of combined differential expression results.

##               table deseq_sigup deseq_sigdown edger_sigup edger_sigdown
## 1 SCN_ko_vs_SCN_het           0             0         276           156
##   limma_sigup limma_sigdown
## 1          17             0

## Only  has information, cannot create an UpSet.

## Plot describing unique/shared genes in a differential expression table.

## NULL

mm_scn_sva_sig <- extract_significant_genes(
  mm_scn_sva_table, excel = "excel/mm_subset_scn_sva_sig.xlsx")

## Unable to discern the methods to test significance, according_to is NULL.

mm_scn_sva_sig

## A set of genes deemed significant according to .

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

## data frame with 0 columns and 0 rows

10.1 dLGN

10.1.1 PCA pre and post SVAseq

mm38_dlgn_subset_norm_batch <- normalize(mm38_dlgn_subset, convert = "cpm",
                                         transform = "log2", batch = "svaseq")

## transform_counts: Found 55001 values less than 0.

## transform_counts: Found 55001 values equal to 0, adding 1 to the matrix.

mm38_dlgn_subset_norm_nobatch <- normalize(mm38_dlgn_subset, convert = "cpm", norm = "quant",
                                           transform = "log2")

## transform_counts: Found 197525 values equal to 0, adding 1 to the matrix.

pca_mm38_dlgn_subset_norm_nobatch <- plot_pca(mm38_dlgn_subset_norm_nobatch)
pca_norm_plot_nobatch <-  ggplot(pca_mm38_dlgn_subset_norm_nobatch$table, aes(x = PC1, y = PC2, label = labels)) +
  geom_point(aes(color = condition),
             size = 5) +
  scale_color_manual(values = color_values,
                     aesthetics = "color") +
  theme_classic() +
  labs(color="Cell Compartment/Genotype")   +
  ggtitle("PCA with No Batch Correction")  +
  labs(x = paste0("PC1: ",  pca_mm38_dlgn_subset_norm_nobatch$prop_var[1], "% variance"),
       y =  paste0("PC2: ",  pca_mm38_dlgn_subset_norm_nobatch$prop_var[2], "% variance"))  +
  geom_label_repel(data=filter(pca_mm38_dlgn_subset_norm_nobatch$table,
                               labels %in% c('S3_Ht_dLGN', 'S4_Ht_dLGN')),
                   nudge_x = 0.1,
                   nudge_y = 0.1)


pca_mm38_dlgn_subset_norm_batch <- plot_pca(mm38_dlgn_subset_norm_batch)
pca_norm_plot_batch <-  ggplot(pca_mm38_dlgn_subset_norm_batch$table, aes(x = PC1, y = PC2, label = labels)) +
  geom_point(aes(color = condition),
             size = 5) +
  scale_color_manual(values = color_values,
                     aesthetics = "color") +
  theme_classic() +
  labs(color="Cell Compartment/Genotype")   +
  ggtitle("PCA with SVA Batch Correction")  +
  labs(x = paste0("PC1: ",  pca_mm38_dlgn_subset_norm_batch$prop_var[1], "% variance"),
       y =  paste0("PC2: ",  pca_mm38_dlgn_subset_norm_batch$prop_var[2], "% variance"))  +
  geom_label_repel(data=filter(pca_mm38_dlgn_subset_norm_batch$table,
                               labels %in% c('S3_Ht_dLGN', 'S4_Ht_dLGN')),
                   nudge_x = 0.1,
                   nudge_y = 0.1)


pca_norm_plot_nobatch

pca_norm_plot_batch

10.1.2 Volcano Plots ad DE table results

keepers <- names(mm_de_dlgn_sva$deseq$all_tables)

## Error in `mm_de_dlgn_sva$deseq$all_tables`:
## ! $ operator is invalid for atomic vectors

res_tbls_dlgn_sva <- c()
for (tbl in keepers) {
  deseq_tbl <- get_de_table(mm_de_dlgn_sva, table = tbl)
  res_tbls_dlgn_sva[[tbl]] <- merge(deseq_tbl, hisat_annot, by = "row.names")
  res_tbls_dlgn_sva[[tbl]] <- set_sig_limma(res_tbls_dlgn_sva[[tbl]], factors = c(paste(gsub("\\_.*","",tbl), "\nEnriched"),
                                                                                  paste(sub('.*\\_', '', tbl), "\nEnriched")))
  res_tbls_dlgn_sva[[paste0(tbl, "volc_plotly")]] <- volc_plot(res_tbls_dlgn_sva[[tbl]],
                                                               title = paste("Volcano Plot:", tbl),
                                                               type = "plotly",
                                                               tbl = "limma",
                                                               gene_name_col = "uniprot_gn_symbol")
  res_tbls_dlgn_sva[[paste0(tbl, "volc_ggplot")]] <- volc_plot(res_tbls_dlgn_sva[[tbl]],
                                                               title = paste("Volcano Plot:", tbl),
                                                               tbl = "limma")
}

## Error:
## ! object 'keepers' not found

res_tbls_dlgn_sva$dLGNko_vs_dLGNhetvolc_ggplot

## NULL

res_tbls_dlgn_sva$dLGNko_vs_dLGNhetvolc_plotly

## NULL

res_tbls_dlgn_sva$dLGNko_vs_dLGNhet |>
  filter(adj.P.Val <= 0.05) |>
  filter(abs(logFC) >= 1) |>
  select("ensembl_gene_id", "uniprot_gn_symbol", "Significance", "description", "baseMean", "logFC", "adj.P.Val") |>
  arrange(logFC) |>
  DT::datatable(rownames = FALSE)

## Error in `UseMethod()`:
## ! no applicable method for 'filter' applied to an object of class "NULL"

Samd11 and Dmrt3 are still both upregulated in dLGN KOs (as they were without including the high-expressors).

Pttg1, Sema3g, Ifih1, Dppa5a, and Prg4 are all upregulated in dLGN Hets including and excluding the high expressors. We do lose a few genes when including high expressors that were enriched in dLGN hets before.

10.2 SCN

10.2.1 PCA pre and post SVAseq

mm38_scn_subset_norm_batch <- normalize(mm38_scn_subset, convert = "cpm",
                                        transform = "log2", batch = "svaseq")

## transform_counts: Found 49815 values less than 0.

## transform_counts: Found 49815 values equal to 0, adding 1 to the matrix.

mm38_scn_subset_norm_nobatch <- normalize(mm38_scn_subset, convert = "cpm", norm = "quant",
                                          transform = "log2")

## transform_counts: Found 176602 values equal to 0, adding 1 to the matrix.

pca_mm38_scn_subset_norm_nobatch <- plot_pca(mm38_scn_subset_norm_nobatch)
pca_norm_plot_nobatch_scn <-  ggplot(pca_mm38_scn_subset_norm_nobatch$table,
                                     aes(x = PC1, y = PC2, label = labels)) +
  geom_point(aes(color = condition),
             size = 5) +
  scale_color_manual(values = color_values[5:6],
                     aesthetics = "color") +
  theme_classic() +
  labs(color="Cell Compartment/Genotype")   +
  ggtitle("PCA with No Batch Correction")  +
  labs(x = paste0("PC1: ",  pca_mm38_scn_subset_norm_nobatch$prop_var[1], "% variance"),
       y =  paste0("PC2: ",  pca_mm38_scn_subset_norm_nobatch$prop_var[2], "% variance"))  +
  geom_label_repel(data=filter(pca_mm38_scn_subset_norm_nobatch$table,
                               labels %in% c('S3_Het_SCN', 'S4_Het_SCN')),
                   nudge_x = 0.2,
                   nudge_y = 0.2)


pca_mm38_scn_subset_norm_batch <- plot_pca(mm38_scn_subset_norm_batch)
pca_norm_plot_batch_scn <-  ggplot(pca_mm38_scn_subset_norm_batch$table, aes(x = PC1, y = PC2, label = labels)) +
  geom_point(aes(color = condition),
             size = 5) +
  scale_color_manual(values = color_values[5:6],
                     aesthetics = "color") +
  theme_classic() +
  labs(color="Cell Compartment/Genotype")   +
  ggtitle("PCA with SVA Batch Correction")  +
  labs(x = paste0("PC1: ",  pca_mm38_scn_subset_norm_batch$prop_var[1], "% variance"),
       y =  paste0("PC2: ",  pca_mm38_scn_subset_norm_batch$prop_var[2], "% variance"))  +
  geom_label_repel(data=filter(pca_mm38_scn_subset_norm_batch$table,
                               labels %in% c('S3_Het_SCN', 'S4_Het_SCN')),
                   nudge_x = 0.2,
                   nudge_y = 0.2)


pca_norm_plot_nobatch_scn

pca_norm_plot_batch_scn

10.2.2 Volcano Plots ad DE table results

keepers <- names(mm_de_scn_sva$deseq$all_tables)

## Error in `mm_de_scn_sva$deseq$all_tables`:
## ! $ operator is invalid for atomic vectors

res_tbls_scn_sva <- c()
for (tbl in keepers) {
  deseq_tbl <- get_de_table(mm_de_scn_sva, table = tbl)
  res_tbls_scn_sva[[tbl]] <- merge(deseq_tbl, hisat_annot, by = "row.names")
  res_tbls_scn_sva[[tbl]] <- set_sig_limma(res_tbls_scn_sva[[tbl]], factors = c(paste(gsub("\\_.*","",tbl), "\nEnriched"),
                                                                                paste(sub('.*\\_', '', tbl), "\nEnriched")))
  res_tbls_scn_sva[[paste0(tbl, "volc_plotly")]] <- volc_plot(res_tbls_scn_sva[[tbl]],
                                                              title = paste("Volcano Plot:", tbl),
                                                              type = "plotly",
                                                              tbl = "limma",
                                                              gene_name_col = "uniprot_gn_symbol")
  res_tbls_scn_sva[[paste0(tbl, "volc_ggplot")]] <- volc_plot(res_tbls_scn_sva[[tbl]],
                                                              title = paste("Volcano Plot:", tbl),
                                                              tbl = "limma")
}

## Error:
## ! object 'keepers' not found

res_tbls_scn_sva$SCNko_vs_SCNhetvolc_ggplot

## NULL

res_tbls_scn_sva$SCNko_vs_SCNhetvolc_plotly

## NULL

res_tbls_scn_sva$SCNko_vs_SCNhet |>
  filter(adj.P.Val <= 0.05) |>
  filter(abs(logFC) >= 1) |>
  select("ensembl_gene_id", "uniprot_gn_symbol", "Significance", "description", "baseMean", "logFC", "adj.P.Val") |>
  arrange(logFC) |>
  DT::datatable(rownames = FALSE)

## Error in `UseMethod()`:
## ! no applicable method for 'filter' applied to an object of class "NULL"

10.3 Retina

10.3.1 PCA pre and post SVAseq

mm38_ret_subset_norm_batch <- normalize(mm38_ret_subset, convert = "cpm",
                                        transform = "log2", batch = "svaseq")

## transform_counts: Found 48664 values less than 0.

## transform_counts: Found 48664 values equal to 0, adding 1 to the matrix.

mm38_ret_subset_norm_nobatch <- normalize(mm38_ret_subset, convert = "cpm", norm = "quant",
                                          transform = "log2")

## transform_counts: Found 137516 values equal to 0, adding 1 to the matrix.

pca_mm38_ret_subset_norm_nobatch <- plot_pca(mm38_ret_subset_norm_nobatch)
pca_norm_plot_nobatch_ret <-  ggplot(pca_mm38_ret_subset_norm_nobatch$table, aes(x = PC1, y = PC2, label = labels)) +
  geom_point(aes(color = condition),
             size = 5) +
  scale_color_manual(values = color_values[3:4],
                     aesthetics = "color") +
  theme_classic() +
  labs(color="Cell Compartment/Genotype")   +
  ggtitle("PCA with No Batch Correction")  +
  labs(x = paste0("PC1: ",  pca_mm38_ret_subset_norm_nobatch$prop_var[1], "% variance"),
       y =  paste0("PC2: ",  pca_mm38_ret_subset_norm_nobatch$prop_var[2], "% variance"))  +
  geom_label_repel(data=filter(pca_mm38_ret_subset_norm_nobatch$table,
                               labels %in% c('S3_Het_Ret', 'S4_Het_Ret')),
                   nudge_x = 0.2,
                   nudge_y = 0.2)


pca_mm38_ret_subset_norm_batch <- plot_pca(mm38_ret_subset_norm_batch)
pca_norm_plot_batch_ret <-  ggplot(pca_mm38_ret_subset_norm_batch$table, aes(x = PC1, y = PC2, label = labels)) +
  geom_point(aes(color = condition),
             size = 5) +
  scale_color_manual(values = color_values[3:4],
                     aesthetics = "color") +
  theme_classic() +
  labs(color="Cell Compartment/Genotype")   +
  ggtitle("PCA with SVA Batch Correction")  +
  labs(x = paste0("PC1: ",  pca_mm38_ret_subset_norm_batch$prop_var[1], "% variance"),
       y =  paste0("PC2: ",  pca_mm38_ret_subset_norm_batch$prop_var[2], "% variance"))  +
  geom_label_repel(data=filter(pca_mm38_ret_subset_norm_batch$table,
                               labels %in% c('S3_Het_Ret', 'S4_Het_Ret')),
                   nudge_x = 0.2,
                   nudge_y = 0.2)


pca_norm_plot_nobatch_ret

pca_norm_plot_batch_ret

10.3.2 Volcano Plots ad DE table results

keepers <- names(mm_de_ret_sva$deseq$all_tables)

## Error in `mm_de_ret_sva$deseq$all_tables`:
## ! $ operator is invalid for atomic vectors

res_tbls_ret_sva <- c()
for (tbl in keepers) {
  deseq_tbl <- get_de_table(mm_de_ret_sva, table = tbl)
  res_tbls_ret_sva[[tbl]] <- merge(deseq_tbl, hisat_annot, by = "row.names")
  res_tbls_ret_sva[[tbl]] <- set_sig_limma(res_tbls_ret_sva[[tbl]], factors = c(paste(gsub("\\_.*","",tbl), "\nEnriched"),
                                                                                paste(sub('.*\\_', '', tbl), "\nEnriched")))
  res_tbls_ret_sva[[paste0(tbl, "volc_plotly")]] <- volc_plot(res_tbls_ret_sva[[tbl]],
                                                              title = paste("Volcano Plot:", tbl),
                                                              type = "plotly",
                                                              tbl = "limma",
                                                              gene_name_col = "uniprot_gn_symbol")
  res_tbls_ret_sva[[paste0(tbl, "volc_ggplot")]] <- volc_plot(res_tbls_ret_sva[[tbl]],
                                                              title = paste("Volcano Plot:", tbl),
                                                              tbl = "limma")
}

## Error:
## ! object 'keepers' not found

res_tbls_ret_sva$Retinako_vs_Retinahetvolc_ggplot

## NULL

res_tbls_ret_sva$Retinako_vs_Retinahetvolc_plotly

## NULL

res_tbls_ret_sva$Retinako_vs_Retinahet |>
  filter(adj.P.Val <= 0.05) |>
  filter(abs(logFC) >= 1) |>
  select("ensembl_gene_id", "uniprot_gn_symbol", "Significance", "description", "baseMean", "logFC", "adj.P.Val") |>
  arrange(logFC) |>
  DT::datatable(rownames = FALSE)

## Error in `UseMethod()`:
## ! no applicable method for 'filter' applied to an object of class "NULL"

11 Re-Implementation

Let us start with the full set of desired tables. In the above, the data was separated into ~ 9 separate chunks and I think it will benefit from being kept together – also I seem to recall expressing this to Theresa and that she changed the notebook to this kind of analysis. If I am correct, then it is possible she did so in a separate notebook which I have not found and that this will provide the fully reconciled set of results.

I have reasonable hope for this because the various plots etc early on look identical between this containerized notebook, and those notebooks I found in her working tree.

While I wait for these tables to come back, let us count the number of samples in the plot provided by Colenso so that I may doublecheck that nothing got filtered out or anything important like that.

8 pink samples, 17 green samples, 17 blue samples.

4 light pink, 4 dark pink 9 medium green, 4 light green, 4 dark green. 9 medium blue, 4 light blue, 4 dark blue.

With the above in mind, let us make a small change in how the data is handled and see if I cannot get numbers more closely aligned to the tables from Colenso.

## Here is a reminder of Theresa's filter:
## filter(sample_id %in% c("S3_Het_Ret_S1", "S4_Het_Ret_S4", "S5_Het_Ret_S16", "S1_Het_Ret_S1", "S2_Het_Ret_S4", "S3_Het_Ret_S7", "S4_Het_Ret_S10", "S5_Het_Ret_S13")) |>

mm_final <- subset_se(mm38_hisat, subset = "included=='yes'") |>
  set_conditions(fact = "location_genotype")

## The numbers of samples by condition are:

## 
## Retina_het  Retina_ko  Retina_wt    SCN_het     SCN_ko     SCN_wt   dLGN_het 
##          9          4          4          9          4          4          9 
##    dLGN_ko    dLGN_wt 
##          4          4

## Setting colors with no pre-defined colors, using the Dark2 palette.

mm_norm <- normalize(mm_final, transform = "log2", convert = "cpm", filter = TRUE, norm = "quant")

## Removing 130546 low-count genes (54076 remaining).

## transform_counts: Found 1106225 values equal to 0, adding 1 to the matrix.

plot_pca(mm_norm)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by Retina_het, Retina_ko, Retina_wt, SCN_het, SCN_ko, SCN_wt, dLGN_het, dLGN_ko, dLGN_wt
## Shapes are defined by 1, 2, 3.

mm_nb <- normalize(mm_final, transform = "log2", convert = "cpm", filter = TRUE, batch = "svaseq")

## Removing 130546 low-count genes (54076 remaining).

## transform_counts: Found 369880 values less than 0.

## transform_counts: Found 369880 values equal to 0, adding 1 to the matrix.

plot_pca(mm_nb)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by Retina_het, Retina_ko, Retina_wt, SCN_het, SCN_ko, SCN_wt, dLGN_het, dLGN_ko, dLGN_wt
## Shapes are defined by 1, 2, 3.

## Warning in MASS::cov.trob(data[, vars], wt = weight * nrow(data)): Probable
## convergence failure
## Warning in MASS::cov.trob(data[, vars], wt = weight * nrow(data)): Probable
## convergence failure

There remain a couple of weirdo samples, but things look quite a bit easier to interpret to my eyes.

final_keepers <- list(
  ## Three comparisons among the dLGN samples.
  "dlgn_kw" = c("dLGN_ko", "dLGN_wt"),
  "dlgn_kh" = c("dLGN_ko", "dLGN_het"),
  "dlgn_hw" = c("dLGN_het", "dLGN_wt"),
  ## Three comparisons among the SCN samples.
  "scn_kw" = c("SCN_ko", "SCN_wt"),
  "scn_kh" = c("SCN_ko", "SCN_het"),
  "scn_hw" = c("SCN_het", "SCN_wt"),
  ## and the retina samples
  "retina_kw" = c("Retina_ko", "Retina_wt"),
  "retina_kh" = c("Retina_ko", "Retina_het"),
  "retina_hw" = c("Retina_het", "Retina_wt")
)

genoloc_de_sva <- all_pairwise(mm_final, model_fstring = "~ 0 + condition", filter = TRUE,
                               model_svs = "svaseq", keepers = final_keepers)

## Retina_het  Retina_ko  Retina_wt    SCN_het     SCN_ko     SCN_wt   dLGN_het 
##          9          4          4          9          4          4          9 
##    dLGN_ko    dLGN_wt 
##          4          4

## Removing 130546 low-count genes (54076 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## Basic step 0/3: Transforming data.

## Setting 1297894 entries to zero.

## This received a matrix of SVs.

## Error in DESeqDataSet(se, design = design, ignoreRank) : 
##   some values in assay are not integers

## Warning in variancePartition::voomWithDreamWeights(counts = data, formula = model_fstring, : The maximum precision weight is 4.654586e+13, suggesting a poor smoothing fit
## on the mean-variance plot for large expression values. Such large weights can
## have unexpected effects downstream.  Consider examining the mean-variance plot
## and reducing the span parameter.

## conditions
## Retina_het  Retina_ko  Retina_wt    SCN_het     SCN_ko     SCN_wt   dLGN_het 
##          9          4          4          9          4          4          9 
##    dLGN_ko    dLGN_wt 
##          4          4

## conditions
## Retina_het  Retina_ko  Retina_wt    SCN_het     SCN_ko     SCN_wt   dLGN_het 
##          9          4          4          9          4          4          9 
##    dLGN_ko    dLGN_wt 
##          4          4

## conditions
## Retina_het  Retina_ko  Retina_wt    SCN_het     SCN_ko     SCN_wt   dLGN_het 
##          9          4          4          9          4          4          9 
##    dLGN_ko    dLGN_wt 
##          4          4

## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The deseq results are not of class hpgltools::deseq_pairwise.

genoloc_table_sva <- combine_de_tables(genoloc_de_sva, excel = "excel/genotype_location_tables.xlsx")

## Looking for subscript invalid names, end of extract_keepers.

## Adding the coefficient plot to dLGN_ko_vs_dLGN_wt.
## Adding the coefficient plot to dLGN_ko_vs_dLGN_wt.
## Adding the coefficient plot to dLGN_ko_vs_dLGN_wt.
## Adding the coefficient plot to dLGN_ko_vs_dLGN_wt.
## Adding the coefficient plot to dLGN_ko_vs_dLGN_wt.
## Adding the coefficient plot to dLGN_ko_vs_dLGN_wt.

## Adding the coefficient plot to dLGN_ko_vs_dLGN_het.
## Adding the coefficient plot to dLGN_ko_vs_dLGN_het.
## Adding the coefficient plot to dLGN_ko_vs_dLGN_het.
## Adding the coefficient plot to dLGN_ko_vs_dLGN_het.
## Adding the coefficient plot to dLGN_ko_vs_dLGN_het.
## Adding the coefficient plot to dLGN_ko_vs_dLGN_het.

## Adding the coefficient plot to dLGN_het_vs_dLGN_wt.
## Adding the coefficient plot to dLGN_het_vs_dLGN_wt.
## Adding the coefficient plot to dLGN_het_vs_dLGN_wt.
## Adding the coefficient plot to dLGN_het_vs_dLGN_wt.
## Adding the coefficient plot to dLGN_het_vs_dLGN_wt.
## Adding the coefficient plot to dLGN_het_vs_dLGN_wt.

## Adding the coefficient plot to SCN_ko_vs_SCN_wt.
## Adding the coefficient plot to SCN_ko_vs_SCN_wt.
## Adding the coefficient plot to SCN_ko_vs_SCN_wt.
## Adding the coefficient plot to SCN_ko_vs_SCN_wt.
## Adding the coefficient plot to SCN_ko_vs_SCN_wt.
## Adding the coefficient plot to SCN_ko_vs_SCN_wt.

## Adding the coefficient plot to SCN_ko_vs_SCN_het.
## Adding the coefficient plot to SCN_ko_vs_SCN_het.
## Adding the coefficient plot to SCN_ko_vs_SCN_het.
## Adding the coefficient plot to SCN_ko_vs_SCN_het.
## Adding the coefficient plot to SCN_ko_vs_SCN_het.
## Adding the coefficient plot to SCN_ko_vs_SCN_het.

## Adding the coefficient plot to SCN_het_vs_SCN_wt.
## Adding the coefficient plot to SCN_het_vs_SCN_wt.
## Adding the coefficient plot to SCN_het_vs_SCN_wt.
## Adding the coefficient plot to SCN_het_vs_SCN_wt.
## Adding the coefficient plot to SCN_het_vs_SCN_wt.
## Adding the coefficient plot to SCN_het_vs_SCN_wt.

## Adding the coefficient plot to Retina_ko_vs_Retina_wt.
## Adding the coefficient plot to Retina_ko_vs_Retina_wt.
## Adding the coefficient plot to Retina_ko_vs_Retina_wt.
## Adding the coefficient plot to Retina_ko_vs_Retina_wt.
## Adding the coefficient plot to Retina_ko_vs_Retina_wt.
## Adding the coefficient plot to Retina_ko_vs_Retina_wt.

## Adding the coefficient plot to Retina_ko_vs_Retina_het.
## Adding the coefficient plot to Retina_ko_vs_Retina_het.
## Adding the coefficient plot to Retina_ko_vs_Retina_het.
## Adding the coefficient plot to Retina_ko_vs_Retina_het.
## Adding the coefficient plot to Retina_ko_vs_Retina_het.
## Adding the coefficient plot to Retina_ko_vs_Retina_het.

## Adding the coefficient plot to Retina_het_vs_Retina_wt.
## Adding the coefficient plot to Retina_het_vs_Retina_wt.
## Adding the coefficient plot to Retina_het_vs_Retina_wt.
## Adding the coefficient plot to Retina_het_vs_Retina_wt.
## Adding the coefficient plot to Retina_het_vs_Retina_wt.
## Adding the coefficient plot to Retina_het_vs_Retina_wt.

genoloc_sig_sva <- extract_significant_genes(genoloc_table_sva, according_to = "deseq",
                                             excel = "excel/genotype_location_sig.xlsx")

## Unable to discern the methods to test significance, according_to is NULL.

location_keepers <- list(
  ## ko across locations
  "ko_dlgn_scn" = c("dLGN_ko", "SCN_ko"),
  "ko_dlgn_ret" = c("dLGN_ko", "Retina_ko"),
  "ko_scn_ret" = c("SCN_ko", "Retina_ko"),
  ## het across locations
  "het_dlgn_scn" = c("dLGN_het", "SCN_het"),
  "het_dlgn_ret" = c("dLGN_het", "Retina_het"),
  "het_scn_ret" = c("SCN_het", "Retina_het"),
  ## wt across locations
  "wt_dlgn_scn" = c("dLGN_wt", "SCN_wt"),
  "wt_dlgn_ret" = c("dLGN_wt", "Retina_wt"),
  "wt_scn_ret" = c("SCN_wt", "Retina_wt"))
location_comparisons <- all_pairwise(mm_final, model_fstring = "~ 0 + condition", filter = TRUE,
                               model_svs = "svaseq", keepers = location_keepers)

## Retina_het  Retina_ko  Retina_wt    SCN_het     SCN_ko     SCN_wt   dLGN_het 
##          9          4          4          9          4          4          9 
##    dLGN_ko    dLGN_wt 
##          4          4

## Removing 130546 low-count genes (54076 remaining).

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## Basic step 0/3: Transforming data.

## Setting 1297894 entries to zero.

## This received a matrix of SVs.

## Error in DESeqDataSet(se, design = design, ignoreRank) : 
##   some values in assay are not integers

## Warning in variancePartition::voomWithDreamWeights(counts = data, formula = model_fstring, : The maximum precision weight is 4.654586e+13, suggesting a poor smoothing fit
## on the mean-variance plot for large expression values. Such large weights can
## have unexpected effects downstream.  Consider examining the mean-variance plot
## and reducing the span parameter.

## conditions
## Retina_het  Retina_ko  Retina_wt    SCN_het     SCN_ko     SCN_wt   dLGN_het 
##          9          4          4          9          4          4          9 
##    dLGN_ko    dLGN_wt 
##          4          4

## conditions
## Retina_het  Retina_ko  Retina_wt    SCN_het     SCN_ko     SCN_wt   dLGN_het 
##          9          4          4          9          4          4          9 
##    dLGN_ko    dLGN_wt 
##          4          4

## conditions
## Retina_het  Retina_ko  Retina_wt    SCN_het     SCN_ko     SCN_wt   dLGN_het 
##          9          4          4          9          4          4          9 
##    dLGN_ko    dLGN_wt 
##          4          4

## Warning in correlate_de_tables(results, annot_df = annot_df, extra_contrasts =
## extra_contrasts): The deseq results are not of class hpgltools::deseq_pairwise.

location_table <- combine_de_tables(location_comparisons, excel = "excel/locations.xlsx")

## Looking for subscript invalid names, end of extract_keepers.

## Adding the coefficient plot to dLGN_ko_vs_SCN_ko.
## Adding the coefficient plot to dLGN_ko_vs_SCN_ko.
## Adding the coefficient plot to dLGN_ko_vs_SCN_ko.
## Adding the coefficient plot to dLGN_ko_vs_SCN_ko.
## Adding the coefficient plot to dLGN_ko_vs_SCN_ko.
## Adding the coefficient plot to dLGN_ko_vs_SCN_ko.

## Adding the coefficient plot to dLGN_ko_vs_Retina_ko.
## Adding the coefficient plot to dLGN_ko_vs_Retina_ko.
## Adding the coefficient plot to dLGN_ko_vs_Retina_ko.
## Adding the coefficient plot to dLGN_ko_vs_Retina_ko.
## Adding the coefficient plot to dLGN_ko_vs_Retina_ko.
## Adding the coefficient plot to dLGN_ko_vs_Retina_ko.

## Adding the coefficient plot to SCN_ko_vs_Retina_ko.
## Adding the coefficient plot to SCN_ko_vs_Retina_ko.
## Adding the coefficient plot to SCN_ko_vs_Retina_ko.
## Adding the coefficient plot to SCN_ko_vs_Retina_ko.
## Adding the coefficient plot to SCN_ko_vs_Retina_ko.
## Adding the coefficient plot to SCN_ko_vs_Retina_ko.

## Adding the coefficient plot to dLGN_het_vs_SCN_het.
## Adding the coefficient plot to dLGN_het_vs_SCN_het.
## Adding the coefficient plot to dLGN_het_vs_SCN_het.
## Adding the coefficient plot to dLGN_het_vs_SCN_het.
## Adding the coefficient plot to dLGN_het_vs_SCN_het.
## Adding the coefficient plot to dLGN_het_vs_SCN_het.

## Adding the coefficient plot to dLGN_het_vs_Retina_het.
## Adding the coefficient plot to dLGN_het_vs_Retina_het.
## Adding the coefficient plot to dLGN_het_vs_Retina_het.
## Adding the coefficient plot to dLGN_het_vs_Retina_het.
## Adding the coefficient plot to dLGN_het_vs_Retina_het.
## Adding the coefficient plot to dLGN_het_vs_Retina_het.

## Adding the coefficient plot to SCN_het_vs_Retina_het.
## Adding the coefficient plot to SCN_het_vs_Retina_het.
## Adding the coefficient plot to SCN_het_vs_Retina_het.
## Adding the coefficient plot to SCN_het_vs_Retina_het.
## Adding the coefficient plot to SCN_het_vs_Retina_het.
## Adding the coefficient plot to SCN_het_vs_Retina_het.

## Adding the coefficient plot to dLGN_wt_vs_SCN_wt.
## Adding the coefficient plot to dLGN_wt_vs_SCN_wt.
## Adding the coefficient plot to dLGN_wt_vs_SCN_wt.
## Adding the coefficient plot to dLGN_wt_vs_SCN_wt.
## Adding the coefficient plot to dLGN_wt_vs_SCN_wt.
## Adding the coefficient plot to dLGN_wt_vs_SCN_wt.

## Adding the coefficient plot to dLGN_wt_vs_Retina_wt.
## Adding the coefficient plot to dLGN_wt_vs_Retina_wt.
## Adding the coefficient plot to dLGN_wt_vs_Retina_wt.
## Adding the coefficient plot to dLGN_wt_vs_Retina_wt.
## Adding the coefficient plot to dLGN_wt_vs_Retina_wt.
## Adding the coefficient plot to dLGN_wt_vs_Retina_wt.

## Adding the coefficient plot to SCN_wt_vs_Retina_wt.
## Adding the coefficient plot to SCN_wt_vs_Retina_wt.
## Adding the coefficient plot to SCN_wt_vs_Retina_wt.
## Adding the coefficient plot to SCN_wt_vs_Retina_wt.
## Adding the coefficient plot to SCN_wt_vs_Retina_wt.
## Adding the coefficient plot to SCN_wt_vs_Retina_wt.

location_sig <- extract_significant_genes(location_table, according_to = "deseq",
                                             excel = "excel/location_sig.xlsx")

## Unable to discern the methods to test significance, according_to is NULL.

Save some of the DE results for comparison agains the TRAP data

de_list <- c("mm_de_dlgn_sva", "mm_de_ret_sva", "mm_de_scn_sva",
             "mm_de_dlgn", "mm_de_ret", "mm_de_scn",
             "genoloc_de_sva", "genoloc_table_sva")
saved <- save(list = de_list, file = "bulk_results.rda")