args = commandArgs(trailingOnly=TRUE)

FC_CUTOFF=100 # maximal displayed enrichments (higher values, either in enrichment or in depletion, will be roofed to that)
P_CUTOFF=1e-30 # minimal displayed FDR-adjusted p-value (values lower than that will be floored to that)
PRETTY_P_DISPLAY=c(1,1e-10,1e-20,1e-30) # position of ticks on y-axis (for FDR-adjusted p-values)
PVAL_DISPLAY=0.05 # significance cutoff on FDR-adjusted p-value
HIGHLIGHT=3 # highlight (in red) that number of top (i.e.: lowest) p-values and top (i.e.: highest) fold-changes having a p-value<0.05

data=read.table(args[1],sep='\t',skip=11,header=T)
enrichment=as.numeric(sub(' < 0.01',0,data$upload_1..fold.Enrichment.)) # Because https://pantherdb.org replaces all null values (dividing 0 by whatever number) by " < 0.01"
enrichment[enrichment>FC_CUTOFF]=FC_CUTOFF
enrichment[enrichment<1/FC_CUTOFF]=1/FC_CUTOFF
adj_pval=data$upload_1..FDR.
adj_pval[adj_pval<P_CUTOFF]=P_CUTOFF

select=c(1:nrow(data))[adj_pval<PVAL_DISPLAY & enrichment>=tail(sort(enrichment[adj_pval<PVAL_DISPLAY]),n=HIGHLIGHT)[1]] # selecting GO terms with adj_pval<PVAL_DISPLAY and with enrichment being in the top HIGHLIGHT values
select=append(select,c(1:nrow(data))[adj_pval<PVAL_DISPLAY & enrichment<=head(sort(enrichment[adj_pval<PVAL_DISPLAY]),n=HIGHLIGHT)]) # selecting GO terms with adj_pval<PVAL_DISPLAY and with depletion being in the top HIGHLIGHT values
select=append(select,c(1:nrow(data))[adj_pval<=head(sort(adj_pval),n=HIGHLIGHT)]) # selecting GO terms with the HIGHLIGHT lowest adjusted p-values
unselect=setdiff(1:nrow(data),select)

svg(paste('Volcano_plot_from_',args[1],'.svg',sep=''),5,5)
plot(1,ty='n',log='xy',xlab='GO term enrichment',ylab='FDR-adjusted p-value',xlim=c(1/FC_CUTOFF,FC_CUTOFF),ylim=c(P_CUTOFF,1),axes=F)
lines(c(1,1),c(P_CUTOFF,1),col='gray',lty=2)
lines(c(1/FC_CUTOFF,FC_CUTOFF),c(PVAL_DISPLAY,PVAL_DISPLAY),col='gray',lty=2)
par(new=T)
plot(enrichment[unselect],adj_pval[unselect],log='xy',xlab='',ylab='',xlim=c(1/FC_CUTOFF,FC_CUTOFF),ylim=c(P_CUTOFF,1),axes=F)
par(new=T)
plot(enrichment[select],adj_pval[select],log='xy',xlab='',ylab='',xlim=c(1/FC_CUTOFF,FC_CUTOFF),ylim=c(P_CUTOFF,1),axes=F,col='red')
axis(1)
axis(2,labels=PRETTY_P_DISPLAY,at=PRETTY_P_DISPLAY)
#text(enrichment[select],adj_pval[select],data$GO.biological.process.complete[select],col='red',cex=0.15)
dev.off()
out=list(data$GO.biological.process.complete[select],enrichment[select],adj_pval[select])
names(out)=c('GO term','Enrichment','FDR-adj. p-value')
write.csv(as.data.frame(out),paste('Annotations_',args[1],'.csv',sep=''))

