Hello!

I am trying to understand how SNPRelate operates under the hood when samples have missing values.

Specifically, in my VCF I have 150 samples, split into 6 groups, 25 samples each (for each group, 10 samples were sequenced at 30x and 15 at 5x). When I conduct PCA (`snpgdsPCA`

), I see samples cluster according to their groups, as follows:

However, there is also a substructure within groups according to the sequence coverage (and therefore the missing rate).

This gives the PCA plot a star-like shape, where the 30x-samples form a blob towards the margins of the plot and the 5x-samples form a line towards the center of plot.

There is one group where this blob-line pattern does not occur (black-group). Samples in this group have much less missingness.

I try to replicate the PCA using the `prcomp`

function in `R`

, using the genotype count matrix I extract as follows:

```
seqVCF2GDS(vcf, "vcf.gds") #convert VCF to gds
genofile <- seqOpen("vcf.gds") # open gds
geno_cts <- read.gdsn(index.gdsn(genofile, "genotype/data")) # get allele counts (separated by chromosomes)
# Since for each SNP, there are two vectors with allele counts, add them up to get:
# homref = 0, het = 1, homalt = 2
# missing (.) encoded as 3, therefore, 6 = ./.
gt_cts_ma <- matrix(nrow = max(dim(geno_cts)), ncol = 150)
for(i in 1:max(dim(geno_cts))){
gt_cts_ma[i,] <- geno_cts[,,i] %>% apply(2, sum)
}
colnames(gt_cts_ma) <- read.gdsn(index.gdsn(genofile, "sample.id"))
gt_cts_ma_pca <- prcomp(t(gt_cts_ma))
plot(gt_cts_ma_pca$x[,"PC1"], gt_cts_ma_pca$x[,"PC2"])
```

...but the plot does not resemble the SNPRelate PCA plot at all.

Looking at the function `snpgdsPCA`

, the function `gnrPCA`

is invoked, but I cannot find what it does anywhere.

Could you please tell me where I could learn more about the detailed steps used by `snpgdsPCA`

in order to figure out what is behind the pattern I observe in my data?