#1  packages installation

install.packages('vcfR')
library(vcfR)
install.packages('dplyr')
library(dplyr)
install.packages('tidyr')
library(tidyr)
install.packages('ggplot2')
library(ggplot2)

#2  computation of two data frames containing homozygosity and heterozygosity numbers along samples
HomoHeterozygosity(read.vcfR('VCF file goes here')) #upload VCF file

#3 plot results
homo %>% 
  gather(key="Populations", value="hom.deletions") %>%
  ggplot( aes(x=Populations, y=hom.deletions, fill=Populations)) +
  scale_x_discrete(limits=P) + geom_violin() + coord_flip() + theme(legend.position = "none") +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))
hete %>% 
  gather(key="Populations", value="het.deletions") %>%
  ggplot( aes(x=Populations, y=het.deletions, fill=Populations)) +
  scale_x_discrete(limits=P) + geom_violin() + coord_flip()+ theme(legend.position = "none")+
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))




#function HomoHeterozygosity: compute numbers of homozygous and heterozygous events along all samples and creates two separate data frames 
#                             input- vcf...vcf file containing variants 
#                            output- homo(data.frame)...table containing numbers of homozygous events along samples
#                                    hete(data.frame)...table containing numbers of heterozygous events along samples
HomoHeterozygosity<-function(vcf){
  GT<-vcf@gt #get genotype table 
  info<-INFO2df(vcf)
  SV_type<-info$SVTYPE
  samples<-read.csv('samples.csv') #ulpoad file containing samples with annotated population, continental group and gender
  DEL<-GT[which(SV_type=='DEL'), 2:2505] #get deletion genotypes
  max_length<-max(table(samples$Population))
  p<-0 #counter for column number in data frame 
  homo<-matrix(nrow = 113,ncol = 26)
  hete<-matrix(nrow = 113,ncol = 26)
  P<-rev(c('STU','PJL','ITU','GIH','BEB','TSI','IBS','GBR','FIN','CEU','KHV','JPT','CHS','CHB','CDX','PUR','PEL','MXL','CLM','ASW','ACB','YRI','MSL','LWK','GWD','ESN'))
  assign("P",P,envir = globalenv()) #upload P vector into environment
  for (j in P){
    p<-p+1
    GT_DEL_het<-c()
    GT_DEL_homo<-c()
    pop<-as.character(samples$ï..ID[which(samples$Population==j)])
    GT_pop<-DEL[,pop] #get genotypes only for a specific population j
    for (i in 1:length(pop)){
      x<-GT_pop[,i] #vector of genotypes from a single sample i from population j
      GT_DEL_homo<-c(GT_DEL_homo,length(which(x=='1|1' | x=='2|2' | x=='3|3' ))) #count number of homozygous events per sample
      GT_DEL_het<-c(GT_DEL_het,length(which(x=='0|1' | x=='0|2' | x=='0|3' | x=='0|4' | x=='1|0' | x=='1|2' | x=='2|0' | x=='2|1' | x=='3|0' | x=='3|2'  ) )) #count number of heterozygous events per sample
    } 
    
    
    pom<-c(GT_DEL_homo,rep(NA,max_length-length(GT_DEL_homo)))
    homo[,p]<-pom
    pom2<-c(GT_DEL_het,rep(NA,max_length-length(GT_DEL_het)))
    hete[,p]<-pom2
  }
  homo<-as.data.frame(homo)
  colnames(homo)<-P
  hete<-as.data.frame(hete)
  colnames(hete)<-P
  
  assign("homo",homo,envir = globalenv()) #upload data frame homo into environment
  assign("hete",hete,envir = globalenv()) #upload data frame hete into environment
  
}

