其他
R语言微生物群落数据分析:从原始数据到结果（No2）
library("phyloseq")
library("ggplot2")
library("gridExtra")
ps0 = readRDS("ps.rds")
ps0ASV表格过滤
##ASV表格过滤
# Define prevalence of each taxa
# (in how many samples did each taxa appear at least once)
prev0 = apply(X = otu_table(ps0),
MARGIN = ifelse(taxa_are_rows(ps0), yes = 1, no = 2),
FUN = function(x){sum(x > 0)})
prevdf = data.frame(Prevalence = prev0,
TotalAbundance = taxa_sums(ps0),
tax_table(ps0))
keepPhyla = table(prevdf$Phylum)[(table(prevdf$Phylum) > 5)]
prevdf1 = subset(prevdf, Phylum %in% names(keepPhyla))
# Define prevalence threshold as 5% of total samples
prevalenceThreshold = 0.05 * nsamples(ps0)
prevalenceThreshold
# [1] 18
# Execute prevalence filter, using `prune_taxa()` function
ps1 = prune_taxa((prev0 > prevalenceThreshold), ps0)
ps1
# phyloseq-class experiment-level object
# otu_table() OTU Table: [ 353 taxa and 360 samples ]
# sample_data() Sample Data: [ 360 samples by 14 sample variables ]
# tax_table() Taxonomy Table: [ 353 taxa by 6 taxonomic ranks ]
# phy_tree() Phylogenetic Tree: [ 353 tips and 351 internal nodes ]
# Filter entries with unidentified Phylum.
ps2 = subset_taxa(ps1, Phylum %in% names(keepPhyla))
ps2
ggplot(prevdf1, aes(TotalAbundance, Prevalence, color = phylum)) +
geom_hline(yintercept = prevalenceThreshold, alpha = 0.5, linetype = 2) +
geom_point(size = 2, alpha = 0.7) +
scale_y_log10() + scale_x_log10() +
xlab("Total Abundance") +
facet_wrap(~phylum)
# phyloseq-class experiment-level object
# otu_table() OTU Table: [ 341 taxa and 360 samples ]
# sample_data() Sample Data: [ 360 samples by 14 sample variables ]
# tax_table() Taxonomy Table: [ 341 taxa by 6 taxonomic ranks ]
# phy_tree() Phylogenetic Tree: [ 341 tips and 339 internal nodes ]
提取属等级
taxGlomRank = "genus"
length(get_taxa_unique(ps2, taxonomic.rank = taxGlomRank))
#表明一共有47个属
## [1] 47
ps3 = tax_glom(ps2, taxrank = taxGlomRank)
# Phylogenetic agglomeration
# How many genera are present after filtering?
h1 = 0.4
ps4 = tip_glom(ps2, h = h1)
可视化上面操作的进化树
multiPlotTitleTextSize = 8
p2tree = plot_tree(ps2,
method = "treeonly",
ladderize = "left",
title = "Before Agglomeration") +
theme(plot.title = element_text(size = multiPlotTitleTextSize))
p3tree = plot_tree(ps3,
method = "treeonly",
ladderize = "left",
title = paste0("tax_glom(..., taxrank='", taxGlomRank, "')")) +
theme(plot.title = element_text(size = multiPlotTitleTextSize))
p4tree = plot_tree(ps4,
method = "treeonly",
ladderize = "left",
title = paste0("tip_glom(..., h=", h1, ")")) +
theme(plot.title = element_text(size = multiPlotTitleTextSize))
# group plots together
grid.arrange(nrow = 1, p2tree, p3tree, p4tree)
##进行count丰度和相对丰度的比较
plot_abundance = function(physeq, ylabn = "",
Facet = "order",
Color = "phylum"){
mphyseq = psmelt(physeq)
mphyseq <- subset(mphyseq, Abundance > 0)
ggplot(data = mphyseq,
mapping = aes_string(x = "sex", y = "Abundance",
color = Color, fill = Color)) +
geom_violin(fill = NA) +
geom_point(size = 1, alpha = 0.3,
position = position_jitter(width = 0.3)) +
facet_wrap(facets = Facet) + ylab(ylabn) +
scale_y_log10()
}
# Transform to relative abundance. Save as new object.
ps3ra = transform_sample_counts(ps3, function(x){x / sum(x)})
plot_abundance(ps3,"Original Abundances")
plot_abundance(ps3ra,"Relative Abundances")
挑选指定的门类去展示
#挑选指定的分类等级去展示
psOrd = subset_taxa(ps3ra, order == "Lactobacillales")
plot_abundance(psOrd, Facet = "genus", Color = NULL)
qplot(sample_data(ps0)$age, geom = "histogram") + xlab("age")
qplot(log10(rowSums(otu_table(ps0)))) +
xlab("Logged counts-per-sample")
#对小鼠年龄进行统计
qplot(sample_data(ps0)$age, geom = "histogram") + xlab("age")
#统计样品的序列条数
qplot(log10(rowSums(otu_table(ps0)))) +
xlab("Logged counts-per-sample")
pslog <- transform_sample_counts(ps0, function(x) log(1 + x))
sample_data(pslog)$age_binned <- cut(sample_data(pslog)$age,
breaks = c(0, 100, 200, 400))
out.wuf.log <- ordinate(pslog, method = "MDS", distance = "unifrac")
evals <- out.wuf.log$values$Eigenvalues
plot_ordination(pslog, out.wuf.log, color = "age_binned") +
labs(col = "Binned Age") +
coord_fixed(sqrt(evals[2] / evals[1]))
setup_example(c("phyloseq", "ggplot2", "plyr", "dplyr", "reshape2",
"ade4", "ggrepel"))
out.bc.log <- ordinate(pslog, method = "MDS", distance = "bray")
evals <- out.bc.log$values$Eigenvalues
plot_ordination(pslog, out.bc.log, color = "age_binned") +
coord_fixed(sqrt(evals[2] / evals[1])) +
labs(col = "Binned Age")
out.dpcoa.log <- ordinate(pslog, method = "DPCoA")
evals <- out.dpcoa.log$eig
plot_ordination(pslog, out.dpcoa.log, color = "age_binned",
shape = "family_relationship") +
coord_fixed(sqrt(evals[2] / evals[1])) +
labs(col = "Binned Age", shape = "Litter")
plot_ordination(pslog, out.dpcoa.log, type = "species", color = "Phylum") +
coord_fixed(sqrt(evals[2] / evals[1]))
### pca
out.wuf.log <- ordinate(pslog, method = "PCoA", distance = "wunifrac")
abund <- otu_table(pslog)
abund_ranks <- t(apply(abund, 1, rank))
evals <- out.wuf.log$values$Eigenvalues
plot_ordination(pslog, out.wuf.log, color = "age_binned",
shape = "family_relationship") +
coord_fixed(sqrt(evals[2] / evals[1])) +
labs(col = "Binned Age", shape = "Litter")
plot_ordination(pslog, out.wuf.log, type = "species", color = "Phylum") +
coord_fixed(sqrt(evals[2] / evals[1]))
abund_df <- melt(abund, value.name = "abund") %>%
left_join(melt(abund_ranks, value.name = "rank"))
colnames(abund_df) <- c("sample", "seq", "abund", "rank")
abund_df <- melt(abund, value.name = "abund") %>%
left_join(melt(abund_ranks, value.name = "rank"))
colnames(abund_df) <- c("sample", "seq", "abund", "rank")
sample_ix <- sample(1:nrow(abund_df), 8)
ggplot(abund_df %>%
filter(sample %in% abund_df$sample[sample_ix])) +
geom_point(aes(x = abund, y = rank, col = sample),
position = position_jitter(width = 0.2), size = .7) +
labs(x = "Abundance", y = "Thresholded rank") +
scale_color_brewer(palette = "Set2")
abund_ranks <- abund_ranks - 329
abund_ranks[abund_ranks < 1] <- 1
ranks_pca <- dudi.pca(abund_ranks, scannf = F, nf = 3)
row_scores <- data.frame(li = ranks_pca$li,
SampleID = rownames(abund_ranks))
col_scores <- data.frame(co = ranks_pca$co,
seq = colnames(abund_ranks))
tax <- tax_table(ps)@.Data %>%
data.frame(stringsAsFactors = FALSE)
tax$seq <- rownames(tax)
main_orders <- c("Clostridiales", "Bacteroidales", "Lactobacillales",
"Coriobacteriales")
tax$order[!(tax$Order %in% main_orders)] <- "Other"
tax$order <- factor(tax$order, levels = c(main_orders, "Other"))
tax$otu_id <- seq_len(ncol(otu_table(ps)))
row_scores <- row_scores %>%
left_join(sample_data(pslog))
col_scores <- col_scores %>%
left_join(tax)
evals_prop <- 100 * (ranks_pca$eig / sum(ranks_pca$eig))
ggplot() +
geom_point(data = row_scores, aes(x = li.Axis1, y = li.Axis2), shape = 2) +
geom_point(data = col_scores, aes(x = 25 * co.Comp1, y = 25 * co.Comp2, col = order),
size = .3, alpha = 0.6) +
scale_color_brewer(palette = "Set2") +
facet_grid(~ age_binned) +
guides(col = guide_legend(override.aes = list(size = 3))) +
labs(x = sprintf("Axis1 [%s%% variance]", round(evals_prop[1], 2)),
y = sprintf("Axis2 [%s%% variance]", round(evals_prop[2], 2))) +
coord_fixed(sqrt(ranks_pca$eig[2] / ranks_pca$eig[1])) +
theme(panel.border = element_rect(color = "#787878", fill = alpha("white", 0)))
ps_ccpna <- ordinate(pslog, "CCA", formula = pslog ~ age_binned + family_relationship)
ps_scores <- vegan::scores(ps_ccpna)
sites <- data.frame(ps_scores$sites)
sites$SampleID <- rownames(sites)
sites <- sites %>%
left_join(sample_data(ps))
species <- data.frame(ps_scores$species)
dim(species)
species$otu_id <- seq_along(colnames(otu_table(ps))[-414])
species <- species %>%
left_join(tax)
library("ggrepel")
evals_prop <- 100 * ps_ccpna$CCA$eig[1:2] / sum(ps_ccpna$CA$eig)
ggplot() +
geom_point(data = sites, aes(x = CCA1, y = CCA2), shape = 2, alpha = 0.5) +
geom_point(data = species, aes(x = CCA1, y = CCA2, col = order), size = 0.5) +
geom_text_repel(data = species %>% filter(CCA2 < -2),
aes(x = CCA1, y = CCA2, label = otu_id),
size = 1.5, segment.size = 0.1) +
facet_grid(. ~ age_binned) +
guides(col = guide_legend(override.aes = list(size = 3))) +
labs(x = sprintf("Axis1 [%s%% variance]", round(evals_prop[1], 2)),
y = sprintf("Axis2 [%s%% variance]", round(evals_prop[2], 2))) +
scale_color_brewer(palette = "Set2") +
coord_fixed(sqrt(ps_ccpna$CCA$eig[2] / ps_ccpna$CCA$eig[1])) +
theme(panel.border = element_rect(color = "#787878", fill = alpha("white", 0)))
ggplot() +
geom_point(data = sites, aes(x = CCA1, y = CCA2), shape = 2, alpha = 0.5) +
geom_point(data = species, aes(x = CCA1, y = CCA2, col = order), size = 0.5) +
geom_text_repel(data = species %>% filter(CCA2 < -2),
aes(x = CCA1, y = CCA2, label = otu_id),
size = 1.5, segment.size = 0.1) +
facet_grid(. ~ family_relationship) +
guides(col = guide_legend(override.aes = list(size = 3))) +
labs(x = sprintf("Axis1 [%s%% variance]", round(evals_prop[1], 2)),
y = sprintf("Axis2 [%s%% variance]", round(evals_prop[2], 2))) +
scale_color_brewer(palette = "Set2") +
coord_fixed(sqrt(ps_ccpna$CCA$eig[2] / ps_ccpna$CCA$eig[1])) +
theme(panel.border = element_rect(color = "#787878", fill = alpha("white", 0)))
pls监督模型分类预测
library(caret)
setup_example(c("phyloseq", "ggplot2", "caret", "plyr", "dplyr"))
sample_data(pslog)$age2 <- cut(sample_data(ps0)$age, c(0, 100, 400))
dataMatrix <- data.frame(age = sample_data(pslog)$age2, otu_table(pslog))
trainingMice <- sample(unique(sample_data(pslog)$host_subject_id), size = 8)
inTrain <- which(sample_data(ps)$host_subject_id %in% trainingMice)
length(inTrain)
inTrain=inTrain[1:228]
training <- dataMatrix[inTrain,]
testing <- dataMatrix[-inTrain,]
plsFit <- train(age ~ ., data = training,
method = "pls", preProc = "center")
plsClasses <- predict(plsFit, newdata = testing)
table(plsClasses, testing$age)
# 换个算法继续尝试
library("randomForest")
rfFit <- train(age ~ ., data = training, method = "rf",
preProc = "center", proximity = TRUE)
rfClasses <- predict(rfFit, newdata = testing)
table(rfClasses, testing$age)
library("pls")
pls_biplot <- list("loadings" = loadings(plsFit$finalModel),
"scores" = scores(plsFit$finalModel))
class(pls_biplot$scores) <- "matrix"
pls_biplot$scores <- data.frame(sample_data(pslog)[inTrain, ],
pls_biplot$scores)
tax <- tax_table(ps)@.Data %>%
data.frame(stringsAsFactors = FALSE)
main_orders <- c("Clostridiales", "Bacteroidales", "Lactobacillales",
"Coriobacteriales")
tax$order[!(tax$order %in% main_orders)] <- "Other"
tax$order <- factor(tax$order, levels = c(main_orders, "Other"))
class(pls_biplot$loadings) <- "matrix"
pls_biplot$loadings <- data.frame(tax, pls_biplot$loadings)
ggplot() +
geom_point(data = pls_biplot$scores,
aes(x = Comp.1, y = Comp.2), shape = 2) +
geom_point(data = pls_biplot$loadings,
aes(x = 25 * Comp.1, y = 25 * Comp.2, col = order),
size = 0.3, alpha = 0.6) +
scale_color_brewer(palette = "Set2") +
labs(x = "Axis1", y = "Axis2", col = "Binned Age") +
guides(col = guide_legend(override.aes = list(size = 3))) +
facet_grid( ~ age2) +
theme(panel.border = element_rect(color = "#787878", fill = alpha("white", 0)))
rf_prox <- cmdscale(1 - rfFit$finalModel$proximity) %>%
data.frame(sample_data(pslog)[inTrain, ])
ggplot(rf_prox) +
geom_point(aes(x = X1, y = X2, col = age_binned),
size = .4, alpha = 0.6) +
scale_color_manual(values = c("#A66EB8", "#238DB5", "#748B4F")) +
guides(col = guide_legend(override.aes = list(size = 3))) +
labs(col = "Binned Age", x = "Axis1", y = "Axis2")
as.vector(tax_table(ps)[which.max(importance(rfFit$finalModel)), c("family", "genus")])
## [1] "Lachnospiraceae" "Roseburia"
impOtu <- as.vector(otu_table(pslog)[,which.max(importance(rfFit$finalModel))])
maxImpDF <- data.frame(sample_data(pslog), abund = impOtu)
ggplot(maxImpDF) +
geom_histogram(aes(x = abund)) +
facet_grid(age2 ~ .) +
labs(x = "Abundance of discriminative microbe", y = "Number of samples")
setup_example(c("igraph", "phyloseq", "phyloseqGraphTest", "ggnetwork", "intergraph"))
net <- make_network(ps, max.dist=0.35)
sampledata <- data.frame(sample_data(ps))
V(net)$id <- sampledata[names(V(net)), "host_subject_id"]
V(net)$litter <- sampledata[names(V(net)), "family_relationship"]
class(V(net)$aaa)
V(net)$litter = as.character(V(net)$aaa )
ggplot(net, aes(x = x, y = y, xend = xend, yend = yend), layout = "fruchtermanreingold") +
geom_edges(color = "darkgray") +
geom_nodes(aes(color = id, shape = litter)) +
theme(axis.text = element_blank(), axis.title = element_blank(),
legend.key.height = unit(0.5,"line")) +
guides(col = guide_legend(override.aes = list(size = .25)))
gt <- graph_perm_test(ps, "family_relationship",
grouping = "host_subject_id",
distance = "jaccard", type = "mst")
gt$pval
plot_test_network(gt) +
theme(legend.text = element_text(size = 8),
legend.title = element_text(size = 9))
plot_permutations(gt)
gt <- graph_perm_test(ps, "family_relationship",
grouping = "host_subject_id",
distance = "jaccard", type = "knn", knn = 1)
gt$pval
plot_test_network(gt) +
theme(legend.text = element_text(size = 8),
legend.title = element_text(size = 9))
plot_permutations(gt)
gt <- graph_perm_test(ps, "family_relationship",
grouping = "host_subject_id",
distance = "bray", type = "knn", knn = 2)
gt$pval
plot_test_network(gt) +
theme(legend.text = element_text(size = 8),
legend.title = element_text(size = 9))
plot_permutations(gt)
gt <- graph_perm_test(ps, "family_relationship", grouping = "host_subject_id",
distance = "bray", type = "threshold.nedges", nedges = 720,
keep.isolates = FALSE)
gt$pval
plot_test_network(gt) +
theme(legend.text = element_text(size = 8),
legend.title = element_text(size = 9))
plot_permutations(gt)
gt <- graph_perm_test(ps, "family_relationship", grouping = "host_subject_id",
distance = "bray", type = "threshold.nedges", nedges = 2000,
keep.isolates = FALSE)
gt$pval
plot_test_network(gt) +
theme(legend.text = element_text(size = 8),
legend.title = element_text(size = 9))
plot_permutations(gt)