Pseudobulk_Limma_Harmony.BatchIndividual_ClusterRes0.1_minPCT0.2.Rmd

---
title: "Differential Expression with Limma"
output: html_document
---

```{r, message=F}
library(Seurat)
library(Matrix)
library(dplyr)
library(edgeR)
library(limma)
library(reshape2)
library(ggplot2)
library(UpSetR)
```


choose parameters (integration type, clustering res, min pct threshold)
```{r}
f<- 'Harmony.Batchindividual'
pct<-0.2
res<- 'SCT_snn_res.0.1'
```

```{r}
path<- here::here("output/DGELists/")
dge<- readRDS(paste0(path,"Pseudobulk_dge_",f, "_", res,"_minPCT",pct,".rds"))
```

```{r}
cpm<- cpmByGroup(dge, group=dge$samples$cluster)
lcpm<- cpmByGroup(dge, group=dge$samples$cluster, log=TRUE)
```

```{r}
hist(lcpm)
```



```{r}
L<- mean(dge$samples$lib.size) *1e-6
M<- median(dge$samples$lib.size) *1e-6
```


```{r}
genes.ribo <- grep('^RP',rownames(dge),value=T)
genes.no.ribo <- rownames(dge)[which(!(rownames(dge) %in% genes.ribo))]
dge$counts <- dge$counts[which(rownames(dge$counts) %in% genes.no.ribo),] #remove ribosomal genes
```

```{r}
dge<- calcNormFactors(dge, method="TMM")

summary(dge$samples$norm.factors)
```

```{r}
design<- model.matrix(~0+ dge$samples$cluster + dge$samples$batch + dge$samples$ind)
```


```{r}
v<- voom(dge, design, plot=TRUE)
v
```


```{r}
fit<- lmFit(v,design)
```

```{r}
dim(fit)
```

```{r}
head(fit)
```

```{r}
#all pairwise cluster comparisons
nclust<- length(unique(dge$samples$cluster))
nterms<- ncol(fit)
contrasts<- NULL
for (b in 1:nclust){
for (i in 1:nclust){
    c<- rep(0,nterms)
    c[b]<- 1
    c[i]<- -1
    contrasts<- cbind(contrasts, c)
}
}

selfcols<- c()
el<- 1
while (length(selfcols) <= nclust){
  selfcols<- c(selfcols, el)
  el<- el + nclust + 1
}

contrasts<- contrasts[,-selfcols]
```

```{r}
contrasts
```


```{r}
#first, testing all pairwise cluster comparisons. 
fit<- contrasts.fit(fit, contrasts= contrasts)
```

```{r}
efit<- eBayes(fit)
```

```{r}
plotSA(efit)
```


```{r}
summary(decideTests(efit))
```


```{r}
qqt(efit$t, df=efit$df.prior+efit$df.residual, pch=16, cex=0.2)
abline(0,1)
```


```{r}
topTable(efit, coef=1, sort.by= "P")
```



```{r}
volcanoplot(efit, coef=13, highlight=5, names=rownames(dge))
```
```{r}
vol<- topTable(efit, coef=13, n=nrow(fit))
labsig<- vol[abs(vol$logFC) >=9.5 | vol$adj.P.Val < 1e-29,]
labsiggenes<- rownames(labsig)
thresh<- vol$adj.P.Val < 0.05
vol<-cbind(vol, thresh)
ggplot(vol, aes(x=logFC, y= -log10(adj.P.Val))) +
  geom_point(aes(colour=thresh), show.legend = FALSE) +
  scale_colour_manual(values = c("TRUE" = "red", "FALSE" = "black")) +
  geom_text(data=labsig, aes(label=labsiggenes))
```

```{r}
vol<- topTable(efit, coef=1, n=nrow(fit))
labsig<-  vol[(vol$adj.P.Val < 5e-29) | (vol$logFC> 0 & vol$adj.P.Val < 1e-22) |vol$logFC> 5 |vol$logFC < -7,]
labsiggenes<- rownames(labsig)
thresh<- vol$adj.P.Val < 0.05
vol<-cbind(vol, thresh)
ggplot(vol, aes(x=logFC, y= -log10(adj.P.Val))) +
  geom_point(aes(colour=thresh), show.legend = FALSE) +
  scale_colour_manual(values = c("TRUE" = "red", "FALSE" = "black")) +
  geom_text(data=labsig, aes(label=labsiggenes))
```

```{r}
dt<- decideTests(efit)

#which genes are always significantly upreg in cluster 0
c0.common<- which(dt[,1]>0 & dt[,2]>0 & dt[,3]>0 & dt[,4]>0 & dt[,5]>0 & dt[,6]>0)

length(c0.common)
```


```{r}
genes.common<- rownames(dt[c0.common,])
```


```{r}
tops0<- topTable(efit, coef=c(1:9), n=nrow(fit))
tops0.common<- tops0[rownames(tops0) %in% genes.common,]

head(tops0.common, n=40)
```

```{r}
"POU5F1" %in% rownames(tops0.common)
```

```{r}
"NANOG" %in% rownames(tops0.common)
```



