## Executive Summary

During the **[DIYtranscriptomics](https://diytranscriptomics.com/)** curriculum, we analyze human skin biopsies from healthy donors and patients suffering from cutaneous leishmaniasis ([Amorim et al., 2019, *Sci Transl Med*](https://doi.org/10.1126/scitranslmed.aax4204)).

This document synthesizes the condensed **'the essentials'** code chunks from each step of the course (Steps 1, 2, 3, 5, and 6) into a single, publication-grade reproducible report powered by **Quarto**. Every data transformation, quality metric, statistical test, and visualization is generated dynamically from the full underlying count data (`myDGEList`, containing all **34,935** human genes across all 10 libraries).

---

## Step 1 Essentials: Study Design & Data Ingestion

In Step 1, we read in the study design and import our mapped RNA-seq count data. The study includes 10 biological samples: **5 healthy skin controls** (`HS01`–`HS05`) and **5 cutaneous leishmaniasis lesion biopsies** (`CL08`, `CL10`, `CL11`, `CL12`, `CL13`).

We ingest the complete transcriptomic count matrix (`myDGEList`), which was established by mapping with Kallisto against Ensembl GRCh38 and collapsed to gene level with `tximport`.

```{r step1-ingestion, message=FALSE, warning=FALSE}
suppressPackageStartupMessages({
  library(tidyverse)
  library(edgeR)
  library(matrixStats)
  library(cowplot)
  library(DT)
  library(plotly)
  library(limma)
  library(gplots)
  library(RColorBrewer)
})

study_file <- if (file.exists("studydesign.txt")) {
  "studydesign.txt"
} else if (file.exists("Data/files/studydesign.txt")) {
  "Data/files/studydesign.txt"
} else {
  "../../studyDesign.txt"
}

targets <- read_tsv(study_file, show_col_types = FALSE)

if (file.exists("myDGEList")) {
  load("myDGEList")
} else if (file.exists("Data/files/myDGEList")) {
  load("Data/files/myDGEList")
}

targets
```

---

## Step 2 Essentials: Data Wrangling & Normalization

In Step 2, we evaluate expression distributions before and after filtering lowly expressed transcripts, and normalize library composition using the **Trimmed Mean of M-values (TMM)** method.

We retain genes with at least 1 Count Per Million (CPM) in 5 or more libraries (the size of the smallest biological group). We then generate a multi-panel violin plot (`plot_grid(p1, p2, p3)`) comparing raw, filtered non-normalized, and TMM-normalized distributions across all 10 samples.

```{r step2-wrangling, fig.width=12, fig.height=5, fig.cap="Figure 1: Gene expression distributions: (A) unfiltered non-normalized, (B) filtered non-normalized, and (C) filtered TMM normalized."}
sampleLabels <- targets$sample
log2.cpm <- cpm(myDGEList, log=TRUE)

log2.cpm.df <- as_tibble(log2.cpm, rownames = "geneID")
colnames(log2.cpm.df) <- c("geneID", sampleLabels)
log2.cpm.df.pivot <- pivot_longer(log2.cpm.df, # dataframe to be pivoted
                                  cols = -1, # column names to be stored as a SINGLE variable
                                  names_to = "samples", # name of that new variable (column)
                                  values_to = "expression") # name of new variable (column) storing all the values (data)

p1 <- ggplot(log2.cpm.df.pivot) +
  aes(x=samples, y=expression, fill=samples) +
  geom_violin(trim = FALSE, show.legend = FALSE) +
  stat_summary(fun = "median", 
               geom = "point", 
               shape = 95, 
               size = 10, 
               color = "black", 
               show.legend = FALSE) +
  labs(y="log2 expression", x = "sample",
       title="Log2 Counts per Million (CPM)",
       subtitle="unfiltered, non-normalized",
       caption=paste0("produced on ", Sys.time())) +
  theme_bw()

cpm <- cpm(myDGEList)
keepers <- rowSums(cpm>1)>=5 # user defined: retain genes with CPM > 1 in >= 5 samples
myDGEList.filtered <- myDGEList[keepers,]

log2.cpm.filtered <- cpm(myDGEList.filtered, log=TRUE)
log2.cpm.filtered.df <- as_tibble(log2.cpm.filtered, rownames = "geneID")
colnames(log2.cpm.filtered.df) <- c("geneID", sampleLabels)
log2.cpm.filtered.df.pivot <- pivot_longer(log2.cpm.filtered.df, # dataframe to be pivoted
                                           cols = -1, # column names to be stored as a SINGLE variable
                                           names_to = "samples", # name of that new variable (column)
                                           values_to = "expression") # name of new variable (column) storing all the values (data)

p2 <- ggplot(log2.cpm.filtered.df.pivot) +
  aes(x=samples, y=expression, fill=samples) +
  geom_violin(trim = FALSE, show.legend = FALSE) +
  stat_summary(fun = "median", 
               geom = "point", 
               shape = 95, 
               size = 10, 
               color = "black", 
               show.legend = FALSE) +
  labs(y="log2 expression", x = "sample",
       title="Log2 Counts per Million (CPM)",
       subtitle="filtered, non-normalized",
       caption=paste0("produced on ", Sys.time())) +
  theme_bw()

myDGEList.filtered.norm <- calcNormFactors(myDGEList.filtered, method = "TMM")
log2.cpm.filtered.norm <- cpm(myDGEList.filtered.norm, log=TRUE)
log2.cpm.filtered.norm.df <- as_tibble(log2.cpm.filtered.norm, rownames = "geneID")
colnames(log2.cpm.filtered.norm.df) <- c("geneID", sampleLabels)
log2.cpm.filtered.norm.df.pivot <- pivot_longer(log2.cpm.filtered.norm.df, # dataframe to be pivoted
                                                cols = -1, # column names to be stored as a SINGLE variable
                                                names_to = "samples", # name of that new variable (column)
                                                values_to = "expression") # name of new variable (column) storing all the values (data)

p3 <- ggplot(log2.cpm.filtered.norm.df.pivot) +
  aes(x=samples, y=expression, fill=samples) +
  geom_violin(trim = FALSE, show.legend = FALSE) +
  stat_summary(fun = "median", 
               geom = "point", 
               shape = 95, 
               size = 10, 
               color = "black", 
               show.legend = FALSE) +
  labs(y="log2 expression", x = "sample",
       title="Log2 Counts per Million (CPM)",
       subtitle="filtered, TMM normalized",
       caption=paste0("produced on ", Sys.time())) +
  theme_bw()

plot_grid(p1, p2, p3, labels = c('A', 'B', 'C'), label_size = 12)
```

---

## Step 3 Essentials: Multivariate Analysis

In Step 3, we evaluate global transcriptomic variance using Principal Component Analysis (PCA) and construct an interactive data table of normalized expression with group averages (`healthy.AVG`, `disease.AVG`) and `LogFC`.

```{r step3-pca}
group <- targets$group
group <- factor(group)

pca.res <- prcomp(t(log2.cpm.filtered.norm), scale.=F, retx=T)
pc.var <- pca.res$sdev^2 # sdev^2 captures these eigenvalues from the PCA result
pc.per <- round(pc.var/sum(pc.var)*100, 1) 
pca.res.df <- as_tibble(pca.res$x)
pca.plot <- ggplot(pca.res.df) +
  aes(x=PC1, y=PC2, label=sampleLabels, color = group) +
  geom_point(size=4) +
  stat_ellipse() +
  xlab(paste0("PC1 (",pc.per[1],"%",")")) + 
  ylab(paste0("PC2 (",pc.per[2],"%",")")) +
  labs(title="PCA plot",
       caption=paste0("produced on ", Sys.time())) +
  coord_fixed() +
  theme_bw()

ggplotly(pca.plot)
```

```{r step3-datatable}
mydata.df <- mutate(log2.cpm.filtered.norm.df,
                    healthy.AVG = (HS01 + HS02 + HS03 + HS04 + HS05)/5, 
                    disease.AVG = (CL08 + CL10 + CL11 + CL12 + CL13)/5,
                    # now make columns comparing each of the averages above that you're interested in
                    LogFC = (disease.AVG - healthy.AVG)) %>%
  mutate_if(is.numeric, round, 2)

datatable(mydata.df[,c(1,12:14)], 
          extensions = c('KeyTable'), 
          filter = 'top',
          options = list(keys = TRUE, 
                         searchHighlight = TRUE, 
                         pageLength = 10, 
                         lengthMenu = c("10", "25", "50", "100")))
```

---

## Step 5 Essentials: Differential Gene Expression (Limma-Voom)

In Step 5, we identify differentially expressed genes using **Limma-Voom**. We fit linear models with `voom` precision weights, establish the `infection = disease - healthy` contrast, compute empirical Bayes statistics with `eBayes`, and generate an interactive volcano plot across all evaluated genes (`number = 40000`).

We then identify statistically significant DEGs (FDR < 0.01, |logFC| >= 2) using `decideTests` and display the DEG table.

```{r step5-diffgenes}
design <- model.matrix(~0 + group)
colnames(design) <- levels(group)

v.DEGList.filtered.norm <- voom(myDGEList.filtered.norm, design, plot = FALSE)
fit <- lmFit(v.DEGList.filtered.norm, design)
contrast.matrix <- makeContrasts(infection = disease - healthy,
                                 levels=design)

fits <- contrasts.fit(fit, contrast.matrix)
ebFit <- eBayes(fits)
myTopHits <- topTable(ebFit, adjust ="BH", coef=1, number=40000, sort.by="logFC")
myTopHits.df <- myTopHits %>%
  as_tibble(rownames = "geneID")

vplot <- ggplot(myTopHits.df) +
  aes(y=-log10(adj.P.Val), x=logFC, text = paste("Symbol:", geneID)) +
  geom_point(size=2) +
  geom_hline(yintercept = -log10(0.01), linetype="longdash", colour="grey", linewidth=1) +
  geom_vline(xintercept = 1, linetype="longdash", colour="#BE684D", linewidth=1) +
  geom_vline(xintercept = -1, linetype="longdash", colour="#2C467A", linewidth=1) +
  labs(title="Volcano plot",
       subtitle = "Cutaneous leishmaniasis",
       caption=paste0("produced on ", Sys.time())) +
  theme_bw()

ggplotly(vplot, tooltip = "text")
```

```{r step5-deg-table}
results <- decideTests(ebFit, method="global", adjust.method="BH", p.value=0.01, lfc=2)
colnames(v.DEGList.filtered.norm$E) <- sampleLabels
diffGenes <- v.DEGList.filtered.norm$E[results[,1] !=0,]
diffGenes.df <- as_tibble(diffGenes, rownames = "geneID")

datatable(diffGenes.df,
          extensions = c('KeyTable'),
          caption = 'Table 1: DEGs in cutaneous leishmaniasis',
          options = list(keys = TRUE, searchHighlight = TRUE, pageLength = 10, lengthMenu = c("10", "25", "50", "100"))) %>%
  formatRound(columns=c(2:11), digits=2)
```

---

## Step 6 Essentials: Co-expression Modules & Heatmaps

In Step 6, we perform hierarchical clustering on the significant DEGs (`diffGenes`). Rows are clustered using Pearson correlation distance (`1 - cor(t(diffGenes), method = "pearson")`) and columns are clustered using Spearman correlation distance. We partition the genes into two primary co-expression modules (`cutree(clustRows, k = 2)`) and visualize global patterns using `heatmap.2`.

```{r step6-heatmap, fig.width=9, fig.height=8, fig.cap="Figure 4: Clustered heatmap of significant DEGs across healthy controls and leishmaniasis lesions."}
myheatcolors <- rev(brewer.pal(name="RdBu", n=11))
clustRows <- hclust(as.dist(1-cor(t(diffGenes), method="pearson")), method="complete") # cluster rows by pearson correlation
clustColumns <- hclust(as.dist(1-cor(diffGenes, method="spearman")), method="complete")
module.assign <- cutree(clustRows, k=2)
module.color <- rainbow(length(unique(module.assign)), start=0.1, end=0.9) 
module.color <- module.color[as.vector(module.assign)] 

heatmap.2(diffGenes, 
          Rowv=as.dendrogram(clustRows), 
          Colv=as.dendrogram(clustColumns),
          RowSideColors=module.color,
          col=myheatcolors, scale='row', labRow=NA,
          density.info="none", trace="none",  
          cexRow=1, cexCol=1, margins=c(8,20))
```

### Module 2: Upregulated in Disease

```{r step6-module-up, fig.width=9, fig.height=7, fig.cap="Figure 5: Heatmap of Module 2 (Upregulated in Disease)."}
modulePick <- 2 
myModule_up <- diffGenes[names(module.assign[module.assign %in% modulePick]),] 
hrsub_up <- hclust(as.dist(1-cor(t(myModule_up), method="pearson")), method="complete") 

heatmap.2(myModule_up, 
          Rowv=as.dendrogram(hrsub_up), 
          Colv=NA, 
          labRow = NA,
          col=myheatcolors, scale="row", 
          density.info="none", trace="none", 
          RowSideColors=module.color[module.assign%in%modulePick], margins=c(8,20))
```

### Module 1: Downregulated in Disease

```{r step6-module-down, fig.width=9, fig.height=7, fig.cap="Figure 6: Heatmap of Module 1 (Downregulated in Disease)."}
modulePick <- 1 
myModule_down <- diffGenes[names(module.assign[module.assign %in% modulePick]),] 
hrsub_down <- hclust(as.dist(1-cor(t(myModule_down), method="pearson")), method="complete") 

heatmap.2(myModule_down, 
          Rowv=as.dendrogram(hrsub_down), 
          Colv=NA, 
          labRow = NA,
          col=myheatcolors, scale="row", 
          density.info="none", trace="none", 
          RowSideColors=module.color[module.assign%in%modulePick], margins=c(8,20))
```

---

## Reproducibility & Session Info

```{r session-info}
sessionInfo()
```