Section 4: Amnion Markers

setwd("~/github/BAPvsTrophoblast_Amnion")
library(sva)
library(tidyverse)
library(DESeq2)
library(vsn)
library(pheatmap)
library(ggrepel)
library(RColorBrewer)
library(reshape2)
library(scales)
require(biomaRt)
library(EnhancedVolcano)
library(TissueEnrich)
library(corrplot)
library(rstatix)
library(ggpubr)
library(car)
library(FSA)
library(PupillometryR)
load("assets/myCorrplot.Rdata")

# Roost and Suzuki data used as different amnion tissues
cts.full <-
  as.matrix(read.csv(
    "assets/counts-batch-corrected-v6_tpm.txt",
    sep = "\t",
    row.names = "gene"
  ))
md.full <-
  read.csv(
    "assets/coldata-batch-corrected-v6.txt",
    sep = "\t",
    header = TRUE,
    row.names = 1
  )
mycts.t <- as.data.frame(t(cts.full))
metadata <- md.full[2]
merged.table <-
  merge(metadata, mycts.t, by = "row.names", all.x = TRUE)
merged.table <-
  merged.table %>% remove_rownames %>% column_to_rownames(var = "Row.names")
mycts.mean <- aggregate(. ~ group, merged.table, mean)
ordered <-
  as.data.frame(t(
    mycts.mean %>% remove_rownames %>% column_to_rownames(var = "group")
  ))
se <-
  SummarizedExperiment(
    assays = SimpleList(as.matrix(ordered)),
    rowData = row.names(ordered),
    colData = colnames(ordered)
  )


# results for TE enrichment, max tissues =2; fc threshold = 10; and expression threshold 1
output <- teGeneRetrieval(
  se,
  foldChangeThreshold = 10,
  maxNumberOfTissues = 2,
  expressedGeneThreshold = 1
)

# save
te.results.merged <- as.data.frame(assay(output))
write.table(
  te.results.merged,
  file = "enrichment-results_fc10_merged.tsv",
  sep = "\t",
  row.names = FALSE,
  quote = FALSE
)

full.ec.results <- as.data.frame(assay(output))
amnion.roost <- subset(full.ec.results, Tissue=='amnion' & Group=='Tissue-Enriched')

ensembl = useMart("ENSEMBL_MART_ENSEMBL")
listDatasets(ensembl) %>%
  filter(str_detect(description, "Human"))
#>                 dataset              description    version
#> 1 hsapiens_gene_ensembl Human genes (GRCh38.p13) GRCh38.p13
ensembl = useDataset("hsapiens_gene_ensembl", mart = ensembl)
listFilters(ensembl) %>%
  filter(str_detect(name, "ensembl"))
#>                            name
#> 1               ensembl_gene_id
#> 2       ensembl_gene_id_version
#> 3         ensembl_transcript_id
#> 4 ensembl_transcript_id_version
#> 5            ensembl_peptide_id
#> 6    ensembl_peptide_id_version
#> 7               ensembl_exon_id
#>                                                      description
#> 1                       Gene stable ID(s) [e.g. ENSG00000000003]
#> 2       Gene stable ID(s) with version [e.g. ENSG00000000003.15]
#> 3                 Transcript stable ID(s) [e.g. ENST00000000233]
#> 4 Transcript stable ID(s) with version [e.g. ENST00000000233.10]
#> 5                    Protein stable ID(s) [e.g. ENSP00000000233]
#> 6     Protein stable ID(s) with version [e.g. ENSP00000000233.5]
#> 7                              Exon ID(s) [e.g. ENSE00000000003]
filterType <- "ensembl_gene_id_version"
filterValues <- rownames(cts)
attributeNames <- c('ensembl_gene_id_version',
                    'ensembl_gene_id',
                    'external_gene_name')
annot <- getBM(
  attributes = attributeNames,
  filters = filterType,
  values = filterValues,
  mart = ensembl
)
isDup <- duplicated(annot$ensembl_gene_id)
dup <- annot$ensembl_gene_id[isDup]
annot <- annot[!annot$ensembl_gene_id %in% dup, ]

addGeneIds <- function(df) {
  merge(
    df,
    annot,
    by.x = "Gene",
    by.y = "ensembl_gene_id",
    all.x = TRUE,
    all.y = FALSE
  ) %>% drop_na() %>% filter(external_gene_name != "")
}

