Heritability of skin microbiome method confusion


I have 16S data from a twin study and want to calculate the heritability of the microbiome, both individual taxa and broad (alpha, beta). I have red a few methods from similar studies but have a few questions about my approach, not sure if the QIIME forum is the right place to post but hoping somebody here would be able to give me some advice!

  1. Filtering. I think I should be filtering out features that are not present in a certain amount of my samples from my feature table. Other papers use 50%, in low biomass samples like mine (skin) is filtering out taxa not present in >30% samples acceptable? I think this helps to reduce false positives, but bacteria only present in some samples may be important to answer my research question (I am interested in whether there are heritable bacteria that are involved in human attractiveness to mosquitoes) so not sure if I should filter at all?

  2. How should I normalise my frequency table? I have my count table, can I just convert to relative frequency or should I be using DESeq2/etc to get a normalised table? The papers I have read convert to relative frequency and then rarefy but I have convinced myself against rarefaction… I don’t think converting to proportions is enough?

  3. Should I calculate heritability at every level: feature, species, genus… I think this would be possible using my phyloseq object I have in R.

  4. Heritability of beta. Is it ok to use the first three PC from the beta PCoA I get out of QIIME2. If not how to get the rarefied beta diversity for QIIME2 into R nicely?

  5. Multiple testing correction. Benjamin-Hochberg to correct for multiple testing as I will be testing lots of taxa?

Have other looked at microbiome heritability before? Any general advice on microbiome heritability would be much appreciated!!

1 Like

Hi @alicias1,

You have a lot going on! Its sounds like a super cool study. I think my favorite (and its relatively old) paper on Heritability and the microbiome is Human Genetics Shape the Gut Microbiome. However, it’s a 2014 publication and we’ve figured out more about the microbiome since. So, it might be a starting place. I’d recommend balancing it with Microbiome Datasets Are Compositional: And This Is Not Optional, the DOI of which will probably be on my wedding announcement because of how often I cite it. I would also recommend your own search of the literature because I know there are a lot of heritability papers that I’m missing.

I’m also going to recommend straight away that you look into qiime2R for all your artifact to phyloseq needs.

With that caveat, let me see if I can try and answer your questions.

I agree with this completely! But, I think it depends on two factors. First, how do you plan to model, and second, how many features do you lose? I tend to work in the gut, where I like to filter to present in at least 10% of my samples. I do this because I sometimes to run prevalence (presence/absense) based models and those are happy with a 90/10 split but throw a fit when I try to push it to 95/5. Sparsity is an important feature fo the microbiome and so I think you should keep even somewhat sparse features. At the same time, a feature present in a single individual is insuffeciently powered to allow analysis.

I’m going to refer you to Microbiome Datasets are Compositional again. DeSeq2 doesn’t perform well against compositional models, either. I don’t know if your heritability can be modeled using something like an OLS (it’s been ages since I did ACE calculations), but if it can be, then songbird might be a good option to run.

This is so some degree a personal taste thing. I like feature and maybe genus. I abhor species level descriptions in amplicon data because species-level annotation is funky and Im making the political stand here that we should abandon species for amplicon IDs. I also don’t htink you get much out of something above genus to family level. Collapsing makes the data more tractable, but my dog (:dog:), my cat (:cat:), and my theretical ferret (:house: :shark:) play different roles in an ecosystem even though they all belong to the same order.

If you’re doing heritability, I assume you’re working with twins? Could you use paired distances somehow? Maybe look at what’s possible with q2-longitudinal?

It’s what I tend to use on my models.



Wow thank you for such a detailed response! Definitely gives me some more to think about…

From reading this heritability papers (and a few others) and the compositional paper there are other approaches available in 2020 than these heritability studies used and I should try to use a compositional approach? I haven’t seen a heritability study that doesn’t rarefy… I will look harder… and try to understand this compositional approach. Is there a tutorial available for a compositional pipeline?

Ok I am still confused how to normalise my data… I will look into songbird. I have been using ASRreml to model heritability but the data I am putting in is relative frequency and this isn’t enough I don’t think.

I hadn’t thought about that, yes I have twins, I will look into that tool!

Thank you for answering so many of my questions that makes me feel a bit more confident! I didn’t anticipate there not to be a clear approach available for answering some of my questions when I started this project but am very grateful for this forum :slight_smile:


Hi @alicias1,

Ive been thinking about your problem more. And, I actually want to go back to beta diversity (paired analysis). I really encourage you to look there first, for two reasons. First, this will help you determine if any of the elements of the microbiome related to your phenotype are heritable. If there is heritability, i think you would assume that concordant MZ twins should have more similar microbiomes than concordant DZ twins. Im not sure how the discorant twin pairs fit here, but important to think about? I also like this approach (and I use it a lot because if you use multiple metrics, it helps guide your next steps. So, if you see a difference in an unweighted metric, it probably means you’re looking at changes in presence/absence. If you only see it in weighted, there’s an abundance-based difference. In each case, I think the work you need to do to approach the hypothesis in a heritability model becomes different. For example, if there is no difference in beta diversity, then perhaps you should consider that the trait isn’t heritable or isn’t heritable enough for your sample size. (Which, hey, that’s also interesting). If there’s a difference in a presence/absence based metric, suddenly your normalization problem has become a lot easier because you don’t need to normalize, just define “present”. And, if it’s abundance-based, then you need some more thought.

The more I think about this, the more I think you’re going to need a bespoke solution. You could potentially wrap ASReml in an ANCOM-like interface (but that requires writing your own R code), or you may need to solve the heritability model for the trait and convert to an OLS or something that can be passed into a model. These may be things to discuss with a local biostatistican if you have one. But, as far as I know, there’s no good compositionally aware models



Thank you for the explanation, I have had a go at pairwise distance comparisons using q2-longitudinal. Here is the weighted and unweighted UniFrac distances comparing twin pairs between MZ and DZ pairs:

I think the weighted suggests some evidence of a difference between MZ and DZ pairs, more concordance in MZ twin pairs than DZ pairs. So there may be an abundance based difference.

So, my normalisation problem still exists… I will discuss with a statistician and see what he thinks!

1 Like