2_de-analyses_v2.Rmd

---
title: "Section 2: DESeq2 analysis"
date: "`r Sys.Date()`"
author:
  - name: Arun Seetharam
  - name: Ha Vu
    affiliation: Tuteja Lab
    affiliation_url: https://www.tutejalab.org
output:
  rmdformats::readthedown:
    self_contained: true
    thumbnails: false
    lightbox: true
    gallery: true
    highlight: tango
---

```{r setup, include=FALSE}
options(max.print = "75")
knitr::opts_chunk$set(
  echo = TRUE,
  collapse = TRUE,
  comment = "#>",
  fig.path = "assets/",
  fig.width = 8,
  prompt = FALSE,
  tidy = FALSE,
  message = FALSE,
  warning = TRUE
)
knitr::opts_knit$set(width = 75)
```

## DESeq2 analyses steps

This section uses the count data (for selected datasets) generated in Section 1 to do differential expression (DE) analyses using `DESeq2`. Briefly, the count data are imported in R, batch corrected using `ComBat_seq`, then DE analyses were performed for various contrasts using `DESeq2`. Results were visualized as volcano plots, and cell enrichment performed using `PlacentaCellEnrich` (`PCE`).

## Prerequisites

R packages required for this section are loaded

```{r, warnings=TRUE, message=FALSE}
setwd("~/github/BAPvsTrophoblast_Amnion")
# load the modules
library(sva)
library(tidyverse)
library(DESeq2)
library(vsn)
library(pheatmap)
library(ggrepel)
library(RColorBrewer)
library(reshape2)
require(biomaRt)
library(EnhancedVolcano)
library(TissueEnrich)
library(plotly)
library(DT)
library(cowplot)
library(biomaRt)
library(tidytext)
library(ggpubr)
library(scales)
```


## Import datasets

The `counts` data and its associated metadata (`coldata`) are imported for analyses.

```{r dataset, warnings=TRUE, message=FALSE}
counts = 'assets/counts-subset-v5.txt'
groupFile = 'assets/batch-subset-v5.txt'
coldata <-
  read.csv(
    groupFile,
    row.names = 1,
    sep = "\t",
    stringsAsFactors = TRUE
  )
cts <- as.matrix(read.csv(counts, sep = "\t", row.names = "gene.ids"))
```

Inspect the `coldata`.

```{r coldata}
DT::datatable(coldata)
```

Reorder columns of `cts` according to `coldata` rows. Check if samples in both files match.

```{r order, warnings=TRUE, message=FALSE}
colnames(cts)
all(rownames(coldata) %in% colnames(cts))
cts <- cts[, rownames(coldata)]
```

## Batch correction

Using ComBat_seq (SVA package) to run batch correction - using bioproject IDs as variable (dataset origin).

```{r batchcorrect, warnings=TRUE, message=FALSE}
cov1 <- as.factor(coldata$BioProject)
adjusted_counts <- ComBat_seq(cts, batch = cov1, group = NULL)
all(rownames(coldata) %in% colnames(cts))
cts <- cts[, rownames(coldata)]
```

## DESeq2

The batch corrected read counts are then used for running DESeq2 analyses

```{r deseq2, warnings=TRUE, message=FALSE}
dds <- DESeqDataSetFromMatrix(countData = adjusted_counts,
                              colData = coldata,
                              design = ~ condition)
vsd <- vst(dds, blind = FALSE)
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep, ]
dds <- DESeq(dds)
dds
```

Various contrasts are set up as follows (a total of 8 combinations)

