Downstream statisical tests

[lauramason326]

Hi
Loving the tool. Quick Q (or thought experiment). I'd like to use singleM taxonomic profiles in traditional statistical tests - ie PERMANOVA, NMDS plots, Maaslin, and correlations between enriched organisms and soil chem etc. Is it permissible to use the raw coverage values in the taxonomic profile, or would it be better to transform the data somehow? For example, if I were using 16S data for similar purposes, I would probably either rarefy the data or convert it to relative abundances.
Thanks!
Laura

[wwood]

Thanks for the kind words.

You can convert a condensed profile into one with relative abundance etc using

singlem summarise --input-taxonomic-profile doco_example.profile \
    --output-taxonomic-profile-with-extras doco_example.with_extras.tsv

Note that you should only consider one taxonomic level at a time.

How best to approach this probably depends on your dataset I think - for soil maybe comparison at the genus level might be appropriate?

I'm afraid we haven't quite settled on a precise recommendation for this. For instance, you can rarefy an OTU table with summarise as well (not a condensed profile though), but then you can only consider one marker gene at a time. It might be a good idea to do it with multiple markers and hopefully observe similar patterns.

ben

[lauramason326]

Hi Ben - thanks! I ran a couple of tests:

singlem summarise --input-taxonomic-profiles CDA_0524_CopeCD_23096_taxonomic.profile.tsv CDA_0524_CopeCD_23177_taxonomic.profile.tsv CDA_0524_CopeCD_23279_taxonomic.profile.tsv \
--output-taxonomic-profile-with-extras test_profile_with_extras.tsv

this seems to give similar output to Kraken/Braken - ie 100% of the read are classified as "root"; 95% are bacteria + 5% are archaea, etc. this works for OTUs defined to the species level, but not for OTUs identified to the family level only etc, and I am also interested in those

I also compared the above results with this:

singlem summarise --input-taxonomic-profiles CDA_0524_CopeCD_23096_taxonomic.profile.tsv CDA_0524_CopeCD_23177_taxonomic.profile.tsv CDA_0524_CopeCD_23279_taxonomic.profile.tsv \
    --output-species-by-site-relative-abundance test_otu_profile_by_species.tsv \
    --output-species-by-site-level species

This doesn't really work for my needs either - esp at species level as 98% of reads cannot be ID'd to species level

Trying this as well

cat CDA_0524_CopeCD_23096_taxonomic.profile.tsv CDA_0524_CopeCD_23177_taxonomic.profile.tsv CDA_0524_CopeCD_23279_taxonomic.profile.tsv > test_concat_profiles.tsv

head test_concat_profiles.tsv

sample coverage taxonomy
CDA_0524_CopeCD_23096_R1 2.69 Root; d__Archaea
CDA_0524_CopeCD_23096_R1 25.57 Root; d__Bacteria
CDA_0524_CopeCD_23096_R1 6.63 Root; d__Archaea; p__Thermoproteota
CDA_0524_CopeCD_23096_R1 1.49 Root; d__Archaea; p__Thermoplasmatota
CDA_0524_CopeCD_23096_R1 0.6 Root; d__Archaea; p__Halobacteriota
CDA_0524_CopeCD_23096_R1 0.48 Root; d__Archaea; p__Asgardarchaeota
CDA_0524_CopeCD_23096_R1 1.5 Root; d__Bacteria; p__Chloroflexota
CDA_0524_CopeCD_23096_R1 6.32 Root; d__Bacteria; p__Actinomycetota
CDA_0524_CopeCD_23096_R1 4.02 Root; d__Bacteria; p__Planctomycetota

I prefer this option because it incorporates the OTUs not assigned to a species level and includes their coverage in a "pivot longer" format for tidyverse in R. my question is then, how best to use these coverages to represent OTU abundance by sample? In the past, I have used these coverage values in a simple relative abundance calculation: coverage OTU/ sum coverages in sample = "relative abundance" in R, but is this appropriate?

Thanks,
Laura

[wwood]

Hi Laura,

In short:

In the past, I have used these coverage values in a simple relative abundance calculation: coverage OTU/ sum coverages in sample = "relative abundance" in R, but is this appropriate?

Yes, that's the way it is calculated within the singlem code as well.

But some caveats

this seems to give similar output to Kraken/Braken - ie 100% of the read are classified as "root"; 95% are bacteria + 5% are archaea, etc. this works for OTUs defined to the species level, but not for OTUs identified to the family level only etc, and I am also interested in those

For a start, the profile generated is are not made up of "OTUs" - see https://wwood.github.io/singlem/Glossary - I think you might be conflating the two?

If you want family-level, you can use family instead of species here i.e.

singlem summarise --input-taxonomic-profiles CDA_0524_CopeCD_23096_taxonomic.profile.tsv CDA_0524_CopeCD_23177_taxonomic.profile.tsv CDA_0524_CopeCD_23279_taxonomic.profile.tsv \
    --output-species-by-site-relative-abundance test_otu_profile_by_species.tsv \
    --output-species-by-site-level family

Then the output tables will be "sample by family".

HTH