TCGAbiolinks - DESeq2 pipeline

TCGAbiolinks - DESeq2 pipeline

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("TCGAbiolinks")

library("AnnotationDbi")
library("DESeq2")
library("calibrate")
library("dplyr")
library("ggplot2")
library("ggpmisc")
library("ggpubr")
library("ggrepel")
library("gplots")
library("org.Hs.eg.db")
library("pcaExplorer")
library("reshape2")
library("TCGAbiolinks")
library("tidyverse")

# Gene expression data for Lung cancer (LUAD, LUSC)
query.COAD <- GDCquery(project = c("TCGA-COAD"), 
                       data.category="Transcriptome Profiling",
                       data.type="Gene Expression Quantification",
                       workflow.type="HTSeq - Counts",
                       sample.type=c("Primary solid Tumor","Solid Tissue Normal")) 
#
GDCdownload(query.COAD, method = "api", files.per.chunk = 5)

#
dta.exp.COAD <- GDCprepare(query.COAD, save=TRUE,
                           save.filename="COAD_geneExp_dataframe.rda",
                           summarizedExperiment=TRUE)
#
geneExp <- SummarizedExperiment::assay(dta.exp.COAD)
head(geneExp[,1:5])
colnames(geneExp)
typeof(geneExp)

# tumor vs normal
tumor<-geneExp[ , grepl( ".............01............." , colnames( geneExp ) ) ]
dim(tumor)
normal<-geneExp[ , grepl( ".............11............." , colnames( geneExp ) ) ]
dim(normal)

cnts <- data.frame(tumor,normal) 

# any columns you want
l <- list(data.frame(geneExp))

#list1 <- read.csv("sampleList.csv",header = F)
list1 <- c("TCGA-DM-A1D4-01A-21R-A155-07","TCGA-NH-A6GB-01A-11R-A37K-07")
list1 <- gsub("-",".",list1)
dta_list1 <- lapply(l, function(x) x%>% select(list1))

#list2 <- read.csv("sampleList.csv",header = F)
list2 <- c("TCGA-DM-A1D4-01A-21R-A155-07","TCGA-NH-A6GB-01A-11R-A37K-07")
list2 <- gsub("-",".",list2)
dta_list2 <- lapply(l, function(x) x%>% select(list2))

dta <- data.frame(dta_list1,dta_list2)
write.csv(dta,"test.csv")

#########################################
cnts_noZero <- cnts[rowMeans(cnts) > 0,]
dim(cnts_noZero)

#########################################
conds <- c(rep("tumor",478),rep("normal",41))

coldat=DataFrame(conds=factor(conds))
coldat

dds <- DESeqDataSetFromMatrix(cnts_noZero, colData=coldat, design = ~ conds)
dds <- DESeq(dds, fitType = "mean")
dds

# to obtain normalized count data
cntNorm <- counts(dds, normalized=T)
write.csv(cntNorm,"~/00-NGS/cntNorm_COAD_478tumor_41normal.csv")

# PCA
#rld <- rlogTransformation(dds,fitType = "mean")
#pcaplot(rld, intgroup=c("conds"), ntop=10000, text_labels=F, 
#        title="PCA", pcX=1, pcY=2, point_size=3, ellipse=F)

# correlation heatmap
dta.cor <- cor(cntNorm)
#dta.cor <- cor(assay(rld))
my_palette <- colorRampPalette(c("steelblue", "lightyellow", "red"))
heatmap.2(dta.cor, density="none", trace="none", col=my_palette, margins=c(8,8),
          ColSideColors=c(rep("orange",478),rep("darkgreen",41)),
          RowSideColors=c(rep("orange",478),rep("darkgreen",41)))

# DEGs
#########################################
group1="tumor"
mean_group1 <- rowMeans(cntNorm[,1:478])
group2="normal"
mean_group2 <- rowMeans(cntNorm[,479:519])

res <- results(dds, contrast = c("conds", group1, group2))

## add group mean expression level
res$group1_name <- group1
res$mean_group1 <- mean_group1
res$group2_name <- group2
res$mean_group2 <- mean_group2

columns("org.Hs.eg.db")
# add gene symbol to res
res$symbol = mapIds(org.Hs.eg.db,
                    keys=row.names(res), 
                    column="SYMBOL",
                    keytype="ENSEMBL",
                    multiVals="first")
# add entrez ID to res
res$refseq = mapIds(org.Hs.eg.db,
                    keys=row.names(res), 
                    column="REFSEQ",
                    keytype="ENSEMBL",
                    multiVals="first")

res_sig <- subset(res, padj <= 1e-5 & abs(log2FoldChange) >= log2(10))           # sig DEGs, p < 0.05, fc > 2
res_sig <- res_sig[order(-res_sig$log2FoldChange),]                         # sort by fold change
head(res_sig)
write.csv(res_sig,"~/sigDEG_COAD_tumor_normal.csv")

#write.csv(res,paste(type,"DESeq2_result_",group1,"_vs_",group2,".csv",sep=""))
#write.csv(res_sig,paste(type,"_sig_DEGs_",group1,"_vs_",group2,".csv",sep=""))

###################################################
nrow1 <- nrow(subset(res, padj<=1e-5 & log2FoldChange >= log2(10)))
nrow1
nrow2 <- nrow(subset(res, padj<=1e-5 & log2FoldChange <= -log2(10)))
nrow2

###################################################
#res$log2FoldChange[abs(res$log2FoldChange) >=10] <- 10
#res$padj[res$padj <= 1e-100] <- 1e-100                  

###################################################
# Volcano plot - DEGs
with(res, plot(log2FoldChange,-log10(padj),pch=19,cex=.8,
               col=rgb(.5,.5,.5,.5),
               main=paste(group2," (",nrow2,")   < >   ",group1," (",nrow1,")",sep="")))
with(subset(res, padj<=0.01 & log2FoldChange>=1),
     points(log2FoldChange,-log10(padj),pch=19,cex=.7,col=rgb(1,.2,.2,.5)))
with(subset(res, padj<=0.01 & log2FoldChange<=-1),
     points(log2FoldChange,-log10(padj),pch=19,cex=.7,col=rgb(.2,.6,.2,.4)))
abline(h=c(-log10(0.01)),col="blue",lty=2)
abline(v=c(-1,1),col="blue",lty=2)
abline(v=0,col="blue")


###################################################
# Heatmap - DEGs

degs <- cntNorm[match(rownames(res_sig), rownames(cntNorm)),]

fc <- log2(degs[,c(1:478)]+1) - log2(rowMeans(degs[,479:519]+1))

my_palette <- colorRampPalette(c("steelblue", "lightyellow", "red"))
#breaks = unique(c(seq(-2,0,length=75), seq(0,2,length=75)))
heatmap.2(fc,density="none",trace="none", col=my_palette, #breaks=breaks,
          Rowv = T,Colv = F, dendrogram="row", margins = c(8,8),
          ColSideColors = c(rep("orange",478),rep("darkgreen",41)))