DESeq2_gcloud_DEG.R

library(gplots)
library(RColorBrewer)
library(DESeq2)
library(readr)
library(tximport)
library(fdrtool)
options(bitmapType='cairo') #work around for X11 missing on cluster


##############################################################
##############################################################
##############################################################
#functions here
rld_pca <- function (rld, intgroup = "condition", ntop = 500, colors=NULL, legendpos="bottomleft", main="PCA Biplot", textcx=1, ...) {
  require(genefilter)
  require(calibrate)
  require(RColorBrewer)
  rv = rowVars(assay(rld))
  select = order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
  pca = prcomp(t(assay(rld)[select, ]))
  fac = factor(apply(as.data.frame(colData(rld)[, intgroup, drop = FALSE]), 1, paste, collapse = " : "))
  if (is.null(colors)) {
    if (nlevels(fac) >= 3) {
      colors = brewer.pal(nlevels(fac), "Paired")
    }   else {
      colors = c("black", "red")
    }
  }
  pc1var <- round(summary(pca)$importance[2,1]*100, digits=1)
  pc2var <- round(summary(pca)$importance[2,2]*100, digits=1)
  pc1lab <- paste0("PC1 (",as.character(pc1var),"%)")
  pc2lab <- paste0("PC1 (",as.character(pc2var),"%)")
  plot(PC2~PC1, data=as.data.frame(pca$x), bg=colors[fac], pch=21, xlab=pc1lab, ylab=pc2lab, main=main, ...)
  with(as.data.frame(pca$x), textxy(PC1, PC2, labs=rownames(as.data.frame(pca$x)), cex=textcx))
  legend(legendpos, legend=levels(fac), col=colors, pch=20)
  #     rldyplot(PC2 ~ PC1, groups = fac, data = as.data.frame(pca$rld),
  #            pch = 16, cerld = 2, aspect = "iso", col = colours, main = draw.key(key = list(rect = list(col = colours),
  #                                                                                         terldt = list(levels(fac)), rep = FALSE)))
}

volcanoplot <- function (res, lfcthresh=2, sigthresh=0.1, main="Volcano Plot", legendpos="bottomright", labelsig=TRUE, textcx=1, ...) {
  with(res, plot(log2FoldChange, -log10(pvalue), pch=20, main=main, ...))
  with(subset(res, padj<sigthresh ), points(log2FoldChange, -log10(pvalue), pch=20, col="red", ...))
  with(subset(res, abs(log2FoldChange)>lfcthresh), points(log2FoldChange, -log10(pvalue), pch=20, col="orange", ...))
  with(subset(res, padj<sigthresh & abs(log2FoldChange)>lfcthresh), points(log2FoldChange, -log10(pvalue), pch=20, col="green", ...))
  if (labelsig) {
    require(calibrate)
    with(subset(res, padj<sigthresh & abs(log2FoldChange)>lfcthresh), textxy(log2FoldChange, -log10(pvalue), labs=Gene, cex=textcx, ...))
  }
  legend(legendpos, xjust=1, yjust=1, legend=c(paste("FDR<",sigthresh,sep=""), paste("|LogFC|>",lfcthresh,sep=""), "both"), pch=20, col=c("red","orange","green"))
}

#Must be run with resdata, not res! This is because of the lab=Gene
maplot <- function (res, thresh=0.1, labelsig=TRUE, textcx=1, ...) {
  with(res, plot(baseMean, log2FoldChange, pch=20, cex=.5, log="x", ...))
  with(subset(res, padj<thresh), points(baseMean, log2FoldChange, col="red", pch=20, cex=1.5))
  if (labelsig) {
    require(calibrate)
    with(subset(res, padj<thresh), textxy(baseMean, log2FoldChange, labs=Gene, cex=textcx, col=2))
  }
}


entrezid <- function( resdata ) {
  require(EnsDb.Hsapiens.v79)
  a = resdata$Gene #the column to iterate over will be different if I'm using res vs resdata
  tmp=gsub("\\..*","",a)
  tmp <- as.character(tmp)
  txdb <- EnsDb.Hsapiens.v79 
  df <- AnnotationDbi::select(txdb, keys = tmp, keytype = "GENEID", columns = "ENTREZID")
  df2 <- AnnotationDbi::select(txdb, keys = tmp, keytype = "GENEID", columns = "SYMBOL")
  
  ENTREZID <- c() 
  SYMBOL <- c() 
  counter1 <- 0
  for (i in tmp) { 
    counter1 <- counter1 + 1
    j <- match(i,df$GENEID)
    ENTREZID <- c(ENTREZID, toString(df[j,][2]))
    SYMBOL <- c(SYMBOL, toString(df2[j,][2]))}
  resdata$ENTREZID <- ENTREZID
  resdata$SYMBOL <- SYMBOL
  resdata$EnsemblGene <- tmp
  
  resdata
}
##############################################################
##############################################################
##############################################################