cor.mtest <- function(mat, ...) {
  mat <- as.matrix(mat)
  n <- ncol(mat)
  p.mat<- matrix(NA, n, n)
  diag(p.mat) <- 0
  for (i in 1:(n - 1)) {
    for (j in (i + 1):n) {
      tmp <- cor.test(mat[, i], mat[, j], ...)
      p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
    }
  }
  colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
  p.mat
}

data_summary <- function(x) {
  m <- mean(x)
  ymin <- m-sd(x)
  ymax <- m+sd(x)
  return(c(y=m,ymin=ymin,ymax=ymax))
}

# SOURCE: https://ourcodingclub.github.io/tutorials/dataviz-beautification/
theme_niwot <- function() {
  theme_bw() +
    theme(
      axis.text = element_text(size = 16),
      axis.title = element_text(size = 18),
      axis.line.x = element_line(color = "black"),
      axis.line.y = element_line(color = "black"),
      panel.border = element_blank(),
      panel.grid.major.x = element_blank(),
      panel.grid.minor.x = element_blank(),
      panel.grid.minor.y = element_blank(),
      panel.grid.major.y = element_blank(),
      plot.margin = unit(c(1, 1, 1, 1), units = , "cm"),
      plot.title = element_text(
        size = 18,
        vjust = 1,
        hjust = 0
      ),
      legend.text = element_text(size = 12),
      legend.title = element_blank(),
      legend.position = c(0.95, 0.15),
      legend.key = element_blank(),
      legend.background = element_rect(
        color = "black",
        fill = "transparent",
        size = 2,
        linetype = "blank"
      )
    )
}

name1 <-
  c(
    "amnion_9w.1",
    "amnion_9w.2",
    "amnion_16w.2",
    "amnion_18w.1",
    "amnion_18w.2",
    "amnion_22w.1",
    "amnion_22w.2",
    "CT_7wk.1",
    "CT_7wk.2",
    "CT_9wk.1",
    "CT_11wk.1",
    "pBAP_D3.1",
    "pBAP_D3.2",
    "hESC_H9_BMP4_72h.1",
    "hESC_H9_BMP4_72h.2",
    "H1_gt70_D8_BAP.1",
    "H1_gt70_D8_BAP.2",
    "H1_gt70_D8_BAP.3"
  )
filter1 <- as.data.frame(name1)
name2 <-
  c(
    "BT_EVT_Okae.1",
    "BT_SCT_Okae.1",
    "BT_TSC_Okae.1",
    "BT_EVT_Okae.2",
    "BT_SCT_Okae.2",
    "BT_TSC_Okae.2",
    "CT_EVT_Okae.1",
    "CT_SCT_Okae.1",
    "CT_TSC_Okae.1",
    "CT_EVT_Okae.2",
    "CT_SCT_Okae.2",
    "CT_TSC_Okae.2",
    "CT_EVT_Okae.3",
    "CT_SCT_Okae.3",
    "CT_TSC_Okae.3",
    "amnion_18w.1",
    "amnion_9w.1",
    "amnion_18w.2",
    "amnion_16w.2",
    "amnion_22w.1",
    "amnion_9w.2",
    "amnion_22w.2"
  )
filter2 <- as.data.frame(name2)

getTable.main <- function(genelist = df,
                          data.to.plot = name1) {
  norm.filtered.table <-
    normalized_counts %>% filter(gene %in% genelist$ensembl_gene_id_version)
  norm.filtered.long <- melt(norm.filtered.table, id.vars = "gene")
  filter1 <- as.data.frame(data.to.plot)
  de_names <- filter1$data.to.plot
  norm.filtered.subset.long <-
    norm.filtered.long %>% filter (variable %in% de_names)
  colnames(norm.filtered.subset.long) <-
    c("gene", "condition", "norm.expression")
  norm.filtered.subset.long <-
    merge(
      norm.filtered.subset.long,
      annot,
      by.x = "gene",
      by.y = "ensembl_gene_id_version",
      all.x = TRUE,
      all.y = FALSE
    ) %>% drop_na() %>% filter(external_gene_name != "")
  norm.filtered.subset.long
}

