##############################################################
# R SCRIPT
##############################################################

 
###  DGE
library(tximport)
library(DESeq2)

files <- c("sample1.genes.results", "sample2.genes.results", "sample3.genes.results", "sample4.genes.results", "sample5.genes.results", "sample6.genes.results")
names(files) <- paste0("sample", 1:6)
labels <- c("wt_1", "wt_2", "wt_3", "rde-1_1", "rde-1_2", "rde-1_3")
names(files) <- paste0(labels)

txi.rsem <- tximport(files, type = "rsem", txIn = FALSE, txOut = FALSE)
txi.rsem$length[txi.rsem$length == 0] <- 1
sampleTable <- data.frame(condition = factor(rep(c("A", "B"), each = 3)))
rownames(sampleTable) <- colnames(txi.rsem$counts)
dds <- DESeqDataSetFromTximport(txi.rsem, sampleTable, ~condition)
dds <- DESeq(dds)
write.csv(counts(dds, normalized=TRUE), 'counts_table-rde-1.csv')


###  Compare rde-1 to wt
res_AB <- results(dds, contrast=c("condition","B", "A"))
resOrdered_AB <- res_AB[order(res_AB$pvalue),]
write.csv(resOrdered_AB, 'rde-1_v_swt.csv')


### Determine significant mRNAs for plotting: p < 0.05, fold change > 1.3
# A (wt) vs B (rde-1)
res_sig_AB <-  as.data.frame(subset(res_AB, padj < 0.05))
res_sig_fc_AB <-  as.data.frame(subset(res_sig_AB, (log2FoldChange < -0.378511623) | (log2FoldChange > 0.378511623)))


### Average counts across samples
cnts <- counts(dds, normalized=TRUE)
cnts_avg <- cnts[,c(1,4)] # just copy to fill values later & maintain rownames
colnames(cnts_avg) <- c("wt","rde-1")
cnts_avg[,1] <- log2(rowMeans(cnts[,1:3]))
cnts_avg[,2] <- log2(rowMeans(cnts[,4:6]))
cnts_avg <- replace(cnts_avg, cnts_avg=='-Inf', -4)

### Get transposon and 23H-RNA identifiers
transposons <- read.csv("transposon_names.txt", row.names=1)
sRNAs <- read.csv("23H-RNA_gene_list.txt", row.names=1)

### Mean scatter plots based on the number of plots
pdf("081224_means_plots-rde-1.pdf")
par(mar=c(5,5,5,5))

# wt vs A
plot_log2_ticks<-function(min,max){
  
  tick_sep <- seq(from=min, to=max, by=4)
  
  for (i in 1:length(tick_sep)){
    #major ticks
    axis(1,at=(tick_sep[i]), labels=prettyNum(tick_sep[i],digits=2,format="f"),las=1,cex.axis=1.0)
    axis(2,at=(tick_sep[i]), labels=prettyNum(tick_sep[i],digits=2,format="f"),las=1,cex.axis=1.0)
  }
  
}
plot(-4:20,-4:20, tck = 0, xaxt='n', yaxt='n', type="n",main="rde-1 vs wild type", xlab = "Log2 average reads in wild type", ylab = "Log2 average reads in rde-1", font.lab=2, las = "1", tcl = -.3, cex.lab = 2, cex.lab = 2, cex.main = 2, cex.axis = 2)
plot_log2_ticks(-4,20)
points(cnts_avg[,1], cnts_avg[,2], pch=20,col="grey80", cex=0.5)
abline(h=log2(10), v=log2(10), lty=2, col="grey60")
abline(0,1, col="grey60")
abline(log2(0.5),1, col="grey60")
abline(log2(2),1, col="grey60")
points(cnts_avg[rownames(res_sig_fc_AB),1],cnts_avg[rownames(res_sig_fc_AB),2], pch=19, col="tomato1",  cex=0.6)
points(cnts_avg[rownames(sRNAs),1],cnts_avg[rownames(sRNAs),2], pch=19,col="deepskyblue",  cex=0.6)
points(cnts_avg[rownames(transposons),1],cnts_avg[rownames(transposons),2], pch=19,col="darkorchid4",  cex=0.6)
d <- c("p < 0.05", "p > 0.05", "23H-RNAs", "transposons")
legend('topleft', legend = d, fill=c("tomato1", "grey80", "deepskyblue", "darkorchid4"), bty='n', cex=1.5)

dev.off()