#### colors for CDR3-IMGT lengths

colorlistCDR3lengths=c("1"="#F5F5F5","2"="#D3D3D3","3"="#C0C0C0","4"="#000000","5"="#BBACAC",
                       "6"="#A9A9A9","7"="#808080","8"="#696969","9"="#800000","10"="#B22222",
                       "11"="#B15928","12"="#BC8F8F","13"="#E41A1C","14"="#FF4500","15"="#FF7F00",
                       "16"="#FDBF6F","17"="#FFD700","18"="#FFFF00","19"="#FFFF99","20"="#FFB6C1",
                       "21"="#FB9A99","22"="#FF69B4","23"="#FF1493","24"="#FF00FF","25"="#C805FF",
                       "26"="#9932CC", "27"="#984EA3","28"="#CAB2D6", "29"="#A6CEE3","30"="#559EA0",
                       "31"="#377eb8","32"="#4169E1","33"="#0000FF","34"="#00BFFF","35"="#00FFFF",
                       "36"="#ADFF2F","37"="#90EE90","38"="#00FF00", "39"="#32CD32","40"="#006400")

#### column name for the set1 and set2 displays

Newcolname<-c("CDR3-IMGT sequence (AA)","Experimental ID","Representative sequence index",
              "Nb of \'1 copy\'","Nb of \'More than one\'","Total nb of \'1 copy\' and \'More than one\'",
              "\'1 copy\' Indexes","V-gene","V-allele","D-gene","D-allele",
              "J-gene","J-allele","CDR1-IMGT","CDR2-IMGT","CDR1-IMGT gapped sequence (AA)",
              "CDR2-IMGT gapped sequence (AA)","V-REGION % identity","Sequence length","C104","F118 or W118","Anchors","Sequence ID","Functionality","Sequence file number","sequence clonotype file number")

### clonRem the function to remove CDR3-IMGT with outlier lengths

clonRem <-function(set, min=4, max=45){
  set1rem<- subset(set,((nchar(as.character(set[,1]))<=max)==TRUE & (nchar(as.character(set[,1]))>=min)==TRUE))
  return(set1rem)
}

######
NewcolnameTabperGene<-c("CDR3-IMGT sequence (AA)","Experimental ID","Representative sequence index",
                        "Nb of \'1 copy\'","Nb of \'More than one\'","Total nb of \'1 copy\' and \'More than one\'",
                        "\'1 copy\' Indexes","V-gene","V-allele","D-gene","D-allele",
                        "J-gene","J-allele","CDR1-IMGT","CDR2-IMGT","CDR1-IMGT gapped sequence (AA)",
                        "CDR2-IMGT gapped sequence (AA)","V-REGION % identity","Sequence length","C104","F118 or W118","Anchors","Sequence ID","Functionality","Normalized.nb.sequences","CDR1-IMGT length","CDR2-IMGT length","CDR3-IMGT length","ID with V-J alleles and CDR3-IMGT sequence (AA)")

###### Completebase is the function that generate the column ID, valleles, jalleles and IDcomp

Completebase <- function(data,dataAll){
  report <-data
  report$Normalized.nb.sequences<-((report$total)/(sum(dataAll$total)))*10000
  report$CDR1length<-nchar(as.vector(report$cdr1))
  report$CDR2length<-nchar(as.vector(report$cdr2))
  report$CDR3length<-nchar(as.vector(report$cdr3aa))
  report$ID<-rownames(report)
  report$valleles<-sapply(strsplit(as.character(report$vallele), " "), function(x) x[[2]])
  report$jalleles<-sapply(strsplit(as.character(report$jallele), " "), function(x) x[[2]])
  IDcomp <- function(x,y,z,t){paste(x,"- (",y,", ",z,") ",t,sep = "")}
  report$IDC <- mapply(IDcomp,x=report$ID,y=report$valleles,z=report$jalleles,t=report$cdr3aa)
  return(report[c(1:24,27,28,29,30,34)])}