```{r cotnrasts, warnings=TRUE, message=FALSE}
res.PH9vsBAP <-
  results(dds,
          contrast = c(
            "condition", 
            "Primed_H9_hESCs", 
            "pBAP_D3.Primed_H9_hESCs"))

res.K00vsK72 <-
  results(dds,
          contrast = c(
            "condition", 
            "hESC_H9_untr_0h",
            "hESC_H9_BMP4_72h"))

res.UNDvsSTB <-
  results(dds,
          contrast = c(
            "condition",
            "hESC_H1_D8_MEF.CM.and.FGF2",
            "hESC_H1_STB_gt70um_D8_BAP"))


res.K72vsSTB <-
  results(dds,
          contrast = c(
            "condition", 
            "hESC_H9_BMP4_72h",
            "hESC_H1_STB_gt70um_D8_BAP"))

res.BAPvsSTB <-
  results(
    dds,
    contrast = c(
      "condition",
      "pBAP_D3.Primed_H9_hESCs",
      "hESC_H1_STB_gt70um_D8_BAP"))

res.BAPvsK72 <-
  results(dds,
          contrast = c(
            "condition", 
            "pBAP_D3.Primed_H9_hESCs",
            "hESC_H9_BMP4_72h"))


res.K72vsL40 <-
  results(dds,
          contrast = c(
            "condition", 
            "hESC_H9_BMP4_72h",
            "hESC_H1_STB_lt40um_D8_BAP"))

res.BAPvsL40 <-
  results(
    dds,
    contrast = c(
      "condition",
      "pBAP_D3.Primed_H9_hESCs",
      "hESC_H1_STB_lt40um_D8_BAP"
    )
  )
```

The following function is to save DESeq2 results as well as generate variables to hold the gene lists for running PCE later on.

```{r deseq2funx, warnings=TRUE, message=FALSE}
processDE <- function(res.se, string) {
  res.se <- res.se[order(res.se$padj),]
  res.data <-
    merge(as.data.frame(res.se),
          as.data.frame(counts(dds, normalized = TRUE)),
          by = "row.names",
          sort = FALSE)
  names(res.data)[1] <- "Gene"
  write_delim(res.data,
              file = paste0("DESeq2results-", string, "_fc.tsv"),
              delim = "\t")
  res.up <-
    res.data %>% 
    filter(log2FoldChange >= 1) %>% 
    filter(padj <= 0.05) %>% 
    arrange(desc(log2FoldChange)) %>% 
    dplyr::select(Gene)
  res.dw <-
    res.data %>% 
    filter(log2FoldChange <= -1) %>% 
    filter(padj <= 0.05) %>% 
    arrange(desc(log2FoldChange)) %>% 
    dplyr::select(Gene)
  res.up.new <-
    annot[annot$ensembl_gene_id_version %in% res.up$Gene, ]
  res.dw.new <-
    annot[annot$ensembl_gene_id_version %in% res.dw$Gene, ]
  pce.up1 <- paste0(string, ".up.pce", 1)
  pce.dw1 <- paste0(string, ".dw.pce", 1)
  pce.up2 <- paste0(string, ".up.pce", 2)
  pce.dw2 <- paste0(string, ".dw.pce", 2)
  assign(pce.up1, res.up.new$ensembl_gene_id, envir = .GlobalEnv)
  assign(pce.dw1, res.dw.new$ensembl_gene_id, envir = .GlobalEnv)
  assign(pce.up2, res.up.new$external_gene_name, envir = .GlobalEnv)
  assign(pce.dw2, res.dw.new$external_gene_name, envir = .GlobalEnv)
}
```

## Creating gene lists

The gene lists have Ensembl gene-ID-version. We need them as gene-IDs. We also need other metadata later for these lists.
From Ensembl we will download metadata and attach to these lists.

```{r annotations, warnings=TRUE, message=FALSE}
ensembl = useMart("ENSEMBL_MART_ENSEMBL")
listDatasets(ensembl) %>%
  filter(str_detect(description, "Human"))
ensembl = useDataset("hsapiens_gene_ensembl", mart = ensembl)
listFilters(ensembl) %>%
  filter(str_detect(name, "ensembl"))
filterType <- "ensembl_gene_id_version"
filterValues <- rownames(cts)
listAttributes(ensembl) %>%
  head(20)
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, ] #this object will be saved and used later
```

The results are saved as tsv files.

```{r deseq2save, warnings=TRUE, message=FALSE}
processDE(res.PH9vsBAP, "PH9vsBAP")
processDE(res.K00vsK72, "K00vsK72")
processDE(res.UNDvsSTB, "UNDvsSTB")
processDE(res.K72vsSTB, "K72vsSTB")
processDE(res.BAPvsSTB, "BAPvsSTB")
processDE(res.BAPvsK72, "BAPvsK72")
processDE(res.K72vsL40, "K72vsL40")
processDE(res.BAPvsL40, "BAPvsL40")
```


```{r martdata}
mart <-
  read.csv(
    "assets/mart-genes.tsv",
    sep = "\t",
    stringsAsFactors = TRUE,
    header = TRUE
  ) #this object was obtained from Ensembl as we illustrated in "Creating gene lists"
```

