#1  packeges installation

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


#2  set global variables 
vcf_vcfR<-read.vcfR('VCF file goes here') #uploading the vcf file 
samples<-read.csv('samples.csv')
ID<-getID(vcf_vcfR) #get SV IDs
info<-INFO2df(vcf_vcfR)  #extract data frame containing INFO columns from VCF file 
SV_type<-info$SVTYPE
GT<-vcf_vcfR@gt #get genotypes table from VCF file 
vcf_vcfR<-c()  

#3  Get Vst tables for DELETIONS, DUPLICATIONS and mCNV (!annotated as CNV!);
#Vst is estimated for each variant type individualy and a data frame is created  with function Vst_df 

DEL_Vst<-var_type_vst('DEL')
DUP_Vst<-var_type_vst('DUP')
CNV_Vst<-var_type_vst('CNV')

#4  create a data frame to compare total number of SVs with number of SVs with at least one Vst value above 0.2 
df<-Vst_df(DEL_Vst,DUP_Vst,CNV_Vst)




#function Vst_df: creates a data frame containing total number of variants in first row, and numbers of variants, with at least 
#one Vst value bigger than 0.2, in second row; used for data plotting 
#               input- DEL_Vst(data.frame)...data frame containing Vst values for each deletion site 
#                      DUP_Vst(data.frame)...data frame containing Vst values for each duplication site
#                      CNV_Vst(data.frame)...data frame containing Vst values for each mCNV site
#              output- df(data.frame)...data frame;1.row=total number of SVs (DEL,DUP,CNV);
#                                       2.row=number of SVs with at least one Vst value above 0.2 
Vst_df<-function(DEL_Vst,DUP_Vst,CNV_Vst){
        DUP_high_Vst<-sum(((rowSums(DUP_Vst[, 2:11] >= 0.2, na.rm=T) > 0) * 1)) 
        CNV_high_Vst<-sum(((rowSums(CNV_Vst[, 2:11] >= 0.2, na.rm=T) > 0) * 1))
        DEL_high_Vst<-sum(((rowSums(DEL_Vst[, 2:11] >= 0.2, na.rm=T) > 0) * 1))
   
        value<-c(dim(DEL_Vst)[1],DEL_high_Vst,dim(DUP_Vst)[1],DUP_high_Vst,dim(CNV_Vst)[1],CNV_high_Vst)
        type<-c('DEL','DEL' ,'DUP','DUP','CNV','CNV')
        condition<-rep(c('total number of SV','Vst>=0.2'),3)
        Vst<-c('number of structural variants','Vst>=0.2')
        df<-data.frame(value,type,condition)
        return(df)
        }

#function transform: transform genotype coding into 0 (homozygous reference),1 (heterozygous alternate),2 (homozygous alternate)
#                   input- x(matrix)...genotypes table
#                  output- m(matrix)...transformed genotypes 
transform<-function(x){
           x[which(x=='0')]<-0
           x[which(x=='1')]<-1
           x[which(x=='2')]<-2
           x[which(x=='0|0')]<-0
           x[which(x=='1|1')]<-2
           x[which(x=='2|2')]<-2
           x[which(x=='3|3')]<-2
           x[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' )]<-1
           m<-as.numeric(x)
           m<-matrix(data=m,ncol=dim(x)[2], nrow=dim(x)[1])
           return(m)
           }

#function var_type_vst: compute Vst values for specific variant type by calling Vst function 
#                       input- x(character)...variant type
#                      output- Vst_pop(data.frame)...table containing Vst values for each variant site (in row) and among all pairwise populations (in column)
var_type_vst<-function(x){
              Vst_pop<-data.frame(ID[grepl(x,SV_type )]) #initiate a data frame 
              pop_pairs<-combn(x = c('AFR','AMR','EUR','ASN','SAN'),2)
              for (i in 1:dim(pop_pairs)[2]){
                   a<-pop_pairs[1,i]
                   b<-pop_pairs[2,i]
                   pom<-Vst(a,b,x) #compute Vst values for each variant site between any 2 populations
                   Vst_pop<-data.frame(Vst_pop,pom)
                   }
              colnames(Vst_pop)<-c('ID','AFR-AMR','AFR-EUR','AFR-EAS','AFR-SAS','AMR-EUR','AMR-EAS','AMR-SAS','EUR-EAS','EUR-SAS','EAS-SAS')
              return(Vst_pop)
              }

#function Vst: calculate Vst values for each SV between 2 populations a and b
#              input- a(character)...population A name
#                     b(character)...population B name
#                     type(character)...SV type 
#             output- Vst(vector)...vector containing Vst values 
Vst<-function(a,b,type){
     Vst<-c()
     if(type=='DEL') {
       m<-GT[grepl('DEL',SV_type ), 2:2505]
       } 
     else if (type=='DUP'){
             m<-GT[which(SV_type=='DUP' ), 2:2505]
             } 
     else if (type=='CNV'){
             m<-GT[which(SV_type=='CNV' ), 2:2505]
            }
     a<-m[, which(samples$Continental.group==a)] 
     a<-transform(a)
     b<-m[, which(samples$Continental.group==b)]
     b<-transform(b)
     for (i in 1:dim(m)[1]){
         Va<-sum((a[i,] - mean(a[i,],na.rm = T))^2,na.rm = T) / dim(a)[2]
         Vb<-sum((b[i,] - mean(b[i,],na.rm = T))^2,na.rm = T) / dim(b)[2]
         total<-c((a[i,]),(b[i,]))
         Vt<-sum((total - mean(total,na.rm = T))^2,na.rm = T) / length(total)
         Vs<-(Va*dim(a)[2] + Vb*dim(b)[2])/length(total)
         Vst[i]<-(Vt-Vs)/Vt
         }
     return(Vst)
     }

