library(DESeq2)

coldatas <- read.delim('samples.tsv', row.names=1, header=TRUE, sep='\t')
coldatas

allcounts <- read.delim('all.counts.10-MM.tsv', row.names=1, header=TRUE, sep= '\t')
head(allcounts)
dim(allcounts)

comparisons = list(
		c('myb30','Col0'),
		c('myb92','Col0'),
		c('ors1','Col0'),
		c('tcp11','Col0'),
		c('tga1tga4','Col0'),
		c('wrky31','Col0'),
		c('wrky60','Col0')
		)

for (comparison in comparisons){
  comp1 <- comparison[1]
  comp2 <- comparison[2]
  comp <- paste(comp1,comp2, sep='_vs_')
  print(comp)
  compregex <- paste(comp1,comp2, sep='|')
  keepcomp <- grep(compregex,rownames(coldatas))
  counts <- allcounts[,keepcomp]
  head(counts)
  coldata <- coldatas[keepcomp,]
  print(coldata)

  dds <- DESeqDataSetFromMatrix(countData=as.matrix(counts[,rownames(coldata)]), colData = coldata,
                                  design = ~ 1)
  dds$genotype <- factor(dds$genotype)
  dds$treatment <- factor(dds$treatment)
  design(dds) <- ~ genotype + treatment + genotype:treatment
  dds <- DESeq(dds)

  res <- results(dds, name = paste('genotype',comp1,'.treatmentPEG', sep=''))
  write.table(res[order(res$padj),], file = paste('PEG_vs_mock.',comp1,'_vs_Col0-specific.DESeq2.tsv', sep=''),
		sep='\t', quote = FALSE, row.names = TRUE, col.names = NA)

}