Skip to content

Commit 3ddfb24

Browse files
apallejaapalleja
committed
Final code linting of 03 and 04 markdowns
1 parent 89d9dd2 commit 3ddfb24

2 files changed

Lines changed: 56 additions & 72 deletions

File tree

‎tutorial/03_Abundance-overview_tutorial.Rmd‎

Lines changed: 29 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -77,14 +77,18 @@ training_dir <- "/workspaces/dsp_metagenomics_training/"
7777

7878
```{r}
7979
# Abundance table
80-
abundance <- readRDS(file = paste0(training_dir, "data/MetaphlanAbundance_Species.rds"))
80+
abundance <- readRDS(file = paste0(training_dir,
81+
"data/MetaphlanAbundance_Species.rds"))
8182
rownames(abundance) <- gsub("_SRR_db1.metaphlan", "", rownames(abundance))
8283
8384
# Taxonomical annotation
84-
annotation <- readRDS(file = paste0(training_dir, "data/MetaphlanAnnotations_Species.rds"))
85+
annotation <- readRDS(file = paste0(training_dir,
86+
"data/MetaphlanAnnotations_Species.rds"))
8587
8688
# Metadata
87-
metadata <- read.table(file = paste0(training_dir, "data/metadata.tsv"), header = TRUE,
89+
metadata <- read.table(file = paste0(training_dir,
90+
"data/metadata.tsv"),
91+
header = TRUE,
8892
sep = "\t",
8993
quote = "",
9094
row.names = NULL)
@@ -95,8 +99,8 @@ rownames(metadata) <- metadata$Sample
9599
## Take only samples with sequencing and metadata information
96100

97101
```{r}
98-
abundance <- abundance[rownames(abundance) %in%
99-
rownames(metadata),]
102+
abundance <- abundance[rownames(abundance) %in%
103+
rownames(metadata), ]
100104
```
101105

102106

@@ -105,8 +109,8 @@ abundance <- abundance[rownames(abundance) %in%
105109
```{r}
106110
results_dir <- paste0(training_dir, "results/report/03_Abundance-overview/")
107111
108-
dir.create(results_dir,
109-
recursive = TRUE,
112+
dir.create(results_dir,
113+
recursive = TRUE,
110114
showWarnings = FALSE)
111115
```
112116

@@ -115,14 +119,14 @@ dir.create(results_dir,
115119
```{r}
116120
taxa_total <- colSums(abundance)
117121
118-
top_taxa <- names(sort(taxa_total,
122+
top_taxa <- names(sort(taxa_total,
119123
decreasing = TRUE))[1:10]
120124
```
121125

122126
## Filter to top 10 taxa
123127

124128
```{r}
125-
abundance_top10 <- as.data.frame(abundance) %>%
129+
abundance_top10 <- as.data.frame(abundance) %>%
126130
select(all_of(top_taxa))
127131
```
128132

@@ -142,13 +146,13 @@ abundance_long <- abundance_top10 %>%
142146
values_to = "Abundance")
143147
144148
abundance_long$Species_name <- annotation$species[match(abundance_long$Species,
145-
rownames(annotation))]
149+
rownames(annotation))]
146150
147151
abundance_long$Area <- metadata$Area[match(abundance_long$Sample,
148152
metadata$Sample)]
149153
150154
abundance_long$Sampling <- metadata$Sampling[match(abundance_long$Sample,
151-
metadata$Sample)]
155+
metadata$Sample)]
152156
```
153157

154158
## Stack bar plot of the top-10 most dominant species
@@ -167,33 +171,25 @@ abundance_long$Sampling <- metadata$Sampling[match(abundance_long$Sample,
167171
```{r}
168172
abundance_long <- abundance_long %>%
169173
arrange(Sampling, Sample) %>%
170-
mutate(Sample = factor(Sample,
174+
mutate(Sample = factor(Sample,
171175
levels = unique(Sample)))
172176
173-
pbar <- ggplot(abundance_long, aes(x = Sample,
177+
pbar <- ggplot(abundance_long, aes(x = Sample,
174178
y = Abundance,
175179
fill = Species_name)) +
176-
177180
geom_bar(stat = "identity",
178181
position = "stack",
179182
show.legend = TRUE) +
180-
181183
scale_fill_brewer(palette = "Paired") +
182-
183184
labs(title = "Top-10 Most Abundant Species",
184185
x = "Sample",
185186
y = "Relative Abundance (%)") +
186-
187187
theme_bw() +
188-
189-
theme(axis.text.x = element_text(angle = 45,
188+
theme(axis.text.x = element_text(angle = 45,
190189
hjust = 1)) +
191-
192190
scale_y_continuous(expand = c(0, 0)) +
193-
194191
scale_x_discrete(expand = c(0, 0)) +
195-
196-
facet_wrap(. ~ Area,
192+
facet_wrap(. ~ Area,
197193
scales = "free")
198194
199195
pbar
@@ -223,32 +219,27 @@ abundance_long <- abundance_top10 %>%
223219
224220
abundance_summary <- abundance_long %>%
225221
group_by(Area, Sampling, Species) %>%
226-
summarise(MeanAbundance = mean(Abundance),
222+
summarise(MeanAbundance = mean(Abundance),
227223
.groups = "drop") #ungroup
228224
229-
abundance_summary$Species_name <-
225+
abundance_summary$Species_name <-
230226
annotation$species[match(abundance_summary$Species,
231227
rownames(annotation))]
232228
233-
pbar_sum <- ggplot(abundance_summary,
234-
aes(x = interaction(Area,
235-
Sampling,
229+
pbar_sum <- ggplot(abundance_summary,
230+
aes(x = interaction(Area,
231+
Sampling,
236232
sep = " | "),
237-
y = MeanAbundance,
233+
y = MeanAbundance,
238234
fill = Species_name)) +
239-
240-
geom_bar(stat = "identity",
241-
position = "stack",
235+
geom_bar(stat = "identity",
236+
position = "stack",
242237
show.legend = TRUE) +
243-
244238
scale_fill_brewer(palette = "Paired") +
245-
246239
theme_bw() +
247-
248240
labs(title = "Mean Abundance of Top 10 Taxa by Area and Sampling",
249241
x = "Area | Sampling",
250242
y = "Mean Relative Abundance (%)") +
251-
252243
theme(axis.text.x = element_text(hjust = 0.5),
253244
strip.text = element_text(face = "bold")) +
254245
scale_y_continuous(expand = c(0, 0)) +
@@ -257,8 +248,8 @@ pbar_sum <- ggplot(abundance_summary,
257248
pbar_sum
258249
259250
# Save it as a png
260-
ggsave(paste0(results_dir, "02_top10_species_per_area_sampling.png"),
261-
width = 8, height = 6, dpi = 300)
251+
ggsave(paste0(results_dir, "02_top10_species_per_area_sampling.png"),
252+
width = 8, height = 6, dpi = 300)
262253
263254
```
264255

‎tutorial/04_Group-comparison_tutorial.Rmd‎

Lines changed: 27 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -79,21 +79,25 @@ training_dir <- "/workspaces/dsp_metagenomics_training/"
7979

8080
```{r}
8181
# Abundance table
82-
abundance <- readRDS(file = paste0(training_dir, "data/MetaphlanAbundance_Species.rds"))
82+
abundance <- readRDS(file = paste0(training_dir,
83+
"data/MetaphlanAbundance_Species.rds"))
8384
rownames(abundance) <- gsub("_SRR_db1.metaphlan", "", rownames(abundance))
8485
8586
# Taxonomical annotation
86-
annotation <- readRDS(file = paste0(training_dir, "data/MetaphlanAnnotations_Species.rds"))
87+
annotation <- readRDS(file = paste0(training_dir,
88+
"data/MetaphlanAnnotations_Species.rds"))
8789
8890
# Metadata
89-
metadata <- read.table(file = paste0(training_dir, "data/metadata.tsv"), header = TRUE,
91+
metadata <- read.table(file = paste0(training_dir,
92+
"data/metadata.tsv"),
93+
header = TRUE,
9094
sep = "\t",
9195
quote = "",
9296
row.names = NULL)
9397
rownames(metadata) <- metadata$Sample
9498
9599
# Take only samples with sequencing and metadata information
96-
abundance <- abundance[rownames(metadata),]
100+
abundance <- abundance[rownames(metadata), ]
97101
```
98102

99103
## Creating a directory for the related results
@@ -187,7 +191,7 @@ kruskal_results <- lapply(filtered_taxa,
187191
)
188192
})
189193
190-
kruskal_df <- do.call(rbind,
194+
kruskal_df <- do.call(rbind,
191195
kruskal_results)
192196
rownames(kruskal_df) <- kruskal_df$Taxon
193197
```
@@ -204,7 +208,7 @@ method (False Discovery rate)
204208
kruskal_df$adj_p_value <- p.adjust(kruskal_df$p_value,
205209
method = "BH")
206210
207-
significant_taxa <- kruskal_df %>%
211+
significant_taxa <- kruskal_df %>%
208212
filter(adj_p_value < 0.05)
209213
210214
significant_taxa$Species_name <-
@@ -231,26 +235,22 @@ DT::datatable(data = significant_taxa,
231235
## Run Dunn’s test for significant taxa
232236

233237
```{r, warning=FALSE}
234-
dunn_results <- lapply(significant_taxa$Taxon,
238+
dunn_results <- lapply(significant_taxa$Taxon,
235239
function(taxon) {
236-
237-
test_data <- abundance_df[, c(taxon, "Area")]
238-
colnames(test_data) <- c("Abundance", "Area")
239-
test_data$Abundance <- as.numeric(test_data$Abundance)
240+
test_data <- abundance_df[, c(taxon, "Area")]
241+
colnames(test_data) <- c("Abundance", "Area")
242+
test_data$Abundance <- as.numeric(test_data$Abundance)
240243
241244
# Run dunnTest only if we have >1 group with non-zero data
242245
tryCatch({
243246
dunn <- dunnTest(Abundance ~ Area,
244247
data = test_data,
245248
method = "bh")
246-
247249
dunn_df <- dunn$res
248-
249250
dunn_df$Taxon <- taxon
250-
251251
return(dunn_df)
252252
}, error = function(e) {
253-
message(paste("Skipping",
253+
message(paste("Skipping",
254254
taxon, "due to error:", e$message))
255255
})
256256
})
@@ -302,15 +302,15 @@ top_taxa <- kruskal_df %>%
302302
```{r}
303303
abundance_df$Sampling <- metadata$Sampling
304304
abundance_long <- abundance_df %>%
305-
select(Sample,
306-
Area,
307-
Sampling,
305+
select(Sample,
306+
Area,
307+
Sampling,
308308
all_of(top_taxa)) %>%
309309
pivot_longer(
310310
cols = all_of(top_taxa),
311311
names_to = "Taxon",
312312
values_to = "Abundance"
313-
)
313+
)
314314
```
315315

316316
## Fix data-types and add annotation
@@ -321,30 +321,25 @@ abundance_long$Species_name <-
321321
annotation$species[match(abundance_long$Taxon,
322322
rownames(annotation))]
323323
324-
abundance_long$Species_name <- factor(abundance_long$Species_name,
325-
levels = unique(abundance_long$Species_name))
324+
abundance_long$Species_name <- factor(abundance_long$Species_name,
325+
levels = unique(abundance_long$Species_name))
326326
327327
```
328328

329329
## Plot heatmap with ggplot2
330330

331331
```{r}
332-
heatmap_top20 <- ggplot(abundance_long,
333-
aes(x = Sample,
332+
heatmap_top20 <- ggplot(abundance_long,
333+
aes(x = Sample,
334334
y = Species_name,
335-
fill = log10(Abundance+0.000001))) +
336-
337-
geom_tile(color = "white",
335+
fill = log10(Abundance + 0.000001))) +
336+
geom_tile(color = "white",
338337
size = 0.2) +
339-
340-
scale_fill_viridis_c(option = "B",
338+
scale_fill_viridis_c(option = "B",
341339
name = "Abundance\n(log10)") +
342-
343-
facet_grid(. ~ Area + Sampling,
340+
facet_grid(. ~ Area + Sampling,
344341
scales = "free") +
345-
346342
theme_bw() +
347-
348343
theme(
349344
axis.text.x = element_blank(),
350345
axis.ticks.x = element_blank(),
@@ -361,8 +356,6 @@ ggsave(paste0(results_dir, "04_significants_heatmap.png"),
361356
width = 8, height = 8, dpi = 300)
362357
```
363358

364-
365-
366359
# Printing all package versions (good practice to ensure reproducibility)
367360
```{r session-info, echo=FALSE}
368361
sessionInfo()

0 commit comments

Comments
 (0)