library(stringi)

####### RNA-iPOP3 #######
####locate the directory containing the files####pri
dir <- "~/RNA/rsem-results/iPOP3"

###create vector of filenames from table
samples3 <- read.table(file.path(dir, "RNA_Fiber_ipop_3.csv"), sep=",", header = TRUE)

IDs <- substring(as.character(samples3$ID), 2, 9)
IDs_last_8 <- substring(as.character(samples3$ID), 10)
library(Biostrings)
dna = DNAStringSet(stri_reverse(IDs_last_8))
IDs_last_8_rc <- complement(dna)

##import transcript-level estimates from RSEM files
files3 <- file.path(dir, paste0(IDs, IDs_last_8_rc, "_R1_concatenated_", IDs, IDs_last_8_rc, "_R2_concatenated_rsem.genes.results"))

###
#tx2gene <- read.csv(file.path(dir, "tx2gene.csv"))
#head(tx2gene)


####### RNA-iPOP4 #######
####locate the directory containing the files####pri
dir <- "~/RNA/rsem-results/iPOP4"

###create vector of filenames from table
samples4 <- read.table(file.path(dir, "RNA_Fiber_ipop_4.csv"), sep=",", header = TRUE)

IDs <- substring(as.character(samples4$Barcode), 2, 9)
IDs_last_8 <- substring(as.character(samples4$Barcode), 10)
library(Biostrings)
dna = DNAStringSet(stri_reverse(IDs_last_8))
IDs_last_8_rc <- complement(dna)

##import transcript-level estimates from RSEM files
files4 <- file.path(dir, paste0(IDs, IDs_last_8_rc, "_R1_concatenated_", IDs, IDs_last_8_rc, "_R2_concatenated_rsem.genes.results"))

###
#tx2gene <- read.csv(file.path(dir, "tx2gene.csv"))
#head(tx2gene)

####### RNA-iPOP1_2 #######
####locate the directory containing the files####pri
dir <- "~/RNA/rsem-results/rsem_genes_results_all_trimmed"

###create vector of filenames from table
samples <- read.table(file.path(dir, "samples.txt"), header = TRUE)

##import transcript-level estimates from RSEM files
files1_2 <- file.path(dir, samples$run, paste0(samples$run, ".genes.results"))

samplesa <- samples[,-1]
samples4a <- cbind(samples4[,7],samples4[,8:9],samples4[,6])
colnames(samples4a) <- colnames(samplesa)
samples3a <- cbind(samples3[,7],samples3[,8:9],samples3[,6])
colnames(samples3a) <- colnames(samplesa)
###
samplesa <- rbind(samplesa, samples3a, samples4a)
rownames(samplesa) <- paste0("sample", 1:nrow(samplesa))

samplesa <- samplesa[which(samplesa$week %in% c("Baseline","10","20","30","WashoutD3","WashoutD10")),] #getting rid of the weeks that aren't a part of our expiriment.
samplesa <- samplesa[which(samplesa$fiber %in% c("LCInulin","Arabinoxylan")),] #getting rid of the weeks that aren't a part of our expiriment.


files <- c(files1_2,files3, files4)
names(files) <- paste0("sample", 1:length(files))
files <- files[which(names(files) %in% rownames(samplesa))]
txi <- tximport(files, type = "rsem")
names(txi)

library(BiocParallel)
bpparam <- SnowParam(15, log = TRUE, stop.on.error = FALSE)
register(MulticoreParam(15))


###remove zero length transcripts
txi$length[txi$length == 0] <- 1

###create dataframe with proper rownames
sampleTable <- data.frame(condition = factor(samplesa$fiber))
sampleTable[,2] <- data.frame(week = factor(samplesa$week))
rownames(sampleTable) <- colnames(txi$counts)

###construct a DESeqDataSet from the txi object and sample info in samples.txt file
dds <- DESeqDataSetFromTximport(txi, sampleTable, ~ condition + week + condition:week)
###prefilter low count genes
dds <- dds[ rowMeans(counts(dds)) > 0.8, ]

########################################
#fiber_week
########################################

#dds <- DESeqDataSetFromTximport(txi, sampleTable, ~ week:condition)
#sampleTable$fiber_week <- interaction(samples$fiber, samples$week)