```{r volcano}
volcanoPlots <-
  function(res.se,
           string,
           first,
           second,
           color1,
           color2,
           color3,
           ChartTitle) {
    res.se <- res.se[order(res.se$padj), ]
    res.se <-
      rownames_to_column(as.data.frame(res.se[order(res.se$padj), ]))
    names(res.se)[1] <- "Gene"
    res.data <-
      merge(res.se,
            mart,
            by.x = "Gene",
            by.y = "ensembl_gene_id_version")
    res.data <- res.data %>% mutate_all(na_if, "")
    res.data <- res.data %>% mutate_all(na_if, " ")
    res.data <-
      res.data %>% mutate(gene_symbol = coalesce(gene_symbol, Gene))
    res.data$diffexpressed <- "other.genes"
    res.data$diffexpressed[res.data$log2FoldChange >= 1 &
                             res.data$padj <= 0.05] <-
      paste("Higher expression in", first)
    res.data$diffexpressed[res.data$log2FoldChange <= -1 &
                             res.data$padj <= 0.05] <-
      paste("Higher expression in", second)
    res.data$delabel <- ""
    res.data$delabel[res.data$log2FoldChange >= 1
                     & res.data$padj <= 0.05
                     &
                       !is.na(res.data$padj)] <-
      res.data$gene_symbol[res.data$log2FoldChange >= 1
                           &
                             res.data$padj <= 0.05
                           &
                             !is.na(res.data$padj)]
    res.data$delabel[res.data$log2FoldChange <= -1
                     & res.data$padj <= 0.05
                     &
                       !is.na(res.data$padj)] <-
      res.data$gene_symbol[res.data$log2FoldChange <= -1
                           &
                             res.data$padj <= 0.05
                           &
                             !is.na(res.data$padj)]
    outpath <- "interactive/"
    gg <-
      ggplot(res.data,
             aes(
               x = log2FoldChange,
               y = -log10(padj),
               col = diffexpressed,
               label = delabel
             )) +
      geom_point(alpha = 0.5) +
      xlim(-20, 20) +
      theme_classic() +
      scale_color_manual(name = "Expression", values = c(color1, color2, color3)) +
      ggtitle(ChartTitle) +
      xlab(paste("log2 fold change")) +
      ylab("-log10 pvalue (adjusted)") +
      theme(legend.text.align = 0)
    saveWidget(ggplotly(gg), file = paste0(outpath, "/Figure_volcano_", string, ".html"))
}
```

## Volcano Plots (interactive)

Running Volcano plots for each comparison are shown below.

```{r vol_interactive_all, warnings=TRUE, message=FALSE}
volcanoPlots(
  res.PH9vsBAP,
  "PH9vsBAP",
  "pH9_Io",
  "H9_pBAP_D3_Io",
  "#0571B0",
  "#483D8B",
  "#4d4d4d",
  ChartTitle = "pH9_Io vs. H9_pBAP_D3_Io"
)
volcanoPlots(
  res.K00vsK72,
  "K00vsK72",
  "H9_BMP4.0h_Krendl",
  "H9_BMP4.72h_Krendl",
  "#FF1493",
  "#EE82EE",
  "#4d4d4d",
  ChartTitle = "H9_BMP4.0h_Krendl vs. H9_BMP4.72h_Krendl"
)
volcanoPlots(
  res.UNDvsSTB,
  "UNDvsSTB",
  "H1_Yabe",
  "H1_BAP_D8_>70_Yabe",
  "#598234",
  "#006400",
  "#4d4d4d",
  ChartTitle = "H1_Yabe vs. H1_BAP_D8_>70_Yabe"
)
volcanoPlots(
  res.K72vsSTB,
  "K72vsSTB",
  "H9_BMP4.72h_Krendl",
  "H1_BAP_D8_>70_Yabe",
  "#598234",
  "#EE82EE",
  "#4d4d4d",
  ChartTitle = "H9_BMP4.72h_Krendl vs. H1_BAP_D8_>70_Yabe"
)
volcanoPlots(
  res.BAPvsSTB,
  "BAPvsSTB",
  "H9_pBAP_D3_Io",
  "H1_BAP_D8_>70_Yabe",
  "#598234",
  "#0571B0",
  "#4d4d4d",
  ChartTitle = "H9_pBAP_D3_Io vs. H1_BAP_D8_>70_Yabe"
)
volcanoPlots(
  res.BAPvsK72,
  "BAPvsK72",
  "H9_pBAP_D3_Io",
  "H9_BMP4.72h_Krendl",
  "#EE82EE",
  "#0571B0",
  "#4d4d4d",
  ChartTitle = "H9_pBAP_D3_Io vs. H9_BMP4.72h_Krendl"
)
volcanoPlots(
  res.K72vsL40,
  "K72vsL40",
  "H9_BMP4.72h_Krendl",
  "H1_BAP_D8_<40_Yabe",
  "#AEBD38",
  "#EE82EE",
  "#4d4d4d",
  ChartTitle = "H9_BMP4.72h_Krendl vs. H1_BAP_D8_<40_Yabe"
)
volcanoPlots(
  res.BAPvsL40,
  "BAPvsL40",
  "H9_pBAP_D3_Io",
  "H1_BAP_D8_<40_Yabe",
  "#AEBD38",
  "#0571B0",
  "#4d4d4d",
  ChartTitle = "H9_pBAP_D3_Io vs. H1_BAP_D8_<40_Yabe"
)
```
### Interactive Volcano Plots:

