200415_Australia_HCC_cohort_DEseq.Rmd

---
title: "200415 - Australia HCC cohort DESeq"
output: html_notebook
---

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

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

####Load data:
```{r}
AUS <- read.delim("~/Desktop/170413_htseq_counts_AUS_chr1_22")
cData_AUS <- read.delim("~/Desktop/190822_cData_AUS")
```

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

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

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

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

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

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

####Write the DE analysis results and normalized counts to a table:
```{r}
res.AUS_wco <- results(AUS.deseq, alpha=0.05, lfcThreshold=0,independentFiltering = TRUE, contrast =list("Stagetumor","Stagecontrol"))
write.table(res.AUS_wco, "~/Desktop/170418_deseq_JNS_no07PM0230_wco_LFC0",sep="\t")
write.table(AUS.normalizedCounts, "~/Desktop/170418_libr_size_norm_counts_JNS_no07PM0230_wco_corrected_htseq",sep="\t")
```

####Volcano plot
```{r}
library(ggplot2)
library(ggfortify)

AUSresults <- read.delim("~/Desktop/Livernome_input/170418_deseq_JNS_no07PM0230_wco_LFC0", sep="\t")
RBPs <- read.delim("~/Desktop/Livernome_input/190823_170214_hRBP_final_list.txt")
vp=AUSresults
highlight = read.delim("~/Desktop/Livernome_input/200427_gene_highlight_list.txt", row.names = 1, header=TRUE)

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

vp_sub = subset(vp, rownames(vp) %in% rownames(highlight))
vp_sub = subset(vp_sub, abs(log2FoldChange)>log2(2) & padj<0.05)
vp_sub = merge(x=vp_sub, y=highlight, by.x=0, by.y=0, all.x=TRUE)
#vp_sub = vp_sub[,2:9]

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 = -log10(padj)),
      color = "grey",
      cex = 1
    ) +
    geom_point(
      data = RBP_sub_ns,
      aes(x = log2FoldChange, y = -log10(padj)),
      color = "chartreuse4",
      cex = 1
    ) +
    geom_point(
      data = RBP_sub,
      aes(x = log2FoldChange, y = -log10(padj)),
      fill = "darkgreen",
      color = "black",
      cex = 1.5,
      pch = 21
    ) +
    theme_bw(base_size = 14) +
    ggtitle("Australia 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, 50))

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). 
```