if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(version = "3.15")

BiocManager::install("tximport")

BiocManager::install("readr")

BiocManager::install("magrittr")
BiocManager::install("DESeq2")
BiocManager::install("RColorBrewer")
BiocManager::install("tidyverse")

BiocManager::install("GenomicFeatures", force = TRUE)
BiocManager::install("pheatmap", force = TRUE)
BiocManager::install("ashr")

#Once you have installed the packages once, you won't need to do this again unless you wipe or clear your working environment
library(tximport)
library(DESeq2)
library(tidyverse)
library(ggpubr)
library(GenomicFeatures)
library(magrittr)
library(readr)
library(ashr)
library(RColorBrewer)
library(pheatmap)

AT_TxDb <- makeTxDbFromGFF("/mnt/Kanta/rna_seq_tutorial/Arabidopsis_thaliana.TAIR10.53.mod.gff3", format = "gff")

k <- keys(AT_TxDb, keytype = "TXNAME")
tx2gene <- AnnotationDbi::select(AT_TxDb, k, "GENEID", "TXNAME")

dir_maf4_ABA <- "/mnt/Kanta/rna_seq_tutorial/quant_output"

samplelist_ABA <- read_delim("/mnt/Kanta/rna_seq_tutorial/ABA_samplelist.txt", 
    delim = "\t", escape_double = FALSE, 
    trim_ws = TRUE)

#head(samplelist_ABA) #double check your sample list

Qfiles_ABA <- file.path(dir_maf4_ABA, samplelist_ABA$Run, "quant.sf")

#head(Qfiles_ABA) #make sure the paths to these files exists and is correct

txi_ABA <- tximport(files = Qfiles_ABA, type = "salmon", tx2gene = tx2gene, countsFromAbundance = "lengthScaledTPM")

#Add more informative experiment identifiers to our counts matrix
colnames(txi_ABA$counts) <- samplelist_ABA$Condition

#Lets save this counts matrix 
write.csv(txi_ABA$counts, 
          file="/mnt/Kanta/rna_seq_tutorial/txi_ABA_counts.csv")

##We can do the same for our abundance matrix
colnames(txi_ABA$abundance) <- samplelist_ABA$Condition
write.csv(txi_ABA$abundance, 
          file="/mnt/Kanta/rna_seq_tutorial/txi_ABA_tpm.csv")

dds_txi_ABA <- DESeqDataSetFromTximport(txi_ABA,
                                   colData = samplelist_ABA,
                                   design = ~Condition)

# Interestd in pre-filtering genes before you run the differential expression? You can keep all genes that have at least 10 counts or more 
#keep <- rowSums(counts(dds_txi_ABA)) >= 10 # store the gene names that have at least 10 counts acorss the entire experiment
#dds_txi_ABA_filt <- ddsTxi_maf4_ABA[keep,] # subset the dataframe keeping only the genes that meet the count threshhold

#maf4_ABA
dds_ABA <- DESeq(dds_txi_ABA)

#Let's see how our samples cluster

vst_ABA <- vst(dds_ABA, blind = FALSE)
vst_ass <- assay(vst_ABA)
plotPCA(vst_ABA , intgroup = c("Condition"))

#Looks like the majority of variance is due to time post-ABA treatment. The PCA only takes the top 10% variable genes. Let's see if we can generate a heatmap that visualizes the clustering of samples using all the genes. 

sampleDists <- dist(t(assay(vst_ABA)))
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- paste(vst_ABA$Condition) #assign informative rownames
colnames(sampleDistMatrix) <- paste(vst_ABA$Run) #assign informative column names
colors <- colorRampPalette(rev(brewer.pal(9, "Blues")))(100)
pheatmap(sampleDistMatrix,
          clustering_distance_rows = sampleDists,
          clustering_distance_cols = sampleDists,
          col=colors)

#There are a few samples of different timepoints that cluster together - this something to check on. But, let's continue for now