1. [pH9_Io vs. H9_pBAP_D3_Io](interactive/Figure_volcano_ PH9vsBAP.html){target="_blank"}
2. [H9_BMP4.0h_Krendl vs. H9_BMP4.72h_Krendl](interactive/Figure_volcano_ K00vsK72.html){target="_blank"}
3. [H1_Yabe vs. H1_BAP_D8_>70_Yabe](interactive/Figure_volcano_ UNDvsSTB.html){target="_blank"}
4. [H9_BMP4.72h_Krendl vs. H1_BAP_D8_>70_Yabe](interactive/Figure_volcano_K72vsSTB.html){target="_blank"}
5. [H9_pBAP_D3_Io vs. H1_BAP_D8_>70_Yabe](interactive/Figure_volcano_ BAPvsSTB.html){target="_blank"}
6. [H9_pBAP_D3_Io vs. H9_BMP4.72h_Krendl](interactive/Figure_volcano_ BAPvsK72.html){target="_blank"}
7. [H9_BMP4.72h_Krendl vs. H1_BAP_D8_<40_Yabe](interactive/Figure_volcano_K72vsL40.html){target="_blank"}
8. [H9_pBAP_D3_Io vs. H1_BAP_D8_<40_Yabe](interactive/Figure_volcano_ BAPvsL40.html){target="_blank"}

### Static Volcano Plots:

```{r volcano2}
volcanoPlots2 <-
  function(res.se,
           string,
           first,
           second,
           color1,
           color2,
           color3,
           ChartTitle) {
    res.se <- res.se[order(res.se$padj), ]
    res.se <-
      rownames_to_column(as.data.frame(res.se[order(res.se$padj), ]))
    names(res.se)[1] <- "Gene"
    res.data <-
      merge(res.se,
            mart,
            by.x = "Gene",
            by.y = "ensembl_gene_id_version")
    res.data <- res.data %>% mutate_all(na_if, "")
    res.data <- res.data %>% mutate_all(na_if, " ")
    res.data <-
      res.data %>% mutate(gene_symbol = coalesce(gene_symbol, Gene))
    res.data$diffexpressed <- "other.genes"
    res.data$diffexpressed[res.data$log2FoldChange >= 1 &
                             res.data$padj <= 0.05] <-
      paste("Higher expression in", first)
    res.data$diffexpressed[res.data$log2FoldChange <= -1 &
                             res.data$padj <= 0.05] <-
      paste("Higher expression in", second)
    res.data$delabel <- ""
    res.data$delabel[res.data$log2FoldChange >= 1
                     & res.data$padj <= 0.05
                     &
                       !is.na(res.data$padj)] <-
      res.data$gene_symbol[res.data$log2FoldChange >= 1
                           &
                             res.data$padj <= 0.05
                           &
                             !is.na(res.data$padj)]
    res.data$delabel[res.data$log2FoldChange <= -1
                     & res.data$padj <= 0.05
                     &
                       !is.na(res.data$padj)] <-
      res.data$gene_symbol[res.data$log2FoldChange <= -1
                           &
                             res.data$padj <= 0.05
                           &
                             !is.na(res.data$padj)]
    ggplot(res.data,
             aes(
               x = log2FoldChange,
               y = -log10(padj),
               col = diffexpressed,
               label = delabel
             )) +
      geom_point(alpha = 0.5) +
      xlim(-20, 20) +
      theme_classic() +
      scale_color_manual(name = "Expression", values = c(color1, color2, color3)) +
      geom_text_repel(
        data = subset(res.data, padj <= 0.05),
        max.overlaps  = 15,
        show.legend = F,
        min.segment.length = Inf,
        seed = 42,
        box.padding = 0.5
      ) +
      ggtitle(ChartTitle) +
      xlab(paste("log2 fold change")) +
      ylab("-log10 pvalue (adjusted)") +
      theme(legend.text.align = 0)
}
```