runheatmap <-
  function(get_table = table,
           filename = "test_name",
           num = 20,
           data.to.plot = name1,
           ChartTitle = "Title") {
    norm.filtered.subset.long <- get_table
    norm.filtered.subset.table <-
      dcast(norm.filtered.subset.long,
            external_gene_name ~ condition,
            value.var = "norm.expression")
    write_delim(
      norm.filtered.subset.table,
      file = paste0(filename, "_expValues.tsv"),
      delim = "\t"
    )
    norm.filtered.subset.table <-
      norm.filtered.subset.table %>% column_to_rownames("external_gene_name")
    annotation <- data.frame(Condition = colnames(norm.filtered.subset.table))
    heat_colors <- brewer.pal(9, "YlOrRd")
    hmap.data <-
      as.matrix(norm.filtered.subset.table[, data.to.plot])
    g <- pheatmap(
      hmap.data,
      color = heat_colors,
      main = ChartTitle,
      cluster_rows = F,
      cluster_cols  = F,
      show_rownames = T,
      border_color = NA,
      fontsize = 14,
      scale = "row",
      fontsize_row = 10
    )
    g
  }

corPlot <- function(get_table = table,
                    ChartTitle = "Title",
                    mycolor = "darkorchid") {
  norm.filtered.subset.long <- get_table
  norm.filtered.subset.table <-
    dcast(norm.filtered.subset.long,
          external_gene_name ~ condition,
          value.var = "norm.expression")
  temp <- norm.filtered.subset.table[, -1]
  row.names(temp) <- norm.filtered.subset.table[, 1]
  norm.filtered.subset.table <- temp
  pcor <-
    cor(as.matrix(norm.filtered.subset.table), method = "spearman")
  p.mat <- cor.mtest(as.matrix(norm.filtered.subset.table))
  col <-
    colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))
  my.corrplot(
    pcor,
    method = "color",
    col = col(200),
    type = "upper",
    order = "hclust",
    addCoef.col = mycolor,
    pch.cex = 1.5,
    p.mat = p.mat,
    insig = 'label_sig',
    sig.level = c(0.001, 0.01, 0.05),
    tl.col = "black",
    tl.srt = 45,
    number.cex = 0.8,
    tl.cex = 1,
    pch.col = "tomato",
    diag = FALSE,
    font.main = 4,
    mar = c(0, 0, 1, 0)
  )
}

normTestEach <- function(get_table = table,
                         ChartTitle = "Title") {
  norm.filtered.subset.long <- get_table
  ggplot(norm.filtered.subset.long, aes(x = norm.expression)) +
    stat_density(color = "darkblue", fill = "lightblue") +
    theme_minimal() +
    theme(
      axis.text.x = element_text(
        angle = 45,
        vjust = 1,
        hjust = 1,
        size = 12
      ),
      axis.text.y = element_text(size = 12),
      plot.margin = margin(10, 10, 10, 100),
      legend.position = "none",
      plot.title = element_text(
        color = "black",
        size = 18,
        face = "bold.italic"
      ),
      axis.line.x = element_line(
        colour = 'black',
        size = 0.5,
        linetype = 'solid'
      ),
      axis.ticks.x = element_line(
        colour = 'black',
        size = 1,
        linetype = 'solid'
      ),
      axis.title.x = element_text(
        color = "black",
        size = 14,
        face = "bold"
      ),
      axis.line.y = element_line(
        colour = 'black',
        size = 0.5,
        linetype = 'solid'
      ),
      axis.ticks.y = element_line(
        colour = 'black',
        size = 1,
        linetype = 'solid'
      ),
      axis.title.y = element_text(
        color = "black",
        size = 14,
        face = "bold"
      )
    )  +
    scale_y_continuous(expand = expansion(mult = c(0, .1)), breaks = pretty_breaks()) +
    facet_wrap( ~ condition, ncol = 5) +
    scale_x_log10()
}

