1_DESeq2_batched.Rmd

---
title: "1. RNAseq analysis of FACS-purified ISC/EBs expressing da^RNAi^, da, da:da or sc"
description: "DGE analysis based on DESeq2"
principal investigator: "Joaquín de Navascués"
researchers: "Jerome Korzelius, Aleix Puig-Barbé, Joaquín de Navascués"
output: 
  html_document:
    toc: true
    toc_float: true
    code_folding: hide
    theme: readable
    df_print: paged
    css: doc.css
---
```{r setup, echo=FALSE, cache=FALSE}
ggplot2::theme_set(ggpubr::theme_pubr(base_size=10))
knitr::opts_chunk$set(dev = 'png', 
                      fig.align = 'center', fig.height = 7, fig.width = 8.5, 
                      pdf.options(encoding = "ISOLatin9.enc"),
                      fig.path='notebook_figs/', warning=FALSE, message=FALSE)
```


# 1 Preparation


**Libraries:**
```{r load_libraries, warning=FALSE}
if (!require("librarian")) install.packages("librarian")
librarian::shelf(
  # data
  dplyr, stringr, DESeq2, edgeR, biomaRt, rtracklayer, GenomicFeatures, limma,
  # graphics
  ggplot2, pheatmap, scico, dendsort,
  # convenience
  here, writexl, gzcon)
```

**Set working directory where the script is**
```{r setwd}
if (Sys.getenv("USER")=="JQ") {
  setwd("/Users/JQ/Documents/_CODE REPOS/GitHub/Da_RNAseq")
} else if (Sys.getenv("RSTUDIO")==1) {
  setwd( dirname(rstudioapi::getSourceEditorContext(id = NULL)$path) ) # gets what is in the editor
} else {
  setwd(here::here())
  d <- str_split(getwd(),'/')[[1]][length(str_split(getwd(),'/')[[1]])]
  if (d != 'Da_RNAseq') { stop(
    paste0("Could not set working directory automatically to where this",
           " script resides.\nPlease do `setwd()` manually"))
    }
}
getwd()
```

**Function to calculate Transcripts Per Million counts:**
```{r}
tpm <- function(counts, lengths) {
  return ((counts * 1e6) / (lengths * sum(counts/lengths,na.rm=TRUE)))
}
```

**Path to definitive images (outside repo):**
```{r define_dir2figs}
figdir <- paste0(c(head(str_split(getwd(),'/')[[1]],-1),
                   paste0(tail(str_split(getwd(),'/')[[1]],1), '_figures')),
                 collapse='/')
dir.create(figdir, showWarnings = FALSE)
```


## Experimental conditions


We have five different conditions expressed with the _esg-Gal4, UAS-GFP_ driver stock, in two batches:

* Batch `a`:
  * control
  * _UAS-daughterless_ (_da_)
  * _UAS-da^RNAi^_ (TRiP-line, either JF02488 or JF02092)
* Batch `b`:
  * control
  * _UAS-da:da_
  * _UAS-scute_

The target cells were FAC-sorted, their RNA extracted, reverse-transcribed, amplified and sequenced with Illumina.


## Read genomic features from GTF file


Read the gene coordinates from the GTF file into an R data structure from resources.
Export the transcript lengths for all of the transcripts, from this object.
Use the unique() function to get a vector of all gene IDs.
Get the maximum transcript length for each gene. Maximum gene lengths are required for a gene‐length‐normalisation later in the script.

```{r, echo=FALSE, warning=FALSE, results=FALSE}
gtfurl <- paste0('https://ftp.ensembl.org/pub/release-101/gtf/',
                 'drosophila_melanogaster/',
                 'Drosophila_melanogaster.BDGP6.28.101.gtf.gz')
txdb <- makeTxDbFromGRanges( import( GFFFile(gtfurl) ) ) 
allTranscripts <- transcriptLengths(txdb)
allGeneIDs <- unique(allTranscripts$gene_id)
allGeneLengths <- as.data.frame(allTranscripts %>%
  group_by(gene_id) %>%
  summarize(max.tx_len = max(tx_len)) )
```


## Get gene symbols