```{r vol1, fig.cap="Fig 2.1: pH9_Io vs. H9_pBAP_D3_Io", fig.width=8, fig.height=5}
volcanoPlots2(
  res.PH9vsBAP,
  "PH9vsBAP",
  "pH9_Io",
  "H9_pBAP_D3_Io",
  "#0571B0",
  "#483D8B",
  "#4d4d4d",
  ChartTitle = "pH9_Io vs. H9_pBAP_D3_Io"
)
```

```{r vol2, fig.cap="Fig 2.2: H9_BMP4.0h_Krendl vs. H9_BMP4.72h_Krendl", fig.width=8, fig.height=5}
volcanoPlots2(
  res.K00vsK72,
  "K00vsK72",
  "H9_BMP4.0h_Krendl",
  "H9_BMP4.72h_Krendl",
  "#FF1493",
  "#EE82EE",
  "#4d4d4d",
  ChartTitle = "H9_BMP4.0h_Krendl vs. H9_BMP4.72h_Krendl"
)
```

```{r vol3, fig.cap="Fig 2.3: H1_Yabe vs. H1_BAP_D8_>70_Yabe", fig.width=8, fig.height=5}
volcanoPlots2(
  res.UNDvsSTB,
  "UNDvsSTB",
  "H1_Yabe",
  "H1_BAP_D8_>70_Yabe",
  "#598234",
  "#006400",
  "#4d4d4d",
  ChartTitle = "H1_Yabe vs. H1_BAP_D8_>70_Yabe"
)
```

```{r vol4, fig.cap="Fig 2.4: H9_BMP4.72h_Krendl vs. H1_BAP_D8_>70_Yabe", fig.width=8, fig.height=5}
volcanoPlots2(
  res.K72vsSTB,
  "K72vsSTB",
  "H9_BMP4.72h_Krendl",
  "H1_BAP_D8_>70_Yabe",
  "#598234",
  "#EE82EE",
  "#4d4d4d",
  ChartTitle = "H9_BMP4.72h_Krendl vs. H1_BAP_D8_>70_Yabe"
)
```

```{r vol5, fig.cap="Fig 2.5: H9_pBAP_D3_Io vs. H1_BAP_D8_>70_Yabe", fig.width=8, fig.height=5}
volcanoPlots2(
  res.BAPvsSTB,
  "BAPvsSTB",
  "H9_pBAP_D3_Io",
  "H1_BAP_D8_>70_Yabe",
  "#598234",
  "#0571B0",
  "#4d4d4d",
  ChartTitle = "H9_pBAP_D3_Io vs. H1_BAP_D8_>70_Yabe"
)
```

```{r vol6, fig.cap="Fig 2.6: H9_pBAP_D3_Io vs. H9_BMP4.72h_Krendl", fig.width=8, fig.height=5}
volcanoPlots2(
  res.BAPvsK72,
  "BAPvsK72",
  "H9_pBAP_D3_Io",
  "H9_BMP4.72h_Krendl",
  "#EE82EE",
  "#0571B0",
  "#4d4d4d",
  ChartTitle = "H9_pBAP_D3_Io vs. H9_BMP4.72h_Krendl"
)
```

