190822_TCGA_DESeq.Rmd

---
title: "190820 - TCGA DESeq analysis"
output: html_notebook
---

####Load packages:
```{r}
library(DESeq2)
```

####Load packages and settings for multi-core processing
```{r}
library(BiocParallel)
register(MulticoreParam(4))
```

####Load data:
```{r}
TCGA <- read.delim("/Desktop/170418_TCGA_htseqcount_Chr1_22")
cData_TCGA <- read.delim("/Desktop/cData_TCGA.txt")
colnames(TCGA)=rownames(cData_TCGA)
```

####Apply cutoff. Use at least 1 count in 1 or more samples:
```{r}
TCGA_wco=TCGA[(rowSums(TCGA>0)>=1),]
```

####DE analysis with the DESeq2 package (samples are treated as paired)
```{r}
TCGA.deseq <- DESeqDataSetFromMatrix(countData = TCGA_wco, colData = cData_TCGA, design = ~ ID + Stage)
TCGA.deseq <- DESeq(TCGA.deseq, betaPrior=TRUE, parallel=TRUE)
```

####Normalize counts by library size:
```{r}
TCGA.normalizedCounts <- t( t(counts(TCGA.deseq)) / sizeFactors(TCGA.deseq) )
```

####PCA on normalized counts
```{r}
TCGA.PCA<- prcomp(t(log2(TCGA.normalizedCounts+1)))
col.def<-c("blue","red")
colors <- rep(col.def, each=1)
M_TCGA <- summary(TCGA.PCA)
par(pty="s")
plot(TCGA.PCA$x[,1],TCGA.PCA$x[,2],col=colors,pch=20, cex=2, xlab=paste0("PC1 (", round(M_TCGA$importance[2,1]*100), "%)"),ylab=paste0("PC2 (", round(M_TCGA$importance[2,2]*100), "%)"))
plot(TCGA.PCA$x[,3],TCGA.PCA$x[,4],col=colors,pch=20, cex=2, xlab=paste0("PC3 (", round(M_TCGA$importance[2,3]*100), "%)"),ylab=paste0("PC4 (", round(M_TCGA$importance[2,4]*100), "%)"))
```

####Next we library normalize the normalized counts values
```{r}
TCGA.normalizedCounts_colSum=colSums(TCGA.normalizedCounts)
TCGA.normalizedCounts_LibrNorm = sweep(TCGA.normalizedCounts, 2, TCGA.normalizedCounts_colSum, `/`)*(mean(colSums(TCGA.normalizedCounts)))
```

####PCA on library size normalized counts
```{r}
TCGA.PCA<- prcomp(t(log2(TCGA.normalizedCounts_LibrNorm+1)))
col.def<-c("blue","red")
colors <- rep(col.def, each=1)
M_TCGA<-summary(TCGA.PCA)
par(pty="s")
plot(TCGA.PCA$x[,1],TCGA.PCA$x[,2],col=colors,pch=20, cex=2,xlab=paste0("PC1 (", round(M_TCGA$importance[2,1]*100), "%)"),ylab=paste0("PC2 (", round(M_TCGA$importance[2,2]*100),"%)"))
plot(TCGA.PCA$x[,3],TCGA.PCA$x[,4],col=colors,pch=20, cex=2,xlab=paste0("PC3 (", round(M_TCGA$importance[2,3]*100), "%)"),ylab=paste0("PC4 (", round(M_TCGA$importance[2,4]*100),"%)"))
```

####Write the DE analysis results and normalized counts to a table:
```{r}
res.TCGA_wco <- results(TCGA.deseq, alpha=0.05, lfcThreshold=0, independentFiltering = TRUE, contrast = list("StagePrimary_Tumor","StageSolid_Tissue_Normal"))
write.table(res.TCGA_wco, "~/Desktop/170418_deseq_TCGA_wco_LFC0", sep="\t")
write.table(TCGA.normalizedCounts, "~/Desktop/170418_libr_size_norm_counts_TCGA", sep="\t")
```

####load list of genes to highlight
```{r}
highlight = read.delim("~/Desktop/200427_gene_highlight_list.txt", row.names = 1, header=TRUE)
```

####Volcano plot
```{r}
# set values outside axis limits to Inf
upper_padj_limit = 80

TCGAresults_limits = TCGAresults
TCGAresults_limits$log_padj = -log10(TCGAresults_limits$padj)
TCGAresults_limits$log_padj[TCGAresults_limits$log_padj > upper_padj_limit] = Inf
```

```{r}
library(ggplot2)
library(ggfortify)
library(ggrepel)
TCGAresults <- read.delim("~/Desktop/Livernome_input/170418_deseq_TCGA_wco_LFC0", sep="\t")
RBPs <- read.delim("~/Desktop/Livernome_input/190823_170214_hRBP_final_list.txt")
vp=TCGAresults_limits

setwd("~/Desktop")
pdf("200701_TCGA_volcano.pdf", height=3, width=3)
par(pty="s")
par(mar=c(2,2,2,2))

signif_sub = subset(vp, abs(log2FoldChange)>log2(2) & padj<0.05)
RBP_sub = subset(vp, abs(log2FoldChange)>log2(2) & padj<0.05 & row.names(vp)%in%RBPs[,1])
RBP_sub_ns = subset(vp, padj>0.05 & row.names(vp)%in%RBPs[,1])

vp_plot= ggplot(vp) +
  geom_point(
      data = vp,
      aes(x = log2FoldChange, y = log_padj),
      color = "grey",
      cex = 1
    ) +
    geom_point(
      data = signif_sub,
      aes(x = log2FoldChange, y = log_padj),
      color = "red",
      cex = 1
    ) +
    geom_point(
      data = RBP_sub_ns,
      aes(x = log2FoldChange, y = log_padj),
      color = "chartreuse4",
      cex = 1
    ) +
    geom_point(
      data = RBP_sub,
      aes(x = log2FoldChange, y = log_padj),
      fill = "darkgreen",
      color = "black",
      cex = 1.5,
      pch = 21
    ) +
  geom_text_repel(
      data = vp_sub,
      aes(x = log2FoldChange, y = log_padj, label=gene_name),
      size = 5,
      min.segment.length = 0, #use to always put a line
      box.padding = unit(0.35, "lines"),
      point.padding = unit(0.3, "lines")
    ) +
    theme_bw(base_size = 14) +
    ggtitle("TCGA HCC cohort") +
    theme(plot.title = element_text(hjust = 0.5)) +
    scale_x_continuous(limits = c(-7.5, 7.5)) +
    scale_y_continuous(limits = c(0, 80))


vp_plot
print(vp_plot, newpage = FALSE)
dev.off()
```

####Session Info:
```{r}
sessionInfo()
#DEseq2 version used for final DEG result list: 1.14.1 (in bioconductor v. 3.4). 
```