-
Notifications
You must be signed in to change notification settings - Fork 28
/
Copy pathissue_141--phyloseq_converter.Rmd
115 lines (87 loc) · 3.6 KB
/
issue_141--phyloseq_converter.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
### old way
```{r, eval = FALSE}
source('http://bioconductor.org/biocLite.R')
biocLite('phyloseq')
```
```{r}
phyloseq_to_taxmap <- function(phyloseq_obj) {
tax_data <- data.frame(
taxonomy = vapply(seq_len(nrow(phyloseq_obj@tax_table)), FUN.VALUE = character(1),
function(i) {
taxon_names <- as.vector(phyloseq_obj@tax_table[i, ])
rank_names <- colnames(phyloseq_obj@tax_table)
rank_names <- rank_names[!is.na(taxon_names)]
taxon_names <- taxon_names[!is.na(taxon_names)]
paste(rank_names, taxon_names, sep = "::", collapse = ";;")
}),
phyloseq_id = rownames(phyloseq_obj@tax_table)
)
otu_data <- as.data.frame(t(phyloseq_obj@otu_table))
otu_data$phyloseq_id <- rownames(otu_data)
tax_data <- merge(tax_data, otu_data, by.x = "phyloseq_id", by.y = "phyloseq_id")
parse_taxonomy_table(tax_data, taxon_col = c("class" = 2), other_col_type = "obs_info",
class_regex = "^(.*)::(.*)$",
class_key = c(rank = "taxon_info", "name"),
class_sep = ";;")
}
```
```{r}
library(phyloseq)
phy <- readRDS("example.phyloseq")
```
```{r}
data <- phyloseq_to_taxmap(phy)
```
## New way with taxa taxmap
```{r}
# devtools::install_github("grunwaldlab/metacoder")
# devtools::install_github("ropensci/taxa")
library(metacoder)
library(taxa)
phy <- readRDS("example.phyloseq")
obj <- parse_phyloseq(phy)
obj$filter_taxa(taxon_names == "Bacteria", subtaxa = TRUE)
healthy_samples <- obj$data$sam_data$sample_id[obj$data$sam_data$status == "healthy"]
healthy_counts <- obs_apply(obj, "otu_table", function(i) sum(obj$data$otu_table[i, healthy_samples]), simplify = TRUE)
obj %>%
heat_tree(node_size = healthy_counts,
node_color = healthy_counts,
node_label = taxon_names)
```
## Even newer way
```{r}
devtools::install_github("grunwaldlab/metacoder")
devtools::install_github("ropensci/taxa")
library(metacoder)
phy <- readRDS("example.phyloseq")
obj <- parse_phyloseq(phy)
# Convert counts to proportions
obj$data$otu_table <- calc_obs_props(obj,
dataset = "otu_table",
cols = obj$data$sam_data$sample_ids,
other_cols = TRUE)
# Calculate per-taxon proportions
obj$data$tax_table <- calc_taxon_abund(obj,
dataset = "otu_table",
cols = obj$data$sam_data$sample_ids)
# Calculate difference between treatments
obj$data$diff_table <- compare_treatments(obj,
dataset = "tax_table",
sample_ids = obj$data$sam_data$sample_ids,
treatments = obj$data$sam_data$status,
other_cols = TRUE)
# Plot differneces
color_interval <- c(-5, 5) # The range of values (log 2 ratio of median proportion) to display
obj %>%
taxa::filter_taxa(taxon_names == "Bacteria", subtaxa = TRUE) %>%
heat_tree(node_size_axis_label = "Number of OTUs",
node_size = n_obs,
node_color_axis_label = "Log 2 ratio of median proportions",
node_color = log2_median_ratio,
node_color_range = diverging_palette(),
node_color_trans = "linear",
node_color_interval = color_interval,
edge_color_interval = color_interval,
node_label = taxon_names,
node_label_max = 150)
```