#this is where you want to introduce the contrasts that you would like to make
res_ABA_2hr <- results(dds_ABA, contrast = c("Condition", "2hr_ABA", "0hr_ABA"))
#res_ABA_2hr 

res_ABA_4hr <- results(dds_ABA, contrast = c("Condition", "4hr_ABA", "0hr_ABA"))
#res_ABA_4hr 

res_ABA_6hr <- results(dds_ABA, contrast = c("Condition", "6hr_ABA", "0hr_ABA"))
#res_ABA_6hr 

res_ABA_8hr <- results(dds_ABA, contrast = c("Condition", "8hr_ABA", "0hr_ABA"))
#res_ABA_8hr 

topGene <- rownames(res_ABA_2hr)[which.min(res_ABA_8hr$padj)] #grab the gene that has the smallest adjusted p-value for that time-point
plotCounts(dds_ABA, gene = topGene, intgroup = c("Condition")) #plot it 

resultsNames(dds_ABA) #This will give you all the comparisons that you can look at 

res <- lfcShrink(dds_ABA, coef="Condition_8hr_ABA_vs_0hr_ABA", type="ashr") #normalize lfc at a specific timepoint

DESeq2::plotMA(res, xlim=c(5,1e6), ylim=c(-10,10)) #plot your MA plot to visualize the differentially expressed genes

at_ftc <- read.table("/mnt/Kanta/rna_seq_tutorial/at_drna_ftc.out", header=TRUE, row.names=1) 
at_ftc <- at_ftc[ ,6:ncol(at_ftc)] #Only keep the columns containing raw counts for each gene (columns 6,7,8,9)
colnames(at_ftc) <- gsub("\\.[sb]am$", "", colnames(at_ftc)) #Clean up the column names
colnames(at_ftc) <- gsub("align_to_genome.", "", colnames(at_ftc))

at_ftc <- as.matrix(at_ftc) #make this counts table into a matrix
condition <- factor(c(rep("ctl", 2), rep("stress", 2))) #save the conditions into an object 

coldata <- data.frame(row.names=colnames(at_ftc), condition) # Create a coldata frame that maps each sample/column to conditions for DESeq2

dds_heat_mat <- DESeqDataSetFromMatrix(countData=at_ftc, colData=coldata, design=~condition) #generate the dataframe


# Interestd in pre-filtering genes before you run the differential expression? You can keep all genes that have at least 3 counts or more 
#keep <- rowSums(counts(dds_heat_mat)) >= 3 # store the gene names that have at least 3 counts acorss the entire experiment
#dds_heat_mat_filt <- dds_heat_mat[keep,] # subset the dataframe keeping only the genes that meet the count threshhold

dds_at <- DESeq(dds_heat_mat, parallel = T)

vsd_at <- vst(dds_heat_mat, blind = F) 

res_dds_at <- results(dds_at, contrast = c("condition", "stress", "ctl")) 

write.csv(res_dds_at, file = "/mnt/Kanta/rna_seq_tutorial/res_dds_at.csv")

plotPCA(DESeqTransform(vsd_at))

#Majority of our variance is explained by condition!

sampleDists <- dist(t(assay(vsd_at)))
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- paste(vsd_at$condition)
colors <- colorRampPalette( rev(brewer.pal(9, "Reds")) )(100)
pheatmap(sampleDistMatrix,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         col=colors,
         annotation_names_col = TRUE)

#our samples cluster by conditions!

topGene <- rownames(res_dds_at)[which.min(res_dds_at$padj)] #grab the gene that has the smallest adjusted p-value for that time-point
plotCounts(dds_at, gene = topGene, intgroup = c("condition")) #plot it 

resultsNames(dds_at) #This will give you all the comparisons that you can look at 

res <- lfcShrink(dds_at, coef="condition_stress_vs_ctl", type="ashr") #normalize lfc at a specific timepoint

DESeq2::plotMA(res, xlim=c(5,1e6), ylim=c(-10,10)) #plot your MA plot to visualize the differentially expressed genes