#################Functions and colors to obtain the graph of the presence of AA at each position of CDR-IMGT
AA<- c("A","R","N","D","C","Q","E","G","H","I","L","K","M","F","P","S","T","W","Y","V","X","\\*")
AAcodonStop<- c("A","R","N","D","C","Q","E","G","H","I","L","K","M","F","P","S","T","W","Y","V","X","*")
AAplot<- c('\\*', 'X', 'V', 'Y', 'W', 'T', 'S', 'P', 'F', 'M', 'K', 'L', 'I', 'H', 'G', 'E', 'Q', 'C', 'D', 'N', 'R', 'A')
AAcodonStopplot<- c('*', 'X', 'V', 'Y', 'W', 'T', 'S', 'P', 'F', 'M', 'K', 'L', 'I', 'H', 'G', 'E', 'Q', 'C', 'D', 'N', 'R', 'A')
classPhysico<- c("Aliphatic","Basic","Sulfur","Hydroxyl","Acidic","Amide","G","F","P","W","Y","X","*")
classHydro<- c("Hydrophobic","Neutral","Hydrophilic","X","*")
classVolume<- c("Very small","Small","Medium","Large","Very large","X","*")
classChemical<- c("Aliphatic","Aromatic","Sulfur","Hydroxyl","Basic","Acidic","Amide","X","*")
classCharge <-  c("Positive charged","Negativ charged","Uncharged","X","*")
classHydroDonorAA <- c("Donor","Acceptor","Donor and acceptor","None","X","*")
classPolarity<- c("Polar","Nonpolar","X","*")

colors=c("A"="#CCFFFF","R"="#E60606","N"="#FF9900","D"="#FFCC99","C"="#00FFFF","Q"="#FF6600",
         "E"="#FFCC00","G"="#00FF00","H"="#FFFF99","I"="#000080","L"="#3366FF","K"="#C64200",
         "M"="#99CCFF","F"="#00CCFF","P"="#FFFF00","S"="#CCFF99","T"="#00FF99","W"="#CC99FF",
         "Y"="#CCFFCC","V"="#0000FF","X"="#FFFFFF","*"="#000000")
#,"#"="#DDDDDD" frame shift color

colorsPhysico<- c("Aliphatic"="#1B04AC","Basic"="#EC1504","Sulfur"="#CCECFF","Hydroxyl"="#89F88B",
                  "Acidic"="#FFCC00","Amide"="#CCA504","G"="#00FF00","F"="#00CCFF","P"="#FFFF00",
                  "W"="#CC99FF","Y"="#CCFFCC","X"="#FFFFFF","*"="#000000")
colorsHydro<- c("Hydrophobic"="#3366FF","Neutral"="#FFFF00","Hydrophilic"="#E60606","X"="#FFFFFF","*"="#000000")

colorsVolume<- c("Very small"="#FF9999","Small"="#E60606","Medium"="#FFFF00","Large"="#99CCFF","Very large"="#3366FF","X"="#FFFFFF","*"="#000000")

colorsChemical<- c("Aliphatic"="#1B04AC","Aromatic"="#CC99FF","Sulfur"="#CCECFF","Hydroxyl"="#89F88B","Basic"="#EC1504",
                   "Acidic"="#FFCC00","Amide"="#CCA504","X"="#FFFFFF","*"="#000000")

colorsCharge <-  c("Positive charged"="#EC1504","Negativ charged"="#FFCC00","Uncharged"="#CCECFF","X"="#FFFFFF","*"="#000000")

colorsHydroDonorAA <- c("Donor"="#00CCFF","Acceptor"="#FFFF00","Donor and acceptor"="#89F88B","None"="#D6D6D6","X"="#FFFFFF","*"="#000000")

colorsPolarity<- c("Polar"="#E60606","Nonpolar"="#3366FF","X"="#FFFFFF","*"="#000000")


frequences<-function(y,x){A <- unlist(lapply(gregexpr(y,x),function(x) if(x[[1]]!=-1) length(x) else 0))
return(A)
}

