#1  packages installation

install.packages('vcfR')
library(vcfR)
install.packages('ggplot2')
library(ggplot2)
install.packages('stats')
library(stats)

#2  compute PCA 
pca <- PCA(read.vcfR('SV.vcf') , read.csv('samples.csv'))

#3  PC1 and PC2 plotting 

pca1_loading = pca$x
PC1 <- pca1_loading[, 1]
PC2 <- pca1_loading[, 2]
samples <- read.csv('samples.csv')
Population <- samples$Population
x <- data.frame(PC1, PC2, Population)
a <- c(1, 1, 5, 3, 4, 3, 3, 2, 1, 4, 4, 5, 1, 4, 5, 3, 3, 1, 1, 2, 2, 5, 2, 5, 4, 1) + 1
ggplot(x, aes(x = PC1, y = PC2)) +  geom_point(aes(shape = Population, color = Population)) + scale_shape_manual(values = c(1:10, 1:10, 1:6))  + scale_color_manual(values = a) +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold")) +
  theme(
    legend.title = element_text(color = "black", size = 14),
    legend.text = element_text(color = "black", size = 20)
  ) 




#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 PCA: PCA computation
#              input- vcf_vcfR...vcf file 
#                     samples...file containing samples with annotated population, continental group and gender
#             output- pca(prcomp object):  sdev (standard deviations of the PCs), 
#                     rotation (matrixwhose columns contain the eigenvectors ) and 
#                     x (the coordinates of the samples on the PCs)
PCA <- function(vcf_vcfR,samples) {
       info <- INFO2df(vcf_vcfR)
       GT <- vcf_vcfR@gt
       vcf_vcfR<-c()
       SV_type <- info$SVTYPE
       DEL <- GT[grepl('DEL', SV_type), 2:2505]
       SV_type<-c()
       info<-c()
       GT<-c()
       m <- transform(DEL)
       DEL<-c()
       colnames(m) <- samples$Population
       rowmeans <- rowMeans(m, na.rm = T)
       p <- sqrt((rowmeans / 2) * (1 - rowmeans / 2))
       for (i in 1:2504) {
           m[, i] <- m[, i] - rowmeans
           m[, i] <- m[, i] / p
           }
       m <- na.omit(m)
       pca <- prcomp(t(m), rank. = 2,center = T)
       return(pca)
       }