Genes are now identified as FlyBase IDs (e.g. FBgn0031208). To get the gene names:
```{r}
ensembl = useEnsembl(biomart = "ENSEMBL_MART_ENSEMBL",
                     dataset="dmelanogaster_gene_ensembl",
                     host = "https://oct2022.archive.ensembl.org")
                     # to update this: https://www.ensembl.org/Help/ArchiveRedirect
filters = listFilters(ensembl) # It defines filters in case you have a specific query
attributes = listAttributes(ensembl) #Defines the features that will be showed

dlist <- getBM(attributes=c('ensembl_gene_id', 'external_gene_name'), mart = ensembl)
rownames(dlist) <- dlist$ensembl_gene_id
dlist[1] <- NULL
names(dlist) <- 'gene_symbol'
write.table(dlist, file="resources/gene_symbols.txt", col.names=NA)
```


## Load raw count data


Read `featureCounts` results for all samples:
```{r}
targets <- read.table("input/targets.txt", header=TRUE, sep="\t")
zipfiles <- unzip("input/featurecounts.zip",list=TRUE)
fcountf <- zipfiles$Name[ grepl('.featurecount$', zipfiles$Name) &
                          !(grepl('^_', zipfiles$Name) )]

rawData <- NULL
# each column of rawData will contain the reads per gene of a sample
counter <- 0
for (fcf in targets$File) {
  if (fcf %in% fcountf) {
    fileContents <- read.table(unzip("input/featurecounts.zip", file=fcf), sep="\t", header=T)

    } else { counter <- counter + 1 }
  rawData <- cbind(rawData, fileContents[,7])
}
if (counter>0) {cat("There is/are ", counter, ' missing featureCount file(s) in the zipped directory')}
unlink('*.featurecount') # cleanup - delete uncompressed files
```

Add column and row names to the `rawData` matrix
```{r}
colnames(rawData) <- paste(targets$Condition, targets$Replicate, targets$Batch, sep='_')
rownames(rawData) <- fileContents$Geneid
```

We want to remove genes with low counts, so we do:
```{r}
cpms <- cpm(rawData)
keep <- rowSums(cpms > 1) >= 3 # detected in at least 3 samples
rawData <- rawData[keep,]
```


# 2 Pipeline


## `DESeq2`-DGE analysis


Create an experimental design object that contains the information from `targets`:
```{r}
exptDesign = data.frame(
  row.names = colnames(rawData),
  condition = targets$Condition,
  batch = targets$Batch)
```