# Positions in CDR1-IMGT
dt1=11
bt1=5
first1=27
last1=first1+dt1
num_cdr1<-function(i,l){
  i=i-1
  m=min(((l-1)%/%2),bt1)
  n=l-min(l-m,dt1-bt1)
  if (i<=m){pos=as.character(first1+i)} else if (i>=n){pos=as.character(last1-(l-1-i))} else {t=i>=m+(n-m)/2
  pos=paste(first1+b1+if(t){1} else{0},".", if(t){n-i} else{i-m},sep="")}
  return (pos)
}
# Positions in CDR2-IMGT
dt2=9
bt2=4
first2=56
last2=first2+dt2
num_cdr2<-function(i,l){
  i=i-1
  m=min(((l-1)%/%2),bt2)
  n=l-min(l-m,dt2-bt2)
  if (i<=m){pos=as.character(first2+i)} else if (i>=n){pos=as.character(last2-(l-1-i))} else {t=i>=m+(n-m)/2
  pos=paste(first2+bt2+if(t){1} else{0},".", if(t){n-i} else{i-m},sep="")}
  return (pos)
}
#Positions in CDR3-IMGT
d=12
b=6
first=105
last=first+d
num_cdr<-function(i,l){
  i=i-1
  m=min(((l-1)%/%2),b)
  n=l-min(l-m,d-b)
  if (i<=m){pos=as.character(first+i)} else if (i>=n){pos=as.character(last-(l-1-i))} else {t=i>=m+(n-m)/2
  pos=paste(first+b+if(t){1} else{0},".", if(t){n-i} else{i-m},sep="")}
  return (pos)
}


####### Shanon, Kabat index for AA and CDR-IMGT Positions
Shannon<-function (data){
  H<-numeric(nrow(data))
  len<-numeric(length(data[1,]))
  for(i in 1:nrow(data)){
    prop<-data[i,]/sum(data[i,])
    for(j in 1:ncol(data)){ len[j]<-prop[j]*log(prop[j])    }
    H[i]<- -sum(len, na.rm=T)}
  return(round(H, 3))}

Kabat<-function (data){
  N<-colSums(data[, drop=FALSE])
  k <-c()
  for(i in 1:ncol(data)){k[i] <- length(which(data[,i, drop=FALSE] != 0))}
  n <-c()
  for(j in 1:ncol(data)){n[j]<-max(data[,j])}
  Variability <- (N*k)/n
  round(Variability, 3)
}

Simpson<-function(data){
  SimpsonDiv<-numeric(ncol(data))
  for(i in 1:ncol(data)){
    s<-sum(data[,i])
    prop<-(data[,i]/s)^2
    D<-sum(prop)
    SimpsonDiv[i]<-1-D
  }

  return(round(SimpsonDiv, 3))
}



##################################

