Maximum value for median proportion

[amusheg]

Hello, I'm not sure what the best way to leave a reproducible example is, but possibly this is just a interpretation question. I'm making a heat tree in which node colors represent the median proportion of a taxon in a treatment group. According to the legend, the range of these values is 0-0.17. The central node, for Bacteria, is colored with the "maximum" color, implying that the median proportion of bacteria in a sample is 17%. However, 100% of the reads in the data set are bacteria; so the "median proportion" value for that node should be 1. The samples are dominated by Proteobacteria, so the median proportion value for the Proteobacteria node should also be higher than 17%. Could you help me understand what's going on? Please let me know what other information I should provide.

(As a second issue, the node size is supposed to represent the number of samples that contain that taxon, but there is no range of values shown on the figure legend for that variable.)

Thank you for your help!

devtools::install_github("ropensci/taxa")
devtools::install_github("grunwaldlab/metacoder")
library(metacoder) #0.2.0.9012 # https://github.com/grunwaldlab/metacoder
obj<-parse_phyloseq(adata)

obj$data$tax_abund <- calc_taxon_abund(obj, "otu_table")
obj$data$tax_occ <- calc_n_samples(obj, "tax_abund", groups = sample_data(adata)$treatment)
obj$data$tax_props <- calc_obs_props(obj, "tax_abund")
obj$data$diff_table <- compare_groups(obj, dataset = "tax_abund", cols=obj$sam_data$sample_ids, groups = sample_data(adata)$treatment)
obj$data$medians<-calc_group_median(obj, dataset="tax_props", groups=sample_data(adata)$treatment)

#Show taxonomic trees with presence/absence and median relative abundance 
#of taxa in each treatment
#Figures S3A-C
ht_AUT<-heat_tree(obj,
                  node_size = obj$data$tax_occ$AUT,
                  node_color = obj$data$medians$AUT,
                  node_label = taxon_names,
                  node_size_axis_label = "Samples with reads",
                  node_color_axis_label = "Median proportion")

ht_NET<-heat_tree(obj,
                  node_size = obj$data$tax_occ$NET,
                  node_color = obj$data$medians$NET,
                  node_label = taxon_names,
                  node_size_axis_label = "Samples with reads",
                  node_color_axis_label = "Median proportion")
ht_SED<-heat_tree(obj,
                  node_size = obj$data$tax_occ$SED,
                  node_color = obj$data$medians$SED,
                  node_label = taxon_names,
                  node_size_axis_label = "Samples with reads",
                  node_color_axis_label = "Median proportion")

[zachary-foster]

Hello @amusheg, its hard for me to say for sure without the data to run the code.

What does range(obj$data$medians$AUT), range(obj$data$medians$NET), and range(obj$data$medians$SED) return?

Can you attach an example of the output?

If you want, you can email me the input data and I can try to reproduce the problem. Thanks!

[amusheg]

Thank you for responding! I will email you.

[amusheg]

Hi there - my dataset and code can be found here: https://github.com/amusheg/Daphnia-microbiota-behavior The cleanup and subsetting of the dataset is in lines 1-67 of the script, and then the metacoder part starts on line 385. The output I get is attached. Thank you so much for your help! I love this package but want to make sure I'm using it correctly.

…

On Fri, Oct 5, 2018 at 10:58 AM Zachary Foster ***@***.***> wrote: Hello @amusheg <https://github.com/amusheg>, its hard for me to say for sure without the data to run the code. What does range(obj$data$medians$AUT), range(obj$data$medians$NET), and range(obj$data$medians$SED) return? Can you attach an example of the output? If you want, you can email me the input data and I can try to reproduce the problem. Thanks! — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#249 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AeKjy4rIvneSH2dTVnWbbICgm4LYxrOVks5uh3OpgaJpZM4XHQYE> .

[amusheg]

And to answer the range question - it indeed gives me a range from 0 to 0.17, just like the tree legend says. Does this mean it's doing something different than I think it's doing?

[zachary-foster]

Hello @amusheg, Thanks for the code and data! I have replicated it and the lack of the size scale is a bug I think. I will look into to it and get back to you.