```{r vol7, fig.cap="Fig 2.7: H9_BMP4.72h_Krendl vs. H1_BAP_D8_<40_Yabe", fig.width=8, fig.height=5}
volcanoPlots2(
  res.K72vsL40,
  "K72vsL40",
  "H9_BMP4.72h_Krendl",
  "H1_BAP_D8_<40_Yabe",
  "#AEBD38",
  "#EE82EE",
  "#4d4d4d",
  ChartTitle = "H9_BMP4.72h_Krendl vs. H1_BAP_D8_<40_Yabe"
)
```

```{r vol8, fig.cap="Fig 2.8: H9_pBAP_D3_Io vs. H1_BAP_D8_<40_Yabe", fig.width=8, fig.height=5}
volcanoPlots2(
  res.BAPvsL40,
  "BAPvsL40",
  "H9_pBAP_D3_Io",
  "H1_BAP_D8_<40_Yabe",
  "#AEBD38",
  "#0571B0",
  "#4d4d4d",
  ChartTitle = "H9_pBAP_D3_Io vs. H1_BAP_D8_<40_Yabe"
)
```



## PlacentaCellEnrich (PCE) analyses

The above gene lists are used for running PCE. The function used for running PCE is below.

```{r pceimport}
# Vento-Tormo et al., dataset
l <-
  load(file = "assets/combine-test-expression1.Rdata")
humanGeneMapping <- dataset$GRCH38$humanGeneMapping
d <- dataset$PlacentaDeciduaBloodData
data <- d$expressionData
cellDetails <- d$cellDetails

# Xiang et al., dataset
te.dataset.xiang <- readRDS("assets/te.dataset.xiang.rds")

# Castel et al., dataset
te.dataset.castel <- readRDS("assets/te.dataset.castel.rds")

# full names for cell types
    xi.md <-
  read.csv(
    "assets/md-xi.tsv",
    sep = "\t",
    header = TRUE,
    row.names = 1
  )
vt.md <-
  read.csv(
    "assets/md-vt.tsv",
    sep = "\t",
    header = TRUE,
    row.names = 1
  )
zp.md <-
  read.csv(
    "assets/md-zp.tsv",
    sep = "\t",
    header = TRUE,
    row.names = 1
  )

```


```{r pcefunction}
run.all.PCE <- function(geneList1, geneList2, filename, ChartTitle, barcolor) {
  expressionData <-
    data[intersect(row.names(data), humanGeneMapping$Gene),]
  se <-
    SummarizedExperiment(
      assays = SimpleList(as.matrix(expressionData)),
      rowData = row.names(expressionData),
      colData = colnames(expressionData)
    )
  cellSpecificGenesExp <-
    teGeneRetrieval(se, expressedGeneThreshold = 1)
  print(length(geneList1))
  gs.vt <- GeneSet(geneIds = toupper(geneList1))
  output.vt <- teEnrichmentCustom(gs.vt, cellSpecificGenesExp)
  en.output.vt <-
    setNames(data.frame(assay(output.vt[[1]]), row.names = rowData(output.vt[[1]])[, 1]),
             colData(output.vt[[1]])[, 1])
  row.names(cellDetails) <- cellDetails$RName
  en.output.vt$Tissue <-
    cellDetails[row.names(en.output.vt), "CellName"]
  gs <- GeneSet(unique(geneList2))
  output.xi <- teEnrichmentCustom(gs, te.dataset.xiang)
  output.zp <- teEnrichmentCustom(gs, te.dataset.castel)
  en.output.xi <-
    setNames(data.frame(assay(output.xi[[1]]), row.names = rowData(output.xi[[1]])[, 1]),
             colData(output.xi[[1]])[, 1])
  en.output.xi$Tissue <- rownames(en.output.xi)
  en.output.zp <-
    setNames(data.frame(assay(output.zp[[1]]), row.names = rowData(output.zp[[1]])[, 1]),
             colData(output.zp[[1]])[, 1])
  en.output.zp$Tissue <- rownames(en.output.zp)
  en.output.zp$source <- "ZP"
  en.output.zp <- en.output.zp[order(-en.output.zp$Log10PValue), ]
  en.output.zp <-
    merge(en.output.zp, zp.md, by = "row.names", all.x = TRUE)
  en.output.zp <- rownames_to_column(en.output.zp, var = "Name")
  en.output.vt$source <- "VT"
  en.output.vt <- en.output.vt[order(-en.output.vt$Log10PValue), ]
  en.output.vt <-
    merge(en.output.vt, vt.md, by = "row.names", all.x = TRUE)
  en.output.vt <- rownames_to_column(en.output.vt, var = "Name")
  en.output.xi$source <- "Xi"
  en.output.xi <- en.output.xi[order(-en.output.xi$Log10PValue), ]
  en.output.xi <-
    merge(en.output.xi, xi.md, by = "row.names", all.x = TRUE)
  en.output.xi <- rownames_to_column(en.output.xi, var = "Name")
  en.conbined <- rbind(en.output.vt, en.output.xi, en.output.zp)
  p <- 0.05
  logp <- -log10(p)
  en.conbined <-  en.conbined %>%
    mutate(Log10PValue = replace(Log10PValue, Log10PValue < logp, 0))
  en.conbined %>%
    group_by(source) %>%
    arrange(source, desc(Log10PValue)) %>% dplyr::slice(1:7)  %>%
    ungroup %>%
    mutate(
      source = as.factor(source),
      CellNames = tidytext::reorder_within(CellNames, Log10PValue, source, sep = ":")
    ) %>%
    ggplot(aes(CellNames, Log10PValue)) + geom_bar(stat = 'identity', fill = barcolor) +  theme_minimal() +
    theme(
      axis.text.x = element_text(
        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.title.y = element_blank(),
      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"
      )
    )  +
    scale_y_continuous(expand = expansion(mult = c(0, .1)), breaks = pretty_breaks()) +
    facet_wrap(~ source, scales = "free", ncol = 3) +
    coord_flip() +
    ggtitle(ChartTitle)
}
```
The PCE is run on each of the gene lists as follows (up and down pairs are displayed together).

