<p>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: 
<a href="https://i.sstatic.net/0znDe.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/0znDe.png" alt="enter image description here"></a>
Here the missingess threshold for genotype inclusion is &lt;20%.</p>

<p>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 <code>ape::PCoA()</code> with an Euclidian distance matrix from <code>paralleldist::pardist()</code>. When I run <code>mclust::mclustBIC()</code>, 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: 
<a href="https://i.sstatic.net/gsJl4.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/gsJl4.png" alt="enter image description here"></a></p>

<p>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.</p>

<p>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.</p>

CC BY-SA 4.0
2019-03-08T11:08:16.500Z
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
CC BY-SA 4.0
2019-03-08T14:18:50.723Z
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.
M__
CC BY-SA 4.0
2019-03-09T09:51:07.827Z
@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
CC BY-SA 4.0
2019-03-11T13:55:58.410Z
@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.
M__
CC BY-SA 4.0
2019-03-11T14:03:41.597Z
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.
M__
<r><snp><clustering><radseq><genotyping>
GBS: clustering of PCoA axes with R package mclust to describe identity by missingness
# GBS: clustering of PCoA axes with R package mclust to describe identity by missingness

(Asked by: username_0 on 2019-03-01)

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:
[![enter image description here](https://i.sstatic.net/0znDe.png)](https://i.sstatic.net/0znDe.png)
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:
[![enter image description here](https://i.sstatic.net/gsJl4.png)](https://i.sstatic.net/gsJl4.png)

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.

**Tags:** r, snp, clustering, radseq, genotyping
> **Comments:**
> **username_1**: 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`  
> **username_2**: 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.  
> **username_1**: @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.  
> **username_2**: @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.  
> **username_2**: 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.  

---

## 0 answers