normTestAll <- function(get_table = table,
                        ChartTitle = "Title") {
  norm.filtered.subset.long <- get_table
  ggplot(norm.filtered.subset.long, aes(x = norm.expression)) +
    stat_density(color = "purple", fill = "dodgerblue") +
    theme_minimal() +
    theme(
      axis.text.x = element_text(
        angle = 45,
        vjust = 1,
        hjust = 1,
        size = 12
      ),
      axis.text.y = element_text(size = 12),
      plot.margin = margin(10, 10, 10, 100),
      legend.position = "none",
      plot.title = element_text(
        color = "black",
        size = 18,
        face = "bold.italic"
      ),
      axis.line.x = element_line(
        colour = 'black',
        size = 0.5,
        linetype = 'solid'
      ),
      axis.line.y = element_line(
        colour = 'black',
        size = 0.5,
        linetype = 'solid'
      ),
      axis.ticks.x = element_line(
        colour = 'black',
        size = 1,
        linetype = 'solid'
      ),
      axis.ticks.y = element_line(
        colour = 'black',
        size = 1,
        linetype = 'solid'
      ),
      axis.title.x = element_text(
        color = "black",
        size = 14,
        face = "bold"
      ),
      axis.title.y = element_text(
        color = "black",
        size = 14,
        face = "bold"
      )
    )  +
    scale_y_continuous(expand = expansion(mult = c(0, .1)), breaks = pretty_breaks()) +
    scale_x_log10()
}

addSampleInfo <- function(get_table = table) {
  norm.filtered.subset.long <- get_table
  metadata <- coldata %>%
    rownames_to_column(var = "condition") %>%
    select(condition, group)
  metadata$group <- sub("^(amnion_).*w", "Amnion", metadata$group)
  metadata$group <-
    sub("^(CT_7|CT_9|CT_11)wk", "CT", metadata$group)
  norm.filtered.subset.long <- merge(
    norm.filtered.subset.long,
    metadata,
    by.x = "condition",
    by.y = "condition",
    all.x = TRUE,
    all.y = FALSE
  )
  norm.filtered.subset.long
}

vlnPlot <- function(info_table = Infotable) {
  norm.filtered.subset.long <- info_table
  ggplot(data = norm.filtered.subset.long,
         aes(x = group, y = norm.expression, fill = group)) +
    geom_flat_violin(position = position_nudge(x = 0.2, y = 0),
                     alpha = 0.8,
                     trim = FALSE) +
    geom_point(
      aes(y = norm.expression, color = group),
      position = position_jitter(width = 0.15),
      size = 1,
      alpha = 0.5
    ) +
    geom_boxplot(width = 0.2,
                 outlier.shape = NA,
                 alpha = 0.8) + stat_summary(
                   fun = mean,
                   geom = "point",
                   shape = 23,
                   size = 2
                 ) +
    labs(y = "\nNormalized Expression", x = NULL) +
    guides(fill = "none", color = "none") +
    scale_y_log10()  +
    theme_niwot() + theme(axis.text.x = element_text(
      angle = 45,
      vjust = 1,
      hjust = 1
    ))
}

amnion.roost.names <- addGeneIds(amnion.roost)

file1 <-
  "~/TutejaLab/amnion_for_ms_20210715/heatmaps.v2/Fig6b_and_monkeygenes/genes-Fig6b.txt"
fig6b <-
  read.csv(file1,
           sep = "\t",
           header = TRUE,
           stringsAsFactors = FALSE)
fig6b <- merge(
  fig6b,
  annot,
  by.x = "genes",
  by.y = "external_gene_name",
  all.x = TRUE,
  all.y = FALSE
) %>% drop_na() 

amnion.roost.new.main <- getTable.main(amnion.roost.names, name1)
amnion.roost.old.main <- getTable.main(fig6b, name1)
amnion.roost.new.supp <- getTable.main(amnion.roost.names, name2)
amnion.roost.old.supp <- getTable.main(fig6b, name2)

amnion.roost.new.main.info <- addSampleInfo(amnion.roost.new.main)
amnion.roost.old.main.info <- addSampleInfo(amnion.roost.old.main)
amnion.roost.new.supp.info <- addSampleInfo(amnion.roost.new.supp)
amnion.roost.old.supp.info <- addSampleInfo(amnion.roost.old.supp)

runheatmap(
  amnion.roost.new.main.info,
  data.to.plot = name1,
  num = 38,
  ChartTitle = "New amnion markers - expression across samples"
)

corPlot(amnion.roost.new.main.info)

normTestEach(amnion.roost.new.main.info)
#> Warning: Transformation introduced infinite values in continuous x-axis
#> Warning: Removed 3046 rows containing non-finite values (stat_density).

