# GBS: clustering of PCoA axes with R package mclust to describe identity by missingness

I have a SNP dataset with 56 populations and 430 samples that was produced by GBS sequencing. I need to identify groups of samples that have similar patterns of genotype missingness, since there is substantial structure in missing genotypes, shown here in a PCoA ordination of genotype matrix with 1-0 values, where 0 represents missing gt and 1 is presence of a gt value: Here the missingess threshold for genotype inclusion is <20%.

I want to cluster the samples using the scores on 20-30 PCoA axes, which is 17%-21% of variation, using the R package mclust. I ran ape::PCoA() with an Euclidian distance matrix from paralleldist::pardist(). When I run mclust::mclustBIC(), I get a weird result in that the there is no modal value for any of the better-supported models, here shown in the legend:

I had expected to find a model that would support the existence of a small number of clusters (2-6 maybe). I also don't understand why the VVE model only has values to 3 clusters, where it presents the maximum value. Any ideas? I have tried relaxing the minimum missing level, and using data with differing minimum locus depth levels. I get very similar results.

Could someone please help me understand what is going on here and how I might address it? Once I have clusters, I can look at the assignment of samples to clusters by population, flowcell, and other stratifying variables. I would be happy to edit to include more information, if that would be helpful.

• So, what do you want help with: Why there isn't a single model with a definite number of models? Why VVE model has only values for 3 components? or How to identify subpopulations? How many SNPS do you have? BTW it would help to visualize better the PCoAs if you use scale = "free" in facet_grid – llrs Mar 8 '19 at 11:08
• Hey @llrs following our earlier discussion ML could actually work here, like a random forest. Here, the model assumptions don't conform to the data, so rather than figure out the "magically distribution" you just wrap ML around it. There is a caveat (PCA isn't a nice ML approach in this case because it would be hard to interpret the result), but it would 'solve' it. – Michael Mar 8 '19 at 14:18
• @MichaelG. I'm not sure that model assumptions don't fit on this data, but that the result is not the one expected which is a completely different thing in my opinion. I don't understand why there should be several subpopulation in first place. The random forest has its own assumptions that should be checked too, I don't understand how using it here would help here. – llrs Mar 9 '19 at 9:51
• @llhs that is certainly a possibility, I did think about it, but either way the BIC must reach a plateau. Random forest, if used with a PCA, doesn't mean very much. If you have an a priori cluster and you want to know (in this case) what SNPs contribute to the sub-population, a random forest is one approach that will work (not using PCA). This analysis per se is interesting and if it worked, could useful. Honestly I would I simply use a traditional approach in pop gen first. PCA is when stuff doesn't work, but only in pop gen. – Michael Mar 11 '19 at 13:55
• The final thing to add is "missing genotypes,", if this means that the experiment didn't work for that particular allele, its no biggy. If it means there was an expected genotype and another genotype was found in its place, then that is a big deal. – Michael Mar 11 '19 at 14:03