setwd("/nfs/nfs_wooi/Botas/EA12048_botas/work")
source("/nfs/nfs_wooi/WTools/wrs/plotMultiplePCA.R")
source("/nfs/nfs_wooi/WTools/wrs/PCA.SampleNames.Gen.R")
library(DESeq2)

# Read in raw counts
cnt.orig <- read.delim("../Batch.2016-07-12_24089_Supp047/rsem_gcounts_matrix.txt", sep="\t", check.names=FALSE, stringsAsFactors=F)

# Read in sample information
sampleInfo <- read.delim("SampleInfo.Checks.txt", sep="\t", stringsAsFactors=TRUE)
rownames(sampleInfo) <- sampleInfo[,1]
sampleInfo <- sampleInfo[colnames(cnt.orig)[-1], ]

# Generic
# Extract data (g1/g2)
g1="EN11"
g2="EW11"


dds.outf = paste("DESeq_out/", g1, "vs", g2, "_v2_", "dds.saved", sep="")
dds.rld.outf = paste("DESeq_out/", g1, "vs", g2, "_v2_", "dds.rld.saved", sep="")
dds.res.outf = paste("DESeq_out/", g1, "vs", g2, "_v2_", "dds.res.txt", sep="")
dds.PCA.outf = paste("DESeq_out/", g1, "vs", g2, "_v2_", "PCA.pdf", sep="")
dds.PCA2.outf = paste("DESeq_out/", g1, "vs", g2, "_v2_", "PCA_rlog.pdf", sep="")
  
  
sInfo <- sampleInfo[c(grep(g1, sampleInfo$Sample), grep(g2, sampleInfo$Sample)), ]

# Remove the outlier
sInfo <- sInfo[-6, ]

scnt <- cnt.orig[, rownames(sInfo)]
rownames(scnt) <- cnt.orig[ ,1]
  
scond <- data.frame(cond = factor(gsub("(A|B|C)$", "", sInfo$Sample), levels=c(g2, g1)))
dds <- DESeqDataSetFromMatrix(countData = round(scnt), 
                              colData = scond,
                              design = ~cond)
dds.rlog <- rlogTransformation(dds, blind = TRUE)
dds.deseq <- DESeq(dds)
dds.res <- results(dds.deseq, contrast=c("cond", g1, g2))
dds.res$grp.mean <- sapply(levels(dds.deseq$cond), function(lvl) rowMeans(counts(dds.deseq,normalized=TRUE)[,dds.deseq$cond == lvl] ) )
dds.res$norm.counts <- counts(dds.deseq, normalized=TRUE)
colnames(dds.res$norm.counts) <- colnames(scnt)
  
save(dds, file=dds.outf)
save(dds.rlog, file=dds.rld.outf)
write.table(as.data.frame(dds.res), file=dds.res.outf, sep="\t")
  
  ##pdf(file=dds.PCA.outf, width=6, height=5)
PCA.SampleNames.Gen.Simple(as.matrix(dds.res$norm.counts), groups=scond$cond, ntop=nrow(dds.res$norm.counts), 
                           title=paste("PCA ", g1, "_", g2, sep=""), tsize=3, lpos="right", bsize=13)
ggsave(file=dds.PCA.outf, width=6, height=5)
  
pdat <- assay(dds.rlog)
colnames(pdat) <- colnames(scnt)
PCA.SampleNames.Gen.Simple(pdat, groups=scond$cond, ntop=nrow(pdat), 
                           title=paste("PCA rlog ", g1, "_", g2, sep=""), tsize=3, lpos="right", bsize=13)
ggsave(file=dds.PCA2.outf, width=6, height=5)
  

# --
cmps <- read.delim("ComparisonsList.txt", sep="\t", header=F, stringsAsFactors=F)

# Get neuron comparisons
cmps.neuron <- cmps[7:12, ]


all.neuron <- c()
suff <- paste(cmps.neuron[1,2], "vs", cmps.neuron[1,1], sep="")
temp <- read.delim(paste("DESeq_out/", suff, "_dds.res.txt", sep=""), sep="\t", stringsAsFactors=F)[c(2, 5:8)]
colnames(temp) <- paste(suff, "_", colnames(temp), sep="")
all.neuron <- temp
all.neuron$gene <- rownames(all.neuron)

for(i in 2:nrow(cmps.neuron)){
  suff <- paste(cmps.neuron[i,2], "vs", cmps.neuron[i,1], sep="")
  
  if(suff=="EN11vsEW11"){
    suff = "EN11vsEW11_v2"
  }
  
  temp <- read.delim(paste("DESeq_out/", suff, "_dds.res.txt", sep=""), sep="\t", stringsAsFactors=F)[c(2, 5:8)]
  colnames(temp) <- paste(suff, "_", colnames(temp), sep="")
  temp$gene <- rownames(temp)
  
  all.neuron <- join_all(list(all.neuron,temp), by = 'gene', type = 'full')  
}

write.table(all.neuron, file="DESeq_all.neuron_v2.txt", sep="\t")