###set factor levels
#factor_levels <- unique(sampleTable$condition)
#dds$condition <- factor(dds$condition, levels = factor_levels)
#This loop is slightly clunky, but it works well enough for now. Next time perhaps I will
#create a loop for every factor. That should clean things up nicely.
#The first one will be for "condition" i.e. fiber type

i <- colnames(sampleTable)[3]
  print(i)
j <- "Baseline"
    print(j)
    ddsDESeq <- DESeq(dds, test="LRT", reduced = ~condition + week, parallel = TRUE, BPPARAM=MulticoreParam(15))
    saveRDS(ddsDESeq,file="/home/slancast/RNA/ddstimeseries.rds")
    res <- results(ddsDESeq)
    
      print("LFCShrink")
      #resLFC <- lfcShrink(ddsDESeq, coef=2, res=res,parallel = TRUE, BPPARAM=MulticoreParam(15))
      print("LFCShrinkDone")
      #      names(resLFC)[1] <- "Gene"
      ## Examine independent filtering
      #attr(res, "filterThreshold")
      #plot(attr(res,"filterNumRej"), type="b", xlab="quantiles of baseMean", ylab="number of rejections")
      #png(paste("~/RNA/rsem-results/deseq_results/DEG_resLFC-maplot.png",sep=""), 1500, 1000, pointsize=20)
      #plotMA(resLFC, ylim=c(-1,1), cex=1)
      #dev.off()
      #      res <- subset(res, baseMean > 0.2) #thowing out the extrememly low read count samples.
      
      
      resdata <- merge(as.data.frame(res), as.data.frame(counts(ddsDESeq, normalized=TRUE)), by="row.names", sort=FALSE)
      names(resdata)[1] <- "Gene"
      resdata <- entrezid(resdata)
      ## Write results
      write.csv(resdata, file=paste("~/RNA/rsem-results/deseq_results/DEG_-diffexpr-resultsb.csv",sep=""))
      
      ###Visualization and data collection
      
      # Plot dispersions
      png(paste("~/RNA/rsem-results/deseq_results/DEG_qc-dispersions.png",sep=""), 1000, 1000, pointsize=20)
      plotDispEsts(ddsDESeq, main="Dispersion plot")
      dev.off()
      
      # Regularized log transformation for clustering/heatmaps, etc
      print("rld")
      
      #rld <- bplapply(ddsDESeq, rlogTransformation(ddsDESeq), BPPARAM=MulticoreParam(15)) 
      #head(assay(rld))
      #hist(assay(rld))
      
      # Colors for plots below
      ## Ugly:
      ## (mycols <- 1:length(unique(condition)))
      ## Use RColorBrewer, better looking
      #(mycols <- brewer.pal(8, "Dark2")[1:length(unique(factor_levels))])
      
      # Sample distance heatmap
      #sampleDists <- as.matrix(dist(t(assay(rld))))
      #png(paste("~/RNA/rsem-results/deseq_results/DEG_qc-heatmap-samples.png",sep=""), w=1000, h=1000, pointsize=20)
      #heatmap.2(as.matrix(sampleDists), key=F, trace="none",
                #col=colorpanel(100, "green", "red"),
                #ColSideColors=mycols[ddsDESeq[[i]]], RowSideColors=mycols[ddsDESeq[[i]]],
                #margin=c(10, 10), main="Sample Distance Matrix")
      #dev.off()
      
      # Principal components analysis
      #png(paste("~/RNA/rsem-results/deseq_results/DEG_qc-pca.png",sep=""), 1000, 1000, pointsize=20)
      #rld_pca(rld, colors=mycols, intgroup="condition")
      #dev.off()
      
      
      # Get differential expression results
      ## Examine plot of p-values
      png(paste("~/RNA/rsem-results/deseq_results/DEG_res-pvalue.png",sep=""), 1000, 1000, pointsize=20)
      hist(res$pvalue, breaks=250, col="grey")
      dev.off()
      
      ## MA plot
      png(paste("~/RNA/rsem-results/deseq_results/DEG_diffexpr-maplot-noLFC.png",sep=""), 1500, 1000, pointsize=20)
      maplot(resdata, ylim = c(-1,1), main="MA Plot")
      dev.off()
      
      ## Volcano plot with "significant" genes labeled
      png(paste("~/RNA/rsem-results/deseq_results/DEG_diffexpr-volcanoplot.png",sep=""), 1200, 1000, pointsize=20)
      volcanoplot(resdata,  sigthresh=0.1, textcx=.8, xlim=c(-5, 5))
      dev.off()
    }