## PCE plots

```{r pce1, fig.cap="Fig 2.9: PCE results for pH9_Io vs. H9_pBAP_D3_Io", fig.width=14, fig.height=8}
a <-
  run.all.PCE(
    PH9vsBAP.dw.pce1,
    PH9vsBAP.dw.pce2,
    "PH9vsBAP.down_allPCE.v2",
    "Overexpressed in H9_pBAP_D3_Io",
    "#0571B0"
  )
b <-
  run.all.PCE(
    PH9vsBAP.up.pce1,
    PH9vsBAP.up.pce2,
    "PH9vsBAP.up_allPCE.v2",
    "Overexpressed in pH9_Io",
    "#483F8E"
  )
panel_plot <-
  plot_grid(a,
            b,
            labels = c("A", "B"),
            ncol = 1,
            nrow = 2)
panel_plot
```

```{r pce2, fig.cap="Fig 2.10: PCE results for H9_BMP4.0h_Krendl vs. H9_BMP4.72h_Krendl", fig.width=14, fig.height=8}
a <-
  run.all.PCE(
    K00vsK72.dw.pce1,
    K00vsK72.dw.pce2,
    "K00vsK72.down_allPCE.v2",
    "Overexpressed in H9_BMP4.72h_Krendl",
    "#EE82EE"
  )
b <-
  run.all.PCE(
    K00vsK72.up.pce1,
    K00vsK72.up.pce2,
    "K00vsK72.up_allPCE.v2",
    "Overexpressed in H9_BMP4.0h_Krendl",
    "#FF1493"
  )
panel_plot <-
  plot_grid(a,
            b,
            labels = c("A", "B"),
            ncol = 1,
            nrow = 2)
panel_plot
```

```{r pce3, fig.cap="Fig 2.11: PCE results for H1_Yabe vs. H1_BAP_D8_>70_Yabe", fig.width=14, fig.height=8}
a <-
  run.all.PCE(
    UNDvsSTB.dw.pce1,
    UNDvsSTB.dw.pce2,
    "UNDvsSTB.down_allPCE.v2",
    "Overexpressed in H1_BAP_D8_>70_Yabe",
    "#598234"
  )
b <-
  run.all.PCE(
    UNDvsSTB.up.pce1,
    UNDvsSTB.up.pce2,
    "UNDvsSTB.up_allPCE.v2",
    "Overexpressed in H1_Yabe",
    "#006400"
  )
panel_plot <-
  plot_grid(a,
            b,
            labels = c("A", "B"),
            ncol = 1,
            nrow = 2)
panel_plot
```

