setwd('/var/www/forge2/production/src/html/files/0x454A546E8E9C11EA8198793F5BFE3F98')
results<-read.table('Unnamed.GWAS.erc.chart.tsv.gz', header=TRUE,sep='\t')
# Class splits the data into non-significant, marginally significant and significant according to 0.05 and 0.01 (in -log10 scale)
results$Class <- cut(results$Pvalue, breaks =c(0, 0.01, 0.05, 1)/length(unique(results[,'Tissue'])), labels=FALSE, include.lowest=TRUE)
# Class splits the data into non-significant, marginally significant and significant according to q-value (B-Y FDR adjusted)
results$Class2 <- cut(results$Qvalue, breaks =c(0, 0.01, 0.05, 1), labels=FALSE, include.lowest=TRUE)
# Re-order the entries according to tissue first and then cell type/line
tissue.cell.order <- unique(results[, c('Tissue', 'Cell')])
tissue.cell.order <- tissue.cell.order[order(tissue.cell.order[,1], tissue.cell.order[,2]), ]
# Collapse into a single string (to support same cell type in different tissues)
tissue.cell.order2 <- apply(tissue.cell.order, 1, paste, collapse = ' -- ')
results$TissueCell <- apply(results[, c('Tissue', 'Cell')], 1, paste, collapse = ' -- ')
results$TissueCell <- factor(results$TissueCell, levels=tissue.cell.order2)
# Plot an empty chart first
pdf('Unnamed.GWAS.erc.chart.pdf', width=22.4, height=8)
ymax = max(-log10(results$Pvalue), na.rm=TRUE)*1.1
ymin = -0.1
par(mar=c(15.5,4,3,1)+0.1)
plot(NA,ylab='', xlab='', main='DMPs in DNase I sites (probably TF sites) in cell lines for erc Unnamed',
ylim=c(ymin,ymax), las=2, pch=19, col = results$Class2, xaxt='n', xlim=c(0,length(levels(results$TissueCell))), cex.main=2)
# Add horizontal guide lines for the Y-axis
abline(h=par('yaxp')[1]:par('yaxp')[2],lty=1, lwd=0.1, col='#e0e0e0')
# Add vertical lines and labels to separate the tissues
tissues <- c(0, cumsum(summary(tissue.cell.order[,'Tissue'])))
abline(v=tissues[2:(length(tissues)-1)]+0.5, lty=6, col='grey')
text((tissues[1:(length(tissues)-1)] + tissues[2:length(tissues)]) / 2 + 0.5, ymax, names(tissues[2:length(tissues)]), col='grey', adj=1, srt=90, cex=1.2)
# Add points (internal color first)
palette(c('red', 'palevioletred1', 'white'))
points(results$TissueCell, -log10(results$Pvalue), pch=19, col = results$Class2, xaxt='n')
# Add contour to the points
palette(c('black', 'palevioletred1', 'steelblue3'))
points(results$TissueCell, -log10(results$Pvalue), pch=1, col = results$Class2, xaxt='n')
# Add X-axis (use cell name only and not TissueCell)
axis(1, seq(1,length(tissue.cell.order[,2])), labels=tissue.cell.order[,2], las=2, cex.axis=0.67)
mtext(1, text='Cell', line=14, cex=1.4)
mtext(2, text='-log10 binomial p-value', line=2, cex=1.4)
# Add legend (internal color first)
palette(c('white', 'palevioletred1', 'red'))
legend('topleft', pch=19, legend=c('q < 0.01', 'q < 0.05', 'non-sig'), col = 3:1, cex=0.8, inset=c(0.001, 0.005), box.col='white', title='FDR q-value', text.col='white', bg='white')
# Add contour to the points in the legend
palette(c('steelblue3', 'palevioletred1', 'black'))
legend('topleft', pch=1, legend=c('q < 0.01', 'q < 0.05', 'non-sig'), col = 3:1, cex=0.8, inset=c(0.001, 0.005), box.col='darkgrey', title='FDR q-value')
palette('default')
dev.off()