Create `DESeq2DataSet` object containing this experimental design and rawData:
```{r, warning=FALSE}
exptObject <- DESeqDataSetFromMatrix(countData = rawData,
                                     colData = exptDesign,
                                     design = ~ batch + condition) 
exptObject$condition <- relevel(exptObject$condition, ref = "Control") # specifies 'Control' as the reference level
```
`condition` needs to be last in the formula (see [here](http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#differential-expression-analysis) why).


### Exploratory analysis


#### Principal component analysis and batch-correction


Transform the normalized counts and plot them as PCA:
```{r}
vsd_Object <- vst(exptObject, blind=TRUE)
saveRDS(vsd_Object, 'output/vst_pseudocounts.RDS')
plotPCA(vsd_Object)
```

Remove batch effect and re-plot:
```{r}
assay(vsd_Object) <- removeBatchEffect(
  assay(vsd_Object),
  batch=vsd_Object$batch,
  design=model.matrix(~condition, colData(vsd_Object))
  )
saveRDS(vsd_Object, 'output/vst_pseudocounts_batchCorrected.RDS')
plotPCA(vsd_Object)
```

This makes sense: PC1 mostly separates *scute* overexpression from the rest, and PC2 dominates the effect of knocking-down or overexpressing *daughterless* (or the *da:da* variant).


#### Sample correlations


When performing quality assessment of our count data, we need to transform the normalized counts for better visualization of the variance for unsupervised clustering analyses. To assess the similarity of the samples using hierarchical heatmaps, transform the normalized counts and perform hierarchical clustering analysis (with hidden dendrograms)
```{r, fig.height=6}
# this bit of code is a bit longer than it needs to
# as I was experimenting to get the heatmap closer to publication standard

# Compute the correlation values between samples
vsd_cor_Object <- cor(assay(vsd_Object)) 

# heatmap
main.title <- 'RNAseq sample correlations'
## get sorted clusters
sort_hclust <- function(x) as.hclust(dendsort(as.dendrogram(x)))
mat_cluster_cols <- hclust(dist(t(vsd_cor_Object)))
mat_cluster_cols <- sort_hclust(mat_cluster_cols)
mat_cluster_rows <- hclust(dist(vsd_cor_Object))
mat_cluster_rows <- sort_hclust(mat_cluster_rows)
## mark the batches
annot_batch <- data.frame(batch = ifelse(test = targets$Batch == 'a',
                                         yes = 'batch A',
                                         no = 'batch B'))
rownames(annot_batch) <- rownames(vsd_cor_Object)
## get minimum correlation value, rounded for the legend
bot <- ceiling(min(vsd_cor_Object)*100)/100
## plot
pheatmap(
  # data
    mat               = vsd_cor_Object,
    scale             = "none", # otherwise numbers are changed
    cellwidth         = 15,
    cellheight        = 15,
  # title
    main              = main.title,
    fontsize          = 14,
    annotation        = dplyr::select(exptDesign, condition),
  # rows
    cluster_rows      = mat_cluster_rows,
    treeheight_row    = 25, # default is 50
    show_rownames     = TRUE,
    labels_row        = rownames(exptDesign),
    fontsize_row      = 9,
    annotation_row    = annot_batch,
  # cols
    cluster_cols      = mat_cluster_cols,
    treeheight_col    = 25,
    show_colnames     = TRUE,
    labels_col        = rownames(exptDesign),
    fontsize_col      = 9,
    angle_col         = 45,
  # legends
    legend_breaks     = c(bot, 1),
  # tiles
    color             = scico(255, palette='bamako'),
    border_color      = 'grey80')
```

This makes sense, so we will save the `vsd_Object` for the descriptive visualisations.
```{r}
saveRDS(vsd_Object, 'output/vsd.RDS')
```


### Differential expression analysis


This will normalise the data, correct for dispersion (variance between replicates) and set data up for a differential comparison of any 2 conditions.
```{r}
analysisObject = DESeq(exptObject)
```


#### Dispersion estimates


After fitting the model with the previous command, let us explore the fit of our data to the negative binomial model by plotting the dispersion estimates using the `plotDispEsts()` function. Remember that the dispersion estimates are used to model the raw counts; if the dispersions do not follow the assumptions made by DESeq2, then the variation in the data could be poorly estimated and the DE results could be less accurate.

The assumptions DESeq2 makes are that the dispersions should generally decrease with increasing mean and that they should more or less follow the fitted line.
```{r}
plotDispEsts(analysisObject)
```

This seems reasonably in order, so we move on. Now we just need from `DESeq2` the "`p-values`" and the "`log2 fold changes`" using `results`.
We will first obtain the "`tmp`" values for the gene-based plots, and include them alongside the raw and DESeq2-normalised counts, and the DGE data, in the supplementary data - that should be enough for anyone to do their own light analysis without having to get the reads from GEO.


#### Save counts for visualisation


This looks reasonable, so let us move on by extracting the rawCounts and the normalisedCounts:
```{r}
rawCounts <- as.data.frame(counts(analysisObject, normalized=FALSE))
normalisedCounts <- as.data.frame(counts(analysisObject, normalized=TRUE))
```

Also obtain the TPMs as a normalisation value:
```{r}
# add column with transcript lengths, turn GeneID rownames into column
rawDataWithLengths <- merge(allGeneLengths, rawCounts, by.x="gene_id", by.y="row.names", all=T)
rawCountData <- rawDataWithLengths[,colnames(rawCounts)]
rownames(rawCountData) <- rawDataWithLengths[,1]

# create matrix of TPM values per sample
tpmData <- NULL
for (colName in colnames(rawCountData)) {
    tpmData <- cbind(tpmData, tpm(rawDataWithLengths[,colName], rawDataWithLengths$max.tx_len))
}
# turn TPM into dataframe
tpmData <- as.data.frame(tpmData)
colnames(tpmData) <- colnames(rawCounts)
rownames(tpmData) <- rawDataWithLengths[,1]
# reduce TPM to the genes detected in the RNAseq samples
tpmNormalisedCounts <- tpmData[match(rownames(rawCounts), rownames(tpmData)), ]

# test that the rows are still the same:
if (!identical(rownames(rawCounts), rownames(normalisedCounts))) {
    stop()
}
if (!identical(rownames(tpmNormalisedCounts), rownames(normalisedCounts))) {
    stop()
}
# test that there is information for all of them (e.g. no mismatch between the GTF files used)
cat(
  cat('There are\t',nrow(rawCounts),
      '\tgenes listed in `rawCounts`, and there are data for\t\t',
      nrow(na.omit(rawCounts)), '\tof them.\n'),
  cat('There are\t',nrow(normalisedCounts),
      '\tgenes listed in `normalisedCounts`, and there are data for\t',
      nrow(na.omit(normalisedCounts)), '\tof them.\n'),
  cat('There are\t',nrow(tpmNormalisedCounts),
      '\tgenes listed in `tpmNormalisedCounts`, and there are data for\t',
      nrow(na.omit(tpmNormalisedCounts)), '\tof them.\n')
)
```

All seems fine, so we can store these data. For visualisation:
```{r}
# for downstream use
saveRDS(rawCounts, 'output/rawCounts.RDS')
saveRDS(normalisedCounts, 'output/normalisedCounts.RDS')
saveRDS(tpmNormalisedCounts, 'output/tpmNormalisedCounts.RDS')
```


#### Save DGE data


Loop over the experimental conditions and save `DESeq2::results` as RDS and in Supplementary Table S3:
```{r, warning=FALSE}
# to get a more informative naming for the samples:
targets$sampleIDs <- names(rawCounts)
# conditions to be tested
test_conditions <- unique( targets[targets$Condition != 'Control',]$Condition )
test_names <- paste0(rep('Control_vs_',length(test_conditions)),test_conditions)
tests <- as.list(rep(NA, length(test_names)))
names(tests) <- test_names 

for (condtn in test_conditions) {
  # get the Counts for those conditions
  deData <- as.data.frame(results(analysisObject,
                                  contrast=c("condition", condtn, 'Control'), # Reference goes last!
                                  pAdjustMethod="BH"))
  # add column of ID
  deData <- cbind(data.frame('ensemblGeneID'=rownames(deData)), deData)
  # sort by pval
  deData <- deData[order(deData$pvalue), ] 
  # add gene symbol column and reorder columns
  deData <- merge(deData, dlist, by=0)
  deData <- deData[,c(1,ncol(deData),2:(ncol(deData)-1))]
  # save for later
  saveRDS(deData, file=paste0("output/", 'Control_vs_', condtn, ".RDS"))
  # save as Supplementary data for publication
  cols <- c('gene_symbol', 'ensemblGeneID', 'log2FoldChange', 'padj')
  rownames(deData) <- NULL
  tests[[paste0('Control_vs_', condtn)]] <- deData %>% dplyr::select(all_of(cols))
}
```
Note that in the parameter `contrast`, `'Control'` is last - I cannot find formal justification of this, but even with the `relevel`ling of the `condition` column of the experimental design object to make 'Control' the reference level, this needs to go last, or the log~2~FC values will have the opposite sign.


#### Produce Supplementary material


Supplementary Table S3:
```{r}
# add `rawCounts` to `tests`
tests <- rlist::list.append(tests, `Raw counts per gene per sample`=rawCounts)
# Supplementary data for publication
write_xlsx(tests, path = 'output/Table S3.xlsx')
```

Finally, to have the experimental designed captured in a flexible manner:
```{r}
targets$condition_md <- plyr::mapvalues(
  targets$Condition,
  from=unique(targets$Condition),
  to=c('*da^RNAi^*', '*da*', '*wild-type*', '*da:da*', '*scute*')
  )
targets$condition_md <- factor(
  targets$condition_md,
  c('*wild-type*', '*da*', '*da:da*','*da^RNAi^*', '*scute*')
  )
saveRDS(targets, 'output/targets.RDS')
```

We now move on to preparing descriptive figures. After, to gene set analysis.