MA plots

library(cummeRbund)

# set working directory

setwd('/home/bio3/00-NGS/RNAseq/E-MTAB-3037_human_aging_fibroblasts/diffout')

# read cuffdiff results
cuff_data<-readCufflinks()   

# DEGs
diff <- diffData(genes(cuff_data),'Young','Old') 
nrow(diff)
head(diff)

diff_p05 <- subset(diff, p_value <0.05)
nrow(diff_p05)
head(diff_p05)

diff_q05 <- subset(diff, q_value <0.05)

log2_plus <- subset(diff, log2_fold_change > 0)

avg=((diff$value_1+diff$value_2)/2)
avg_p05=((diff_p05$value_1+diff_p05$value_2)/2)
avg_q05=((diff_q05$value_1+diff_q05$value_2)/2)
avg_p05_log2_plus=((log2_plus$value_1+log2_plus$value_2)/2)

## MA plot
with(diff, plot(log10(avg),log2_fold_change,pch=20,cex=0.6,
    main="Young vs Old",col="grey",xlim=c(-4,5)))  
with(subset(diff,p_value<.05),points(log10(avg_p05),log2_fold_change,
    pch=20,cex=0.7,col="red"))
with(subset(log2_plus,log2_fold_change > 0),
    points(log10(avg_p05_log2_plus),log2_fold_change,pch=20,cex=0.7,col="blue"))
with(subset(diff,q_value<.05),points(log10(avg_q05),log2_fold_change,
    pch=20,cex=0.7,col="chartreuse3"))
abline(h=c(-1,1),col="black",lty = 5)
abline(h=c(0,0),col="grey")
abline(h=c(-5,5),col="grey")
abline(h=c(-10,10),col="grey")
abline(v=c(-2,2),col="grey")
abline(v=c(-4,4),col="grey")