EpiDriver T-cells

library(limma)
library(gplots) 
library(FactoMineR)

##########################################
# Data importing
setwd("D:/00-LabData/Lab/Experiments/2017-01/EpiDriver/T_cell/")
dta<-read.csv("EpiDriver_Tcells_biasCorrection_Final_StdQt.csv", row.names = 1, header = T)

setwd("D:/00-LabData/Lab/Experiments/2017-01/EpiDriver/BRAIN/")
dta<-read.csv("EpiDriver_Brain_Filtered.csv", row.names = 1, header = T)

dim(dta)
names(dta)

#dta.norm <-(t(t(dta)/rowMeans(t(dta))))  # divide by col mean

boxplot(log2(dta+1),col = c(rep("red",4),rep("darkblue",4),rep("hotpink",4),rep("darkturquoise",4)),
        cex = 1.1, las =2)

dta.cor <- cor(dta, method="pearson", use = "complete")
dta.cor

my_color <- colorRampPalette(c("green","black","red")) 
breaks = unique(c(seq(0,0.4,length=75),seq(0.4,1,length=75)))
heatmap.2(as.matrix(dta.cor),col=my_color,key.title = "Pearson", #breaks = breaks,
          trace="none",density.info="density",symm = T,
          cexRow=1.1,cexCol=1.1,margins=c(7,7),
          Rowv = T, Colv = T) #,
          ColSideColors=c(rep("red",9),rep("darkblue",12),rep("hotpink",12),rep("darkturquoise",9)),
          RowSideColors=c(rep("red",9),rep("darkblue",12),rep("hotpink",12),rep("darkturquoise",9)))
          
          
##########################################
# PCA plot
pca<-PCA(t(dta), scale.unit=T)

plot.PCA(pca,cex=0.8,autoLab="yes",label="ind")
plot.PCA(pca,cex=0.5,autoLab="yes",label="none")

##########################################
# Diff Expression
names(dta)

oHD<-dta[,c(1:4)] # old HD
oWT<-dta[,c(5:8)] # old WT
yHD<-dta[,c(9:12)] # young HD
yWT<-dta[,c(13:16)] # young WT

old<-dta[,c(1:8)] # old
young<-dta[,c(9:16)] # young
wt<-dta[,c(5:8,13:16)] # wildtypes
hd<-dta[,c(1:4,9:12)] # HD

# individual FC
diff_ind_all <- log2(1+dta) - log2(1+rowMeans(yWT))
head(diff_ind_all)

diff_ind_old <- log2(1+old) - log2(1+rowMeans(oWT))
diff_ind_young <- log2(1+young) - log2(1+rowMeans(yWT))   
diff_ind_wt <- log2(1+wt) - log2(1+rowMeans(yWT))  
diff_ind_hd <- log2(1+hd) - log2(1+rowMeans(yWT))  


hm <- heatmap.2( as.matrix(diff_ind_all) )
hc <- as.hclust( hm$rowDendrogram )
mycl<-cutree(hc, h=1.5)
table(mycl)
clusterCols <- rainbow(length(unique(mycl)))
myClusterSideBar <- clusterCols[mycl]

breaks = unique(c(seq(-1,0,length=75),seq(0,1,length=75)))
my_color <- colorRampPalette(c("steelblue","white","red")) 
heatmap.2(as.matrix(diff_ind_all),
          col=my_color,symm=F, key.title = "Rel. to young WT mean (log2)", breaks=breaks,
          trace="none",density.info="density",
          cexRow=0.2,cexCol=1,margins=c(9,9),
          Rowv=as.dendrogram(hc), Colv = F, dendrogram=c("row"),
          ColSideColors=c(rep("red",4),rep("darkblue",4),rep("hotpink",4),rep("darkturquoise",4)),
          RowSideColors= myClusterSideBar) 
myCl <- cbind(clusterID=mycl, diff_ind_all)
myCluster <- myCl[rev(hc$order),]
head(myCluster)

write.csv(myCluster,"myCluster.csv")


# pairwise.wilcox.test
head(diff_ind_all)
grp_name = c(rep("oHD",4),rep("oWT",4),rep("yHD",4),rep("yWT",4))
input <- data.frame(t(rbind(grp_name,diff_ind_all)))
head(input,1)

attach(input)
out<-lapply(input, function(x) pairwise.wilcox.test(as.numeric(x), X1,  p.adj = "none"))
names(out) <- names(input)
out 

oHD <- input[1:4,]
oWT <- input[5:8,]
yHD <- input[9:12,]
yWT <- input[13:16,]

plot(as.numeric(oHD$CBX7.1), as.numeric(yHD$CBX7.1), xlim = c(0,20), ylim = c(0,20))

wilcox.test(as.numeric(ex1$CBX7.1), as.numeric(ex3$CBX7.1),p.adj = "none")
pairwise.wilcox.test(as.numeric(input$CBX7.1), X1, paired=F, p.adj = "none")
t.test(as.numeric(ex1$CBX7.1), as.numeric(ex3$CBX7.1),p.adj = "none")

wilcox<-sapply(out, function(x) {
  p <- x$p.value
  n <- outer(rownames(p), colnames(p), paste, sep='_v_')
  p <- as.vector(p)
  names(p) <- n
  p
})

wilcox.t <- tail(t(wilcox), 306)
head(wilcox.t)


diff_grp_old <- log2(1+rowMeans(oHD)) - log2(1+rowMeans(oWT))
diff_grp_old <- cbind(data.frame(diff_grp_old), wilcox.t[,c(1,2,6,9)])
head(diff_grp_old)

with(diff_grp_old, plot(diff_grp_old,-log10(oWT_v_oHD),pch=19,cex=0.8,col="lightgrey"))
with(subset(diff_grp_old, oWTvoHD<0.05 & abs(diff_grp_old)>1),points(diff_grp_old,-log10(oWTvoHD),pch=20,cex=0.7,col="red"))
with(subset(diff_grp_old, abs(diff_grp_old)>1),points(diff_grp_old,-log10(oWTvoHD),pch=20,cex=0.8,col="red"))
abline(v=c(-1,1),col="blue",lty=2)
abline(v=0,col="blue")
abline(h=(-log10(0.05)),col="blue",lty=2)

write.csv(diff_grp_old,"diff_grp_old.csv")