xiangpin
commented
1 year ago I don't think Manhattan plot is a good idea to visualize the result of DAA, because too much information will be missing. The MicrobiotaProcess provided mp_plot_diff_boxplot, mp_plot_diff_cladogram, mp_plot_diff_res functions to visualize the results of DAA directly. But you can also use mp_plot_abundance to visualize the result. You can refer to the following codes.
> library(MicrobiotaProcess)
> data(mouse.time.mpse)
> # You can obtain the result of DAA and extract it with taxonomy information
> mouse.time.mpse %>% mp_rrarefy() %>% mp_diff_analysis(.abundance=RareAbundance, .group=time) %>% mp_extract_taxatree() -> taxa.tree
> taxa.tree
'treedata' S4 object'.
...@ phylo:
Phylogenetic tree with 218 tips and 186 internal nodes.
Tip labels:
OTU_67, OTU_231, OTU_188, OTU_150, OTU_207, OTU_5, ...
Node labels:
r__root, k__Bacteria, p__Actinobacteria, p__Bacteroidetes, p__Cyanobacteria,
p__Deinococcus-Thermus, ...
Rooted; no branch lengths.
with the following features available:
'nodeClass', 'nodeDepth', 'RareAbundanceBySample', 'LDAupper', 'LDAmean',
'LDAlower', 'Sign_time', 'pvalue', 'fdr'.
# The associated data tibble abstraction: 404 × 12
# The 'node', 'label' and 'isTip' are from the phylo tree.
node label isTip nodeCl…¹ nodeD…² RareAb…³ LDAup…⁴ LDAmean LDAlo…⁵ Sign_…⁶
<int> <chr> <lgl> <chr> <dbl> <list> <dbl> <dbl> <dbl> <chr>
1 1 OTU_67 TRUE OTU 8 <tibble> 3.39 3.36 3.32 Late
2 2 OTU_231 TRUE OTU 8 <tibble> NA NA NA NA
3 3 OTU_188 TRUE OTU 8 <tibble> NA NA NA NA
4 4 OTU_150 TRUE OTU 8 <tibble> NA NA NA NA
5 5 OTU_207 TRUE OTU 8 <tibble> NA NA NA NA
6 6 OTU_5 TRUE OTU 8 <tibble> NA NA NA NA
7 7 OTU_1 TRUE OTU 8 <tibble> NA NA NA NA
8 8 OTU_2 TRUE OTU 8 <tibble> NA NA NA NA
9 9 OTU_3 TRUE OTU 8 <tibble> NA NA NA NA
10 10 OTU_4 TRUE OTU 8 <tibble> 4.40 4.38 4.36 Late
# … with 394 more rows, 2 more variables: pvalue <dbl>, fdr <dbl>, and
# abbreviated variable names ¹nodeClass, ²nodeDepth, ³RareAbundanceBySample,
# ⁴LDAupper, ⁵LDAlower, ⁶Sign_time
# ℹ Use `print(n = ...)` to see more rows, and `colnames()` to see all variable names
>
>
> # Then you can extract the names of differential taxa
> taxa.tree %>% dplyr::filter(nodeClass=='Genus' & !is.na(Sign_time),keep.td=FALSE) %>% dplyr::pull(label)
[1] "g__Bifidobacterium" "g__un_f__Muribaculaceae"
[3] "g__Alistipes" "g__Clostridium_sensu_stricto_1"
[5] "g__ASF356" "g__Lachnospiraceae_NK4A136_group"
[7] "g__Lachnospiraceae_UCG-001" "g__Anaerotruncus"
[9] "g__Oscillibacter" "g__Ruminiclostridium"
[11] "g__Turicibacter" "g__un_f__Mitochondria"
[13] "g__Pseudomonas" "g__Anaeroplasma"
>
>
> # Then, calculate the relative abundance if you want to visualize it and filter the genus.
> mouse.time.mpse %>% mp_rrarefy() %>% mp_cal_abundance(.abundance=RareAbundance) %>% dplyr::filter(Genus %in% keep.genus) %>% mp_plot_abundance(.abundance=RelRareAbundanceBySample, force=T, relative=F, taxa.class=Genus, geom='heatmap', topn='all', .group=time) -> p1
> p1[[1]] <- p1[[1]] + scale_fill_viridis_c(option='H')
> p1[[2]] <- p1[[2]] + scale_fill_manual(values=c('orange', 'deepskyblue'))
> p1
If you want to integrate the P-value or FDR or other statistics results.
> taxa.tree %>% dplyr::filter(nodeClass=='Genus' & !is.na(Sign_time),keep.td=FALSE) %>% ggplot(aes(x=-log10(fdr), y=label, fill=Sign_time)) + geom_col() + ylab(NULL) + theme_bw() + theme(axis.text.y=element_blank()) + scale_x_continuous(expand=c(0,0,0, .1)) + scale_fill_manual(values=c('orange', 'deepskyblue')) -> p2
>
> p1 %>% aplot::insert_right(p2, width=.2)
In addition, since the output of the result is tidy-like, you also can use ggplot2 or ggplot2-extension to visualize it. For example. Manhattan plot
> # First, obtain the taxonomy information of OTU.
> taxa.da <- mouse.time.mpse %>% mp_extract_taxonomy()
> # Then, add to the `data.frame` of the significant OTUs and visualize it.
> taxa.tree %>% dplyr::filter(nodeClass=='OTU', keep.td=F) %>% dplyr::left_join(taxa.da, by=c('label'='OTU')) %>% dplyr::mutate(Sign_time=dplyr::case_when(is.na(Sign_time) ~ 'NoSign', TRUE ~ Sign_time)) %>% ggplot(aes(x=Phylum, y=-log10(fdr), color=Phylum, shape=Sign_time)) + geom_jitter(size=3) + geom_hline(yintercept=-log10(0.05)) + theme(axis.text.x=element_text(angle=-45, hjust=0))
Junhao-Guo
Can the results be presented in other ways? such like a Manhattan plot?