normTestAll(amnion.roost.new.main.info)
#> Warning: Transformation introduced infinite values in continuous x-axis
#> Warning: Removed 3046 rows containing non-finite values (stat_density).

shapiro.each <-
  amnion.roost.new.main.info %>% dplyr::group_by(condition) %>% shapiro_test(norm.expression)
DT::datatable(shapiro.each,)

shapiro.test(amnion.roost.new.main.info$norm.expression[0:5000])
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  amnion.roost.new.main.info$norm.expression[0:5000]
#> W = 0.06815, p-value < 2.2e-16

kruskal.test(norm.expression ~ group, data = amnion.roost.new.main.info)
#> 
#>  Kruskal-Wallis rank sum test
#> 
#> data:  norm.expression by group
#> Kruskal-Wallis chi-squared = 414.88, df = 4, p-value < 2.2e-16

dunnTest <-
  dunnTest(norm.expression ~ group, data = amnion.roost.new.main.info, method =
             "bh")
#> Warning: group was coerced to a factor.
g <- as.data.frame(dunnTest$res)
DT::datatable(g)

vlnPlot(amnion.roost.new.main.info)
#> Warning: Transformation introduced infinite values in continuous y-axis

#> Warning: Transformation introduced infinite values in continuous y-axis

#> Warning: Transformation introduced infinite values in continuous y-axis

#> Warning: Transformation introduced infinite values in continuous y-axis
#> Warning: Removed 3046 rows containing non-finite values (stat_ydensity).
#> Warning: Removed 3046 rows containing non-finite values (stat_boxplot).
#> Warning: Removed 3046 rows containing non-finite values (stat_summary).
#> Warning: Removed 3046 rows containing missing values (geom_point).

runheatmap(
  amnion.roost.old.main.info,
  data.to.plot = name1,
  num = 20,
  ChartTitle = "Io amnion markers - expression across samples"
)

corPlot(amnion.roost.old.main.info, mycolor = "yellow")

normTestEach(amnion.roost.old.main.info)
#> Warning: Transformation introduced infinite values in continuous x-axis
#> Warning: Removed 230 rows containing non-finite values (stat_density).

normTestAll(amnion.roost.old.main.info)
#> Warning: Transformation introduced infinite values in continuous x-axis
#> Warning: Removed 230 rows containing non-finite values (stat_density).

shapiro.each <-
  amnion.roost.old.main.info %>% dplyr::group_by(condition) %>% shapiro_test(norm.expression)
DT::datatable(shapiro.each,)

shapiro.test(amnion.roost.old.main.info$norm.expression[0:5000])
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  amnion.roost.old.main.info$norm.expression[0:5000]
#> W = 0.24651, p-value < 2.2e-16

kruskal.test(norm.expression ~ group, data = amnion.roost.old.main.info)
#> 
#>  Kruskal-Wallis rank sum test
#> 
#> data:  norm.expression by group
#> Kruskal-Wallis chi-squared = 81.725, df = 4, p-value < 2.2e-16

dunnTest <-
  dunnTest(norm.expression ~ group, data = amnion.roost.old.main.info, method =
             "bh")
#> Warning: group was coerced to a factor.
g <- as.data.frame(dunnTest$res)
DT::datatable(g)

vlnPlot(amnion.roost.old.main.info)
#> Warning: Transformation introduced infinite values in continuous y-axis

#> Warning: Transformation introduced infinite values in continuous y-axis

#> Warning: Transformation introduced infinite values in continuous y-axis

#> Warning: Transformation introduced infinite values in continuous y-axis
#> Warning: Removed 230 rows containing non-finite values (stat_ydensity).
#> Warning: Removed 230 rows containing non-finite values (stat_boxplot).
#> Warning: Removed 230 rows containing non-finite values (stat_summary).
#> Warning: Removed 230 rows containing missing values (geom_point).

runheatmap(
  amnion.roost.new.supp.info,
  data.to.plot = name2,
  num = 20,
  ChartTitle = "Io amnion markers - expression across samples"
)

corPlot(amnion.roost.new.supp.info)

normTestEach(amnion.roost.new.supp.info)
#> Warning: Transformation introduced infinite values in continuous x-axis
#> Warning: Removed 4396 rows containing non-finite values (stat_density).

