#
# Code fragment (require loading libraries)
# How to get differential expression results 
# from input table of raw counts
#

PROJECT <- read.table("myproject.tsv", header=TRUE, sep="\t")
data <- read.table("rawcounts_complete_fastq.tsv", header=TRUE, sep="\t", row.names=1)
row.names(data) <- sub("gene:","",row.names(data))
dds <- DESeqDataSetFromMatrix(data, colData=PROJECT, design = ~ conditions)
dds <- DESeq(dds)
resP3vsP5 <- results(dds, contrast = c("conditions", "P3C004", "P5C004")) 
finalP3vsP5 <- data.frame(
							GENEID = row.names(resP3vsP5),
							log2BaseMeanP3vsP5 = log2(resP3vsP5$baseMean),
							log2RatioP3vsP5 = resP3vsP5$log2FoldChange,
							Stderr_log2RatioP3vsP5 = resP3vsP5$lfcSE,
							pvalue = resP3vsP5$pvalue,
							padjust = resP3vsP5$padj,
							assays(dds)$counts
						)

annots <- read.table("mart_export.txt", header=TRUE, sep="\t", quote="", check.names=FALSE)
finalP3vsP5 <- merge(finalP3vsP5, annots, by.x = "GENEID", by.y = "Gene stable ID", all.x = TRUE)
write.table(finalP3vsP5, file="FIESTABLE_pH3C004_vs_PH5C004_complete_fastq.tsv", sep="\t", row.names=FALSE) 

#
# Proposed tasks:
# calculate differential expression for "pH 5 high CO2" vs "pH5 low CO2"
# calculate again from counts obtained ignoring duplicated reads
#