```{r}
#upsetr plot showing overlaps in DE genes between cluster 0 and all other clusters
c0v1<- as.vector(rownames(dt[which(dt[,1] != 0),]))
c0v2<-as.vector(rownames(dt[which(dt[,2] != 0),]))
c0v3<-as.vector(rownames(dt[which(dt[,3] != 0),]))
c0v4<-as.vector(rownames(dt[which(dt[,4] != 0),]))
c0v5<-as.vector(rownames(dt[which(dt[,5] != 0),]))
c0v6<-as.vector(rownames(dt[which(dt[,6] != 0),]))
c0v7<-as.vector(rownames(dt[which(dt[,7] != 0),]))
c0v8<-as.vector(rownames(dt[which(dt[,8] != 0),]))
c0v9<-as.vector(rownames(dt[which(dt[,9] != 0),]))

List.0<- list(c0vc1=c0v1,c0vc2=c0v2,c0vc3=c0v3,c0vc4=c0v4,c0vc5=c0v5,c0vc6=c0v6,c0vc7=c0v7,c0vc8=c0v8,c0vc9=c0v9)

#List.0<- list(dt[,1:9])
upset(fromList(List.0), nsets = 9)
```

1 v all contrasts

```{r}
fit<- lmFit(v,design)
```

```{r}
contrasts<- NULL
for (i in 1:nclust){
    c<- c(rep(-1,nclust),0,0,0,0)
    c[i]<- nclust-1
    
    contrasts<- cbind(contrasts, c)
}
```

```{r}
contrasts
```

```{r}
#first, testing all pairwise cluster comparisons. 
fit<- contrasts.fit(fit, contrasts= contrasts)
```

```{r}
efit<- eBayes(fit)
```

```{r}
plotSA(efit)
```
```{r}
qqt(efit$t, df=efit$df.prior+efit$df.residual, pch=16, cex=0.2)
abline(0,1)
```

```{r}
c0vall<- topTable(efit, coef=1, n=nrow(fit))

head(c0vall, n=40)
```

```{r}
summary(decideTests(efit))
```


```{r}
up.0<- which(c0vall$adj.P.Val < 0.05 & c0vall$logFC >0)

length(up.0)
```

```{r}
head(c0vall[up.0,], n=40)
```



```{r}
vol<- topTable(efit, coef=1, n=nrow(fit))
labsig<- vol[(vol$adj.P.Val < 5e-28) | (vol$logFC> 0 & vol$adj.P.Val < 5e-25) |vol$logFC> 32 |vol$logFC < -35,]
labsiggenes<- rownames(labsig)
thresh<- vol$adj.P.Val < 0.05
vol<-cbind(vol, thresh)
ggplot(vol, aes(x=logFC, y= -log10(adj.P.Val))) +
  geom_point(aes(colour=thresh), show.legend = FALSE) +
  scale_colour_manual(values = c("TRUE" = "red", "FALSE" = "black")) +
  geom_text(data=labsig, aes(label=labsiggenes))
```


```{r}
c1vall<- topTable(efit, coef=2, n=nrow(fit))

head(c1vall, n=40)
```

```{r}
c2vall<- topTable(efit, coef=3, n=nrow(fit))

head(c2vall, n=40)
```

```{r}
c3vall<- topTable(efit, coef=4, n=nrow(fit))

head(c3vall, n=40)
```

```{r}
c4vall<- topTable(efit, coef=5, n=nrow(fit))

head(c4vall, n=40)
```

```{r}
c5vall<- topTable(efit, coef=6, n=nrow(fit))

head(c5vall, n=40)
```

```{r}
c6vall<- topTable(efit, coef=7, n=nrow(fit))

head(c6vall, n=40)
```

batch + ind comparisons

```{r}
fit<- lmFit(v,design)
```

```{r}
fit<-contrasts.fit(fit, coefficients = (nclust+1):(nclust+4))
```


```{r}
efit<- eBayes(fit)
```

```{r}
plotSA(efit)
```

```{r}
qqt(efit$t, df=efit$df.prior+efit$df.residual, pch=16, cex=0.2)
abline(0,1)
```

```{r}
topsBatch<- topTable(efit, coef=c(1:2),n=nrow(fit))

head(topsBatch, n=40)
```

```{r}
topsInd<- topTable(efit, coef=c(3:4),n=nrow(fit))

head(topsInd, n=40)
```

```{r}
"UTF1" %in% rownames(topsInd)
```

```{r}
topsInd[rownames(topsInd)== "UTF1",]
```

```{r}
topsInd[rownames(topsInd)== "TAC3",]
```


```{r}
plot(density(topsBatch$adj.P.Val))
lines(density(topsInd$adj.P.Val), col= "red")
```

```{r}
boxplot(-log10(topsBatch$adj.P.Val), -log10(topsInd$adj.P.Val), names= c("Batch", "Individual"), main="Distribution of p values from F tests", ylab="-log10(adjusted.p.val)")
```

```{r}
median(topsBatch$F)
```

```{r}
median(topsInd$F)
```