normTestAll(amnion.roost.new.supp.info)
#> Warning: Transformation introduced infinite values in continuous x-axis
#> Warning: Removed 4396 rows containing non-finite values (stat_density).

shapiro.each <-
  amnion.roost.new.supp.info %>% dplyr::group_by(condition) %>% shapiro_test(norm.expression)
DT::datatable(shapiro.each,)

shapiro.test(amnion.roost.new.supp.info$norm.expression[0:5000])
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  amnion.roost.new.supp.info$norm.expression[0:5000]
#> W = 0.065747, p-value < 2.2e-16

kruskal.test(norm.expression ~ group, data = amnion.roost.new.supp.info)
#> 
#>  Kruskal-Wallis rank sum test
#> 
#> data:  norm.expression by group
#> Kruskal-Wallis chi-squared = 901.67, df = 6, p-value < 2.2e-16

dunnTest <-
  dunnTest(norm.expression ~ group, data = amnion.roost.new.supp.info, method =
             "bh")
#> Warning: group was coerced to a factor.
g <- as.data.frame(dunnTest$res)
DT::datatable(g)

vlnPlot(amnion.roost.new.supp.info)
#> Warning: Transformation introduced infinite values in continuous y-axis

#> Warning: Transformation introduced infinite values in continuous y-axis

#> Warning: Transformation introduced infinite values in continuous y-axis

#> Warning: Transformation introduced infinite values in continuous y-axis
#> Warning: Removed 4396 rows containing non-finite values (stat_ydensity).
#> Warning: Removed 4396 rows containing non-finite values (stat_boxplot).
#> Warning: Removed 4396 rows containing non-finite values (stat_summary).
#> Warning: Removed 4396 rows containing missing values (geom_point).

runheatmap(
  amnion.roost.old.supp.info,
  data.to.plot = name2,
  num = 20,
  ChartTitle = "Io amnion markers - expression across samples"
)

corPlot(amnion.roost.old.supp.info, mycolor = "yellow")

normTestEach(amnion.roost.old.supp.info)
#> Warning: Transformation introduced infinite values in continuous x-axis
#> Warning: Removed 325 rows containing non-finite values (stat_density).

normTestAll(amnion.roost.old.supp.info)
#> Warning: Transformation introduced infinite values in continuous x-axis
#> Warning: Removed 325 rows containing non-finite values (stat_density).

shapiro.each <-
  amnion.roost.old.supp.info %>% dplyr::group_by(condition) %>% shapiro_test(norm.expression)
DT::datatable(shapiro.each,)

shapiro.test(amnion.roost.old.supp.info$norm.expression[0:5000])
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  amnion.roost.old.supp.info$norm.expression[0:5000]
#> W = 0.13351, p-value < 2.2e-16

kruskal.test(norm.expression ~ group, data = amnion.roost.old.supp.info)
#> 
#>  Kruskal-Wallis rank sum test
#> 
#> data:  norm.expression by group
#> Kruskal-Wallis chi-squared = 7.1575, df = 6, p-value = 0.3065

dunnTest <-
  dunnTest(norm.expression ~ group, data = amnion.roost.old.supp.info, method =
             "bh")
#> Warning: group was coerced to a factor.
g <- as.data.frame(dunnTest$res)
DT::datatable(g)

vlnPlot(amnion.roost.old.supp.info)
#> Warning: Transformation introduced infinite values in continuous y-axis

#> Warning: Transformation introduced infinite values in continuous y-axis

#> Warning: Transformation introduced infinite values in continuous y-axis

#> Warning: Transformation introduced infinite values in continuous y-axis
#> Warning: Removed 325 rows containing non-finite values (stat_ydensity).
#> Warning: Removed 325 rows containing non-finite values (stat_boxplot).
#> Warning: Removed 325 rows containing non-finite values (stat_summary).
#> Warning: Removed 325 rows containing missing values (geom_point).