```{r pce4, fig.cap="Fig 2.12: PCE results for H9_BMP4.72h_Krendl vs. H1_BAP_D8_>70_Yabe", fig.width=14, fig.height=8}
a <-
  run.all.PCE(
    K72vsSTB.dw.pce1,
    K72vsSTB.dw.pce2,
    "K72vsSTB.down_allPCE.v2",
    "Overexpressed in H1_BAP_D8_>70_Yabe",
    "#598234"
  )
b <-
  run.all.PCE(
    K72vsSTB.up.pce1,
    K72vsSTB.up.pce2,
    "K72vsSTB.up_allPCE.v2",
    "Overexpressed in H9_BMP4.72h_Krendl",
    "#EE82EE"
  )
panel_plot <-
  plot_grid(a,
            b,
            labels = c("A", "B"),
            ncol = 1,
            nrow = 2)
panel_plot
```

```{r pce5, fig.cap="Fig 2.13: PCE results for H9_pBAP_D3_Io vs. H1_BAP_D8_>70_Yabe", fig.width=14, fig.height=8}
a <-
  run.all.PCE(
    BAPvsSTB.dw.pce1,
    BAPvsSTB.dw.pce2,
    "BAPvsSTB.down_allPCE.v2",
    "Overexpressed in H1_BAP_D8_>70_Yabe",
    "#598234"
  )
b <-
  run.all.PCE(
    BAPvsSTB.up.pce1,
    BAPvsSTB.up.pce2,
    "BAPvsSTB.up_allPCE.v2",
    "Overexpressed in H9_pBAP_D3_Io",
    "#0571B0"
  )
panel_plot <-
  plot_grid(a,
            b,
            labels = c("A", "B"),
            ncol = 1,
            nrow = 2)
panel_plot
```

```{r pce6, fig.cap="Fig 2.14: PCE results for H9_pBAP_D3_Io vs. H9_BMP4.72h_Krendl", fig.width=14, fig.height=8}
a <-
  run.all.PCE(
    BAPvsK72.dw.pce1,
    BAPvsK72.dw.pce2,
    "BAPvsK72.down_allPCE.v2",
    "Overexpressed in H9_BMP4.72h_Krendl",
    "#EE82EE"
  )
b <-
  run.all.PCE(
    BAPvsK72.up.pce1,
    BAPvsK72.up.pce2,
    "BAPvsK72.up_allPCE.v2",
    "Overexpressed in H9_pBAP_D3_Io",
    "#0571B0"
  )
panel_plot <-
  plot_grid(a,
            b,
            labels = c("A", "B"),
            ncol = 1,
            nrow = 2)
panel_plot
```

```{r pce7, fig.cap="Fig 2.15: PCE results for H9_BMP4.72h_Krendl vs. H1_BAP_D8_<40_Yabe", fig.width=14, fig.height=8}
a <-
  run.all.PCE(
    K72vsL40.dw.pce1,
    K72vsL40.dw.pce2,
    "K72vsL40.down_allPCE.v2",
    "Overexpressed in H1_BAP_D8_<40_Yabe",
    "#AEBD38"
  )
b <-
  run.all.PCE(
    K72vsL40.up.pce1,
    K72vsL40.up.pce2,
    "K72vsL40.up_allPCE.v2",
    "Overexpressed in H9_BMP4.72h_Krendl",
    "#EE82EE"
  )
panel_plot <-
  plot_grid(a,
            b,
            labels = c("A", "B"),
            ncol = 1,
            nrow = 2)
panel_plot
```
```{r pce8, fig.cap="Fig 2.16: PCE results for H9_pBAP_D3_Io vs. H1_BAP_D8_<40_Yabe", fig.width=14, fig.height=8}
a <-
  run.all.PCE(
    BAPvsL40.dw.pce1,
    BAPvsL40.dw.pce2,
    "BAPvsL40.down_allPCE.v2",
    "Overexpressed in H1_BAP_D8_<40_Yabe",
    "#AEBD38"
  )
b <-
  run.all.PCE(
    BAPvsL40.up.pce1,
    BAPvsL40.up.pce2,
    "BAPvsL40.up_allPCE.v2",
    "Overexpressed in H9_pBAP_D3_Io",
    "#0571B0"
  )
panel_plot <-
  plot_grid(a,
            b,
            labels = c("A", "B"),
            ncol = 1,
            nrow = 2)
panel_plot
```


## Session Information

```{r sessioninfo}
sessionInfo()
```