normjuxBarsAllele<-function(data,...){

  dataV <- data[which(data$Allele_Type=="V allele"),]
  data1V=melt(dataV[c(1,2,5,8)], id.vars = c("Allele_Name", "Allele_Type"))
  pBV<-with(data1V, ggplot(data=data1V) +theme_grey()+
              geom_bar(aes(x=Allele_Name,y =value,fill =factor(variable, levels = rev(levels(variable)))), position="dodge",stat = "identity", colour="black",show.legend=FALSE)+
              geom_bar(aes(x=Allele_Name,y =value,fill =factor(variable, levels = rev(levels(variable)))), position="dodge",stat = "identity")+
              coord_flip() +xlab("V alleles")+ylab("Normalized proportions")+
              scale_x_discrete(limits= rev(factor(as.factor(dataV$Allele_Name),levels=levels(as.factor(dataV$Allele_Name)))))+              theme(
                panel.grid.major = element_line(size = 0.5),
                panel.border = element_rect(linetype = "solid", colour = "#808080", fill=NA),
                legend.position="bottom",
                plot.title=element_text( face="bold", size = unit(14, "pt")),
                strip.text.x = element_text(size=unit(10, "pt"), colour="black"),
                axis.title.x = element_text(size = unit(10, "pt"),face="bold", colour="black",vjust=-0.5),
                axis.text.x=element_text(size = unit(10, "pt"), colour="black"),
                strip.text.y = element_text(size=unit(10, "pt"), face="bold", colour="black"),
                axis.title.y = element_text(size = unit(10, "pt"),face="bold", colour="black", vjust=1.5),
                axis.text.y=element_text( size = unit(10, "pt"),hjust=0, colour="black"),
                strip.background = element_rect( colour="#808080", fill="#DBDBDB", size=1, linetype="solid"),
                legend.background = element_rect(colour = "black"),
                legend.title = element_text(size=unit(8, "pt"), face="bold"),
                legend.text = element_text(size = unit(8, "pt"), face = "bold",vjust=1.5),
                legend.key.size = unit(0.35, "cm"))+
              guides(fill=guide_legend(title = "Compared sets: ",nrow=2,keywidth =1.5, keyheight = 1,reverse = TRUE)))

  dataJ <- data[which(data$Allele_Type=="J allele"),]
  data1J=melt(dataJ[c(1,2,5,8)], id.vars = c("Allele_Name", "Allele_Type"))
  pBJ<-with(data1J, ggplot(data=data1J) +theme_grey()+
              geom_bar(aes(x=Allele_Name,y =value,fill =factor(variable, levels = rev(levels(variable)))), position="dodge",stat = "identity", colour="black",show.legend=FALSE)+
              geom_bar(aes(x=Allele_Name,y =value,fill =factor(variable, levels = rev(levels(variable)))), position="dodge",stat = "identity")+
              coord_flip() +xlab("J alleles")+ylab("Normalized proportions")+
              scale_x_discrete(limits= rev(factor(as.factor(dataJ$Allele_Name),levels=levels(as.factor(dataJ$Allele_Name)))))+
              theme(
                panel.grid.major = element_line(size = 0.5),
                panel.border = element_rect(linetype = "solid", colour = "#808080", fill=NA),
                legend.position="bottom",
                plot.title=element_text( face="bold", size = unit(14, "pt")),
                strip.text.x = element_text(size=unit(10, "pt"), colour="black"),
                axis.title.x = element_text(size = unit(10, "pt"),face="bold", colour="black",vjust=-0.5),
                axis.text.x=element_text(size = unit(10, "pt"), colour="black"),
                strip.text.y = element_text(size=unit(10, "pt"), face="bold", colour="black"),
                axis.title.y = element_text(size = unit(10, "pt"),face="bold", colour="black", vjust=1.5),
                axis.text.y=element_text( size = unit(10, "pt"),hjust=0, colour="black"),
                strip.background = element_rect( colour="#808080", fill="#DBDBDB", size=1, linetype="solid"),
                legend.background = element_rect(colour = "black"),
                legend.title = element_text(size=unit(8, "pt"), face="bold"),
                legend.text = element_text(size = unit(8, "pt"), face = "bold",vjust=1.5),
                legend.key.size = unit(0.35, "cm"))+
              guides(fill=guide_legend(title = "Compared sets: ",nrow=2,keywidth =1.5, keyheight = 1,reverse = TRUE)))

  dataD <- data[which(data$Allele_Type=="D allele"),]
  data1D=melt(dataD[c(1,2,5,8)], id.vars = c("Allele_Name", "Allele_Type"))
  pBD<-with(data1D, ggplot(data=data1D) +theme_grey()+
              geom_bar(aes(x=Allele_Name,y =value,fill =factor(variable, levels = rev(levels(variable)))), position="dodge",stat = "identity", colour="black",show.legend=FALSE)+
              geom_bar(aes(x=Allele_Name,y =value,fill =factor(variable, levels = rev(levels(variable)))), position="dodge",stat = "identity")+
              coord_flip() +xlab("D alleles")+ylab("Normalized proportions")+
              scale_x_discrete(limits= rev(factor(as.factor(dataD$Allele_Name),levels=levels(as.factor(dataD$Allele_Name)))))+              theme(
                panel.grid.major = element_line(size = 0.5),
                panel.border = element_rect(linetype = "solid", colour = "#808080", fill=NA),
                legend.position="bottom",
                plot.title=element_text( face="bold", size = unit(14, "pt")),
                strip.text.x = element_text(size=unit(10, "pt"), colour="black"),
                axis.title.x = element_text(size = unit(10, "pt"),face="bold", colour="black",vjust=-0.5),
                axis.text.x=element_text(size = unit(10, "pt"), colour="black"),
                strip.text.y = element_text(size=unit(10, "pt"), face="bold", colour="black"),
                axis.title.y = element_text(size = unit(10, "pt"),face="bold", colour="black", vjust=1.5),
                axis.text.y=element_text( size = unit(10, "pt"),hjust=0, colour="black"),
                strip.background = element_rect( colour="#808080", fill="#DBDBDB", size=1, linetype="solid"),
                legend.background = element_rect(colour = "black"),
                legend.title = element_text(size=unit(8, "pt"), face="bold"),
                legend.text = element_text(size = unit(8, "pt"), face = "bold",vjust=1.5),
                legend.key.size = unit(0.35, "cm"))+
              guides(fill=guide_legend(title = "Compared sets: ",nrow=2,keywidth =1.5, keyheight = 1,reverse = TRUE)))
  list(BarGphV=pBV,BarGphD=pBD,BarGphJ=pBJ)}



#######################
diffpropGphAllele<-function(data,...){
  colors=c("Null difference in proportions  "="#FFFFFF", "Non-significant (rawp)  "="#E40000", "Significant (rawp)  "= "cyan",
           "Significant (All_p)  "="#0026FF", "Significant (Min_2p)  "="#FF00DC","Significant (Only_BH)  "= "green3")

  pbrut <-with(data, ggplot(data,aes(x=factor(Allele_Name, levels=rev(factor(Allele_Name))),y =Difference_proportion))+
                 geom_errorbar(aes(col = Test_interpretation, ymin =Lower_bound_IC_diff_prop,
                                   ymax =Upper_bound_IC_diff_prop), cex=0.75, width=0.6)+ theme_grey()+coord_flip()+
                 geom_point(aes(shape = Test_interpretation), size = 1) +ylim( round(min(data$Lower_bound_IC_diff_prop)-0.01,2),  round(max(data$Upper_bound_IC_diff_prop)+0.01,2)) +
                 geom_hline(aes(yintercept = 0)) +
                 theme(
                   panel.grid.major = element_line(size = 0.5),
                   panel.border = element_rect(linetype = "solid", colour = "#808080", fill=NA),
                   legend.position="bottom",
                   plot.title=element_text( face="bold", size = unit(14, "pt")),
                   strip.text.x = element_text(size=unit(10, "pt"), colour="black"),
                   axis.title.x = element_text(size = unit(10, "pt"),face="bold", colour="black",vjust=-0.5),
                   axis.text.x=element_text(size = unit(10, "pt"), colour="black"),
                   strip.text.y = element_text(size=unit(10, "pt"), face="bold", colour="black"),
                   axis.title.y = element_text(size = unit(10, "pt"),face="bold", colour="black", vjust=1.5),
                   axis.text.y=element_text( size = unit(10, "pt"),hjust=0, colour="black"),
                   strip.background = element_rect( colour="#808080", fill="#DBDBDB", size=1, linetype="solid"),
                   legend.background = element_rect(colour = "black"),
                   legend.title = element_text(size=unit(8, "pt"), face="bold"),
                   legend.text = element_text(size = unit(8, "pt"), face = "bold",vjust=1.5),
                   legend.key.size = unit(0.4, "cm"))+
                 scale_color_manual(name="Test interpretation: ",
                                    values=colors,
                                    breaks=c("Null difference in proportions  ", "Non-significant (rawp)  ","Significant (rawp)  ","Significant (All_p)  ","Significant (Min_2p)  ","Significant (Only_BH)  "))+
                 scale_shape_manual(name="Test interpretation: ",values=c(19,19,19,19,19,19),
                                    breaks=c("Null difference in proportions  ","Non-significant (rawp)  ",  "Significant (rawp)  ","Significant (All_p)  ","Significant (Min_2p)  ","Significant (Only_BH)  ")))

  Datav <- data[which(data$Allele_Type=="V allele"),]

  pv <- with(Datav,ggplot(Datav,aes(x=factor(Allele_Name, levels=rev(factor(Allele_Name))), y =Difference_proportion))+
               geom_errorbar(aes(col = Test_interpretation, ymin =Lower_bound_IC_diff_prop, ymax =Upper_bound_IC_diff_prop), cex=0.75, width=0.6)+ theme_grey()+coord_flip()+
               geom_point(aes(shape = Test_interpretation), size = 1) +
               ylab("Difference in proportions & 95% CI")+
               ylim( round(min(data$Lower_bound_IC_diff_prop)-0.01,2),  round(max(data$Upper_bound_IC_diff_prop)+0.01,2)) +
               geom_hline(aes(yintercept = 0)) + xlab("V alleles")+
               theme(
                 panel.grid.major = element_line(size = 0.5),
                 panel.border = element_rect(linetype = "solid", colour = "#808080", fill=NA),
                 legend.position="bottom",
                 plot.title=element_text( face="bold", size = unit(14, "pt")),
                 strip.text.x = element_text(size=unit(10, "pt"), colour="black"),
                 axis.title.x = element_text(size = unit(10, "pt"),face="bold", colour="black",vjust=-0.5),
                 axis.text.x=element_text(size = unit(10, "pt"), colour="black"),
                 strip.text.y = element_text(size=unit(10, "pt"), face="bold", colour="black"),
                 axis.title.y = element_text(size = unit(10, "pt"),face="bold", colour="black", vjust=1.5),
                 axis.text.y=element_text( size = unit(10, "pt"),hjust=0, colour="black"),
                 strip.background = element_rect( colour="#808080", fill="#DBDBDB", size=1, linetype="solid"),
                 legend.background = element_rect(colour = "black"),
                 legend.title = element_text(size=unit(8, "pt"), face="bold"),
                 legend.text = element_text(size = unit(8, "pt"), face = "bold",vjust=1.5),
                 legend.key.size = unit(1.5, "cm"))+
               scale_color_manual(name="Test interpretation: ",
                                  values=colors,
                                  breaks=c("Null difference in proportions  ","Non-significant (rawp)  ", "Significant (rawp)  ", "Significant (All_p)  ",
                                           "Significant (Min_2p)  ","Significant (Only_BH)  "))+
               scale_shape_manual(name="Test interpretation: ",values=c(19,19,19,19,19,19),
                                  breaks=c("Null difference in proportions  ","Non-significant (rawp)  ",  "Significant (rawp)  ",
                                           "Significant (All_p)  ", "Significant (Min_2p)  ", "Significant (Only_BH)  ")))

  Dataj <- data[which(data$Allele_Type=="J allele"),]

  pj <- with(Dataj,ggplot(Dataj,aes(x=factor(Allele_Name, levels=rev(factor(Allele_Name))),
                                    y =Difference_proportion))+
               geom_errorbar(aes(col = Test_interpretation, ymin =Lower_bound_IC_diff_prop,
                                 ymax =Upper_bound_IC_diff_prop), cex=0.75, width=0.6)+ theme_grey()+coord_flip()+
               geom_point(aes(shape = Test_interpretation), size = 1) + ylab("Difference in proportions & 95% CI")+
               ylim( round(min(data$Lower_bound_IC_diff_prop)-0.01,2),  round(max(data$Upper_bound_IC_diff_prop)+0.01,2)) +
               geom_hline(aes(yintercept = 0)) + xlab("J alleles")+
               theme(
                 panel.grid.major = element_line(size = 0.5),
                 panel.border = element_rect(linetype = "solid", colour = "#808080", fill=NA),
                 legend.position="bottom",
                 plot.title=element_text( face="bold", size = unit(14, "pt")),
                 strip.text.x = element_text(size=unit(10, "pt"), colour="black"),
                 axis.title.x = element_text(size = unit(10, "pt"),face="bold", colour="black",vjust=-0.5),
                 axis.text.x=element_text(size = unit(10, "pt"), colour="black"),
                 strip.text.y = element_text(size=unit(10, "pt"), face="bold", colour="black"),
                 axis.title.y = element_text(size = unit(10, "pt"),face="bold", colour="black", vjust=1.5),
                 axis.text.y=element_text( size = unit(10, "pt"),hjust=0, colour="black"),
                 strip.background = element_rect( colour="#808080", fill="#DBDBDB", size=1, linetype="solid"),
                 legend.background = element_rect(colour = "black"),
                 legend.title = element_text(size=unit(8, "pt"), face="bold"),
                 legend.text = element_text(size = unit(8, "pt"), face = "bold",vjust=1.5),
                 legend.key.size = unit(1.5, "cm"))+
               scale_color_manual(name="Test interpretation: ",
                                  values=colors,
                                  breaks=c("Null difference in proportions  ","Non-significant (rawp)  ",  "Significant (rawp)  ",
                                           "Significant (All_p)  ","Significant (Min_2p)  ","Significant (Only_BH)  "))+
               scale_shape_manual(name="Test interpretation: ",values=c(19,19,19,19,19,19),
                                  breaks=c("Null difference in proportions  ","Non-significant (rawp)  ",  "Significant (rawp)  ",
                                           "Significant (All_p)  ", "Significant (Min_2p)  ", "Significant (Only_BH)  ")))
  Datad <- data[which(data$Allele_Type=="D allele"),]
  pd <- with(Datad,ggplot(Datad,aes(x=factor(Allele_Name, levels=rev(factor(Allele_Name))),
                                    y =Difference_proportion))+
               geom_errorbar(aes(col = Test_interpretation, ymin =Lower_bound_IC_diff_prop,
                                 ymax =Upper_bound_IC_diff_prop), cex=0.75, width=0.6)+ theme_grey()+coord_flip()+
               geom_point(aes(shape = Test_interpretation), size = 1) + ylab("Difference in proportions & 95% CI")+
               ylim( round(min(data$Lower_bound_IC_diff_prop)-0.01,2),  round(max(data$Upper_bound_IC_diff_prop)+0.01,2)) +
               geom_hline(aes(yintercept = 0)) + xlab("D alleles")+
               theme(
                 panel.grid.major = element_line(size = 0.5),
                 panel.border = element_rect(linetype = "solid", colour = "#808080", fill=NA),
                 legend.position="bottom",
                 plot.title=element_text( face="bold", size = unit(14, "pt")),
                 strip.text.x = element_text(size=unit(10, "pt"), colour="black"),
                 axis.title.x = element_text(size = unit(10, "pt"),face="bold", colour="black",vjust=-0.5),
                 axis.text.x=element_text(size = unit(10, "pt"), colour="black"),
                 strip.text.y = element_text(size=unit(10, "pt"), face="bold", colour="black"),
                 axis.title.y = element_text(size = unit(10, "pt"),face="bold", colour="black", vjust=1.5),
                 axis.text.y=element_text( size = unit(10, "pt"),hjust=0, colour="black"),
                 strip.background = element_rect( colour="#808080", fill="#DBDBDB", size=1, linetype="solid"),
                 legend.background = element_rect(colour = "black"),
                 legend.title = element_text(size=unit(8, "pt"), face="bold"),
                 legend.text = element_text(size = unit(8, "pt"), face = "bold",vjust=1.5),
                 legend.key.size = unit(1.5, "cm"))+
               scale_color_manual(name="Test interpretation: ",
                                  values=colors,
                                  breaks=c("Null difference in proportions  ","Non-significant (rawp)  ",  "Significant (rawp)  ",
                                           "Significant (All_p)  ","Significant (Min_2p)  ","Significant (Only_BH)  "))+
               scale_shape_manual(name="Test interpretation: ",values=c(19,19,19,19,19,19),
                                  breaks=c("Null difference in proportions  ","Non-significant (rawp)  ",  "Significant (rawp)  ",
                                           "Significant (All_p)  ", "Significant (Min_2p)  ", "Significant (Only_BH)  ")))
  list(VJgenes=pbrut,Vgenes=pv,Dgenes=pd,Jgenes=pj)}