sessionInfo()
#> R version 4.0.5 (2021-03-31)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 19044)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.1252 
#> [2] LC_CTYPE=English_United States.1252   
#> [3] LC_MONETARY=English_United States.1252
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.1252    
#> 
#> attached base packages:
#> [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
#> [8] methods   base     
#> 
#> other attached packages:
#>  [1] PupillometryR_0.0.4         rlang_0.4.11               
#>  [3] FSA_0.9.1                   car_3.0-11                 
#>  [5] carData_3.0-4               ggpubr_0.4.0               
#>  [7] rstatix_0.7.0               corrplot_0.90              
#>  [9] TissueEnrich_1.10.1         GSEABase_1.52.1            
#> [11] graph_1.68.0                annotate_1.68.0            
#> [13] XML_3.99-0.6                AnnotationDbi_1.52.0       
#> [15] ensurer_1.1                 EnhancedVolcano_1.8.0      
#> [17] biomaRt_2.46.3              scales_1.1.1               
#> [19] reshape2_1.4.4              RColorBrewer_1.1-2         
#> [21] ggrepel_0.9.1               pheatmap_1.0.12            
#> [23] vsn_3.58.0                  DESeq2_1.30.1              
#> [25] SummarizedExperiment_1.20.0 Biobase_2.50.0             
#> [27] MatrixGenerics_1.2.1        matrixStats_0.58.0         
#> [29] GenomicRanges_1.42.0        GenomeInfoDb_1.26.7        
#> [31] IRanges_2.24.1              S4Vectors_0.28.1           
#> [33] BiocGenerics_0.36.1         forcats_0.5.1              
#> [35] stringr_1.4.0               dplyr_1.0.7                
#> [37] purrr_0.3.4                 readr_2.1.1                
#> [39] tidyr_1.1.4                 tibble_3.1.1               
#> [41] ggplot2_3.3.5               tidyverse_1.3.1            
#> [43] sva_3.38.0                  BiocParallel_1.24.1        
#> [45] genefilter_1.72.1           mgcv_1.8-35                
#> [47] nlme_3.1-152               
#> 
#> loaded via a namespace (and not attached):
#>  [1] readxl_1.3.1          backports_1.2.1       BiocFileCache_1.14.0 
#>  [4] plyr_1.8.6            splines_4.0.5         crosstalk_1.1.1      
#>  [7] digest_0.6.27         htmltools_0.5.2       fansi_0.4.2          
#> [10] magrittr_2.0.1        memoise_2.0.1         openxlsx_4.2.4       
#> [13] tzdb_0.2.0            limma_3.46.0          modelr_0.1.8         
#> [16] extrafont_0.17        vroom_1.5.7           extrafontdb_1.0      
#> [19] askpass_1.1           rmdformats_1.0.3      prettyunits_1.1.1    
#> [22] colorspace_2.0-1      blob_1.2.2            rvest_1.0.0          
#> [25] rappdirs_0.3.3        haven_2.4.1           xfun_0.29            
#> [28] crayon_1.4.2          RCurl_1.98-1.3        jsonlite_1.7.3       
#> [31] survival_3.2-11       glue_1.4.2            gtable_0.3.0         
#> [34] zlibbioc_1.36.0       XVector_0.30.0        DelayedArray_0.16.3  
#> [37] proj4_1.0-10.1        dunn.test_1.3.5       Rttf2pt1_1.3.8       
#> [40] maps_3.3.0            abind_1.4-5           DBI_1.1.2            
#> [43] edgeR_3.32.1          Rcpp_1.0.8            xtable_1.8-4         
#> [46] progress_1.2.2        foreign_0.8-81        bit_4.0.4            
#> [49] preprocessCore_1.52.1 DT_0.18               htmlwidgets_1.5.3    
#> [52] httr_1.4.2            ellipsis_0.3.2        farver_2.1.0         
#> [55] pkgconfig_2.0.3       sass_0.4.0            dbplyr_2.1.1         
#> [58] locfit_1.5-9.4        utf8_1.2.1            tidyselect_1.1.1     
#> [61] munsell_0.5.0         cellranger_1.1.0      tools_4.0.5          
#> [64] cachem_1.0.5          cli_3.1.1             generics_0.1.1       
#> [67] RSQLite_2.2.9         broom_0.7.6           evaluate_0.14        
#> [70] fastmap_1.1.0         yaml_2.2.1            knitr_1.37           
#> [73] bit64_4.0.5           fs_1.5.0              zip_2.1.1            
#>  [ reached getOption("max.print") -- omitted 40 entries ]