---
title: "RNASeq_Meth_Compare"
author: "HM Putnam"
date: "3/8/2020"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}
library("DESeq2")
library("tidyverse")
library("pheatmap")
library("genefilter")
library("ggplot2")
library("RColorBrewer")
```

#Load Data
```{r}

#treatmentinfo
treatmentinfo <- read.csv("metadata/Ribodepleted_RNASeq_metadata.csv", header = TRUE, sep = ",")
print(treatmentinfo)

MC.treatmentinfo <- subset(treatmentinfo, Species =="Montipora capitata")
MC.gcount <- as.data.frame(read.csv("RNASeq/cleaned_reads/Mcap_gene_count_matrix.csv", row.names="gene_id"))
head(MC.gcount)

#Ensure all sample IDs in colData are also in CountData and match their orders
rownames(MC.treatmentinfo) <- MC.treatmentinfo$Plug.ID
colnames(MC.gcount) <- MC.treatmentinfo$Plug.ID

MCAP.Data <- DESeqDataSetFromMatrix(countData = MC.gcount,
                              colData = MC.treatmentinfo,
                              design = ~1)

PA.treatmentinfo <- subset(treatmentinfo, Species =="Pocillopora acuta")
PA.gcount <- as.data.frame(read.csv("RNASeq/cleaned_reads/Pact_gene_count_matrix.csv", row.names="gene_id"))
head(PA.gcount)

#Ensure all sample IDs in colData are also in CountData and match their orders
rownames(PA.treatmentinfo) <- PA.treatmentinfo$Plug.ID
colnames(PA.gcount) <- PA.treatmentinfo$Plug.ID

PACT.Data <- DESeqDataSetFromMatrix(countData = PA.gcount,
                              colData = PA.treatmentinfo,
                              design = ~1)


```

#### Visualize gene count data
##### Log-transform the count data
First we are going to log-transform the data using a variance stabilizing transforamtion (vst). This is only for visualization purposes. Essentially, this is roughly similar to putting the data on the log2 scale. It will deal with the sampling variability of low counts by calculating within-group variability (if blind=FALSE). Importantly, it does not use the design to remove variation in the data, and so can be used to examine if there may be any variability do to technical factors such as extraction batch effects.

To do this we first need to calculate the size factors of our samples. This is a rough estimate of how many reads each sample contains compared to the others. In order to use VST (the faster log2 transforming process) to log-transform our data, the size factors need to be less than 4. Otherwise, there could be artifacts in our results.

```{r}
MC.SF <- estimateSizeFactors(MCAP.Data) #estimate size factors to determine if we can use vst  to transform our data. Size factors should be less than 4 to use vst
print(sizeFactors(MC.SF)) #View size factors

PA.SF <- estimateSizeFactors(PACT.Data) #estimate size factors to determine if we can use vst  to transform our data. Size factors should be less than 4 to use vst
print(sizeFactors(PA.SF)) #View size factors
```

```{r}
MC.trans.Data <- vst(MCAP.Data, blind=FALSE) #apply a variance stabilizing transforamtion to minimize effects of small counts and normalize wrt library size
head(assay(MC.trans.Data), 3) #view transformed gene count data
MC.gsampleDists<- dist(t(assay(MC.trans.Data))) #calculate distance matix
MC.gsampleDistMatrix <- as.matrix(MC.gsampleDists) #distance matrix
rownames(MC.gsampleDistMatrix) <- colnames(MC.trans.Data) #assign row names
colnames(MC.gsampleDistMatrix) <- NULL #assign col names
colors <- colorRampPalette(rev(brewer.pal(9, "Blues")) )(255) #assign colors
pheatmap(MC.gsampleDistMatrix, #plot matrix
         clustering_distance_rows=MC.gsampleDists, #cluster rows
         clustering_distance_cols=MC.gsampleDists, #cluster columns
         col=colors) #set colors
```

```{r}
PA.trans.Data <- vst(PACT.Data, blind=FALSE) #apply a variance stabilizing transforamtion to minimize effects of small counts and normalize wrt library size
head(assay(PA.trans.Data), 3) #view transformed gene count data
PA.gsampleDists<- dist(t(assay(PA.trans.Data))) #calculate distance matix
PA.gsampleDistMatrix <- as.matrix(PA.gsampleDists) #distance matrix
rownames(PA.gsampleDistMatrix) <- colnames(PA.trans.Data) #assign row names
colnames(PA.gsampleDistMatrix) <- NULL #assign col names
colors <- colorRampPalette(rev(brewer.pal(9, "Blues")) )(255) #assign colors
pheatmap(PA.gsampleDistMatrix, #plot matrix
         clustering_distance_rows=PA.gsampleDists, #cluster rows
         clustering_distance_cols=PA.gsampleDists, #cluster columns
         col=colors) #set colors
```

#Mcap
#### Plot a heatmap of a subset of genes to visualize outliers

```{r}

MC.mat <- assay(MC.trans.Data[c(1:100)]) #make an expression object
MC.mat <- MC.mat - rowMeans(MC.mat) #difference in expression compared to average across all samples
#pdf("MC.GeneExp.Corr.pdf")
pheatmap(MC.mat,
         clustering_distance_rows="euclidean", clustering_method = "average",
         show_rownames =TRUE,
         show_colnames =TRUE,
         cluster_cols = TRUE) 
#dev.off()


```

#Pact

```{r}

PA.mat <- assay(PA.trans.Data[c(1:100)]) #make an expression object
PA.mat <- PA.mat - rowMeans(PA.mat) #difference in expression compared to average across all samples
#pdf("PA.GeneExp.Corr.pdf")
pheatmap(PA.mat,
         clustering_distance_rows="euclidean", clustering_method = "average",
         show_rownames =TRUE,
         show_colnames =TRUE,
         cluster_cols = TRUE) 
#dev.off()


```