cols=c("#1B9E77", "#D95F02", "#7570B3", "#E7298A", "#66A61E", "#E6AB02", "#A6761D",  "#666666")

  hlagenes <- c("HLA-F","HLA-V","HLA-P","HLA-G","HLA-H",
                "HLA-T","HLA-K","HLA-U","HLA-A","HLA-W","HLA-Y",
                "HLA-J","HLA-L","HLA-N","HLA-E","HLA-C","HLA-B",
                "HLA-S","MICA","MICB", "C4A", "C4B",
                "HLA-DRA","HLA-DRB3",
                "HLA-DRB4","HLA-DRB5","HLA-DRB1","HLA-DQA1",
                "HLA-DQB1","HLA-DQB2", "HLA-DQA2","HLA-DOB","TAP2","TAP1",
                "HLA-DMB","HLA-DMA","HLA-DOA", "HLA-DPA1",
                "HLA-DPB1","HLA-DPA2","HLA-DPB2")

  classI <- c("HLA-A","HLA-B","HLA-C","HLA-E","HLA-F","HLA-G",
              "HLA-H","HLA-HFE","HLA-J","HLA-K","HLA-L","HLA-N",
              "HLA-P","HLA-S","HLA-T","HLA-U","HLA-V","HLA-W","HLA-Y")


  classII <- c("HLA-DMA","HLA-DMB","HLA-DOA","HLA-DOB","HLA-DPA1",
               "HLA-DPA2","HLA-DPB1","HLA-DPB2","HLA-DQA1","HLA-DQA2",
               "HLA-DQB1","HLA-DRA","HLA-DRB1","HLA-DRB2","HLA-DRB3",
               "HLA-DRB4","HLA-DRB5","HLA-DRB6","HLA-DRB7",
               "HLA-DRB8","HLA-DRB9", "HLA-DRB345")

  pseudogenes <- c("HLA-DPA2", "HLA-DPA3", "HLA-DPB2", "HLA-DRB2",
                   "HLA-DRB6", "HLA-DRB7", "HLA-DRB8", "HLA-DRB9",
                   "HLA-H", "HLA-J", "HLA-K", "HLA-L", "HLA-N", "HLA-P",
                   "HLA-S", "HLA-T", "HLA-U", "HLA-V", "HLA-W",
                   "HLA-X", "HLA-Y", "HLA-Z", "KIR2DP1", "KIR3DP1")

da<-read.table("div.csv", header=T)

da$stars <- NA

o = da$right + da$mid/2 < 5000
da$stars[o] = -log10((da$right + da$mid/2)/100000 * dim(da)[1])[o]

afr <- da[da$category == 'afr', 2:9]
eas <- da[da$category == 'eas', 2:9]
btw <- da[da$category == 'btw', 2:9]

rownames(afr) <- afr$gene
rownames(eas) <- eas$gene
rownames(btw) <- btw$gene

afr = afr[hlagenes, ]
eas = eas[hlagenes, ]
btw = btw[hlagenes, ]


ylim=c(0,5)
L = dim(afr)[1]
xlim=c(1,L)
pdf("hla.pdf", width=12, height=4)
par(mar=c(5,1,2,5))
plot(NA,NA, xlim=xlim, ylim=ylim, axes=F, xlab="", ylab="", xpd=F)
grid()
points(1:L-0.25, afr$ob*100, type='h', lwd=4, col=cols[1])
points(1:L+0.00, eas$ob*100, type='h', lwd=4, col=cols[2])
points(1:L+0.25, btw$ob*100, type='h', lwd=4, col=cols[3])
axis(4, las=1, cex.axis=1.5)
mtext("Pairwise Difference (%)", side = 4, line = 3, cex=1.5)

intrepret_star <-function(star)
{
  if(is.na(star)) return(NA)
  x = abs(star)
  if(star > 0) {
    if (x >=3) return("***")
    if (x >=2) return("*")
    if (x >=1.3) return("*")
  }
  return(NA)
}

for(i in 1:L){

  btwptxt = intrepret_star(btw$stars[i])
  if(!is.na(btwptxt)) text(btwptxt, x= i+0.25, y = btw$ob[i]*100 + 0.15, adj=0, srt=90, xpd=T,
                           cex=1, col=cols[3], lwd=2)
}
labels = rownames(afr)
for(i in 1:length(labels)){
  if (labels[i] %in% pseudogenes) labels[i] = paste0('*',labels[i])
}
axis(1, labels= F, at=1:L)
text(x = 1:L,
     y = par("usr")[3] - 0.35,
     labels = labels,
     xpd = NA,
     ## Rotate the labels by 35 degrees.
     srt = 45,
     adj = 1,
     cex = 1.2)

legend("topleft", 
       legend=c("HPRC-AFR", "CPC-p1", "Between"), 
       col=cols[1:3], cex=1.2,
       lwd=4)
box()
dev.off()