The medians are being displayed correctly, but the issue is that you calculated proportions on the summed taxon counts instead of the OTU counts. When the proportions where being calculated, counts are divided by column sums and what was being summed is the taxon count columns, so counts for bacteria were divided by the sum of counts for all the phyla, classes, etc of bacteria, which is not what you intended I think. Do something like this instead:

obj$data$otu_props <- calc_obs_props(obj, "otu_table")
obj$data$tax_abund <- calc_taxon_abund(obj, "otu_props")
obj$data$tax_occ <- calc_n_samples(obj, "tax_abund", groups = sample_data(adata)$treatment)
obj$data$diff_table <- compare_groups(obj, dataset = "tax_abund", cols=obj$sam_data$sample_ids, groups = sample_data(adata)$treatment)
obj$data$medians<-calc_group_median(obj, dataset="tax_abund", groups=sample_data(adata)$treatment)

Now bacteria is always 1, as it should be in this case:

> obj$data$medians
# A tibble: 365 x 4
   taxon_id     AUT      NET      SED
   <chr>      <dbl>    <dbl>    <dbl>
 1 aab      1       1        1       
 2 aac      0.852   0.945    0.813   
 3 aad      0       0        0       
 4 aae      0       0        0       
 5 aaf      0       0        0       
 6 aag      0       0        0       
 7 aah      0       0        0       
 8 aai      0       0        0       
 9 aaj      0.00352 0.000354 0.000474
10 aak      0       0        0       
# ... with 355 more rows

[amusheg]

I'm sorry, I thought I had responded to this! Thank you for the explanation. I'll look for the bug fix for the legend labels.

…

On Wed, Oct 10, 2018 at 5:23 PM Zachary Foster ***@***.***> wrote: Hello @amusheg <https://github.com/amusheg>, Thanks for the code and data! I have replicated it and the lack of the size scale is a bug I think. I will look into to it and get back to you. The medians are being displayed correctly, but the issue is that you calculated proportions on the summed taxon counts instead of the OTU counts. When the proportions where being calculated, counts are divided by column sums and what was being summed is the taxon count columns, so counts for bacteria were divided by the sum of counts for all the phyla, classes, etc of bacteria, which is not what you intended I think. Do something like this instead: obj$data$otu_props <- calc_obs_props(obj, "otu_table")obj$data$tax_abund <- calc_taxon_abund(obj, "otu_props")obj$data$tax_occ <- calc_n_samples(obj, "tax_abund", groups = sample_data(adata)$treatment)obj$data$diff_table <- compare_groups(obj, dataset = "tax_abund", cols=obj$sam_data$sample_ids, groups = sample_data(adata)$treatment)obj$data$medians<-calc_group_median(obj, dataset="tax_abund", groups=sample_data(adata)$treatment) Now bacteria is always 1, as it should be in this case: > obj$data$medians # A tibble: 365 x 4 taxon_id AUT NET SED <chr> <dbl> <dbl> <dbl> 1 aab 1 1 1 2 aac 0.852 0.945 0.813 3 aad 0 0 0 4 aae 0 0 0 5 aaf 0 0 0 6 aag 0 0 0 7 aah 0 0 0 8 aai 0 0 0 9 aaj 0.00352 0.000354 0.000474 10 aak 0 0 0 # ... with 355 more rows — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#249 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AeKjy_f1RrNDl2OL_hoDkgW0GLzQv5_Xks5ujmUxgaJpZM4XHQYE> .

[amusheg]

Hi, I'm sorry to bother you, but I'm getting ready to publish and when generating this figure with version 0.3.0.1, the legend is still missing labels for the size scale. Do you know when this might be fixed?

[zachary-foster]

Sorry for the delay. I will work on this now. I expect I will have a solution for you within the next few days.

[zachary-foster]

Hello @amusheg, it should be fixed. Try reinstalling:

devtools::install_github("grunwaldlab/metacoder")

[amusheg]

Thank you!

…

On Mon, Dec 17, 2018 at 1:18 PM Zachary Foster ***@***.***> wrote: Hello @amusheg <https://github.com/amusheg>, it should be fixed. Try reinstalling: devtools::install_github("grunwaldlab/metacoder") — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#249 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AeKjy55iGiD0c2ff1IkmFv8uVoid1faYks5u59_tgaJpZM4XHQYE> .