# [7] "YAbR10.OD1.MTB.1.bam.sorted.bam"   "YAbR18.OD1.PACS2.3.bam.sorted.bam"
# [9] "YAbR24.OD3.C3719.3.bam.sorted.bam" "YAbR32.OD3.PA14.2.bam.sorted.bam" 
#[11] "YAbR4.OD1.C3719.1.bam.sorted.bam"  "YAbR11.OD1.MTB.2.bam.sorted.bam"  
#[13] "YAbR1.OD1.B133.1.bam.sorted.bam"   "YAbR25.OD3.CF5.1.bam.sorted.bam"  
#[15] "YAbR33.OD3.PA14.3.bam.sorted.bam"  "YAbR5r.OD1.C3719..bam.sorted.bam" 
#[17] "YAbR13.OD1.PA14..bam.sorted.bam"   "YAbR20r.OD3.B133.2.bam.sorted.bam"
#[19] "YAbR27.OD3.CF5..bam.sorted.bam"    "YAbR34r.OD3.PACS2..bam.sorted.bam"
#[21] "YAbR8.OD1.CF5.2.bam.sorted.bam"    "YAbR14.OD1.PA14..bam.sorted.bam"  
#[23] "YAbR21.OD3.B133.3.bam.sorted.bam"  "YAbR28r.OD3.MTB.1.bam.sorted.bam" 
#[25] "YAbR36.OD3.PACS2.3.bam.sorted.bam" "YAbR9.OD1.CF5.3.bam.sorted.bam"   
#[27] "YAbR17.OD1.PACS2.2.bam.sorted.bam" "YAbR22.OD3.C3719.1.bam.sorted.bam"
#[29] "YAbR30.OD3.MTB.3.bam.sorted.bam"   "YAbR3.OD1.B133.3.bam.sorted.bam"

require("DESeq2")
require("BiocParallel")
register(MulticoreParam(40))
FDR_Threshold=0.05


grps=rbind(c(17,22),c(7,12),c(13,30),c(21,26),c(11,16),c(8,27))
rownames(grps)=c("PA14_OD1","MTB_OD1","B133_OD1","CF5_OD1","C3719_OD1","PACS2_OD1")

grps=rbind(c(10,15),c(24,29),c(18,23),c(14,19),c(9,28),c(20,25))
rownames(grps)=c("PA14_OD3","MTB_OD3","B133_OD3","CF5_OD3","C3719_OD3","PACS2_OD3")

#grps=rbind(c(10,15,17,22),c(7,12,24,29),c(13,18,23,30),c(14,19,21,26),c(9,11,16,28),c(8,20,25,27))
#rownames(grps)=c("pa14","mtb1","b136","cf5","c3719","pacs2")
ng=dim(grps)[1]
nr=dim(grps)[2]
r= read.table("all-ribominus.csv", header=T)
for (g1 in 1:3){
for (g2 in 4:ng){
fn2=paste(rownames(grps)[g1],"_vs_",rownames(grps)[g2],sep="");
print (fn2)
col_ordering=c(grps[g1,],grps[g2,]);
#conditions = (c(rep(rownames(grps)[g1], nr), rep(rownames(grps)[g2], nr)))
#print (c(conditions,col_ordering))
conditions = factor(c(rep(rownames(grps)[g1], nr), rep(rownames(grps)[g2], nr)))
rnaseqMatrix = r[,col_ordering]
rownames(rnaseqMatrix)=r$Geneid
rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,]
coldata=as.data.frame(colnames(rnaseqMatrix));
coldata$condition <- conditions
dds <- DESeqDataSetFromMatrix(countData = rnaseqMatrix,colData = coldata,design = ~ condition)
dds <- DESeq(dds,parallel = T)


if (1==1){
#diff.exp.
res.Condition <- results(dds,contrast=c("condition",rownames(grps)[g1],rownames(grps)[g2]),parallel = T)
fn=paste("MAplot_",fn2,".png",sep="");
png(fn,width=1024,height=1024,pointsize = 22)
plotMA(res.Condition,alpha=0.05,main=fn2,ylab="log2 fold change",ylim=c(-5,5))
dev.off()
resOrdered <- res.Condition[order(res.Condition$padj),]
resOrdered=na.omit (resOrdered)
# show number of sig. genes
print(dim(resOrdered[resOrdered$padj<FDR_Threshold,])[1])
rn=rownames(resOrdered)
c=as.data.frame(coef(dds))
x=cbind(resOrdered,c[rn,])
fn=paste("deseq2_",fn2,".csv",sep="");
write.table(x,fn,sep = "\t",row.names = T,quote=T,col.names=NA)
}
}
}