bioinformatics

OTUs三元图 + 富集分区着色

2020-05-19  本文已影响0人  吴十三和小可爱的札记

1. 简介

微生物多样性分析中,普通三元图不同的点代表不同的分otus(或分类水平),点的大小代表平均丰度。不仅如此,还可以根据各OTU在各根系微环境中的丰度数据进行统计检验后,得出各OTU分别在哪种微环境中显著富集,并据此在图中以不同颜色的点表示,此时三元图不仅表现出OTU或者物种的分布,还包含统计显著结果。

数据样式:常见OTU丰度表,列为样本,行为OTU,交叉区域为每种OTU在各样本中的丰度。

分组信息:

2.显著性检验

结合edgeR{}统计组间差异OTU,获得不同分组的富集信息。

rm(list = ls())

file <-  "otu.file"
otu <- read.delim(file = paste0(file, "otu_table.txt"),
                              row.names = 1, sep = '\t', 
                              stringsAsFactors = FALSE, 
                              check.names = FALSE)

group_data = read.delim(file = group.txt,
                              row.names = 1, sep = '\t', 
                              stringsAsFactors = FALSE, 
                              check.names = FALSE)

library(edgeR)

# Construct DGEList
dge_list <- DGEList(counts = otu, group = group_data$groups)


# Remove the lower abundance/(cpm, rpkm)
keep <- rowSums(dge_list$counts) >= 0
dge_keep <- dge_list[keep, ,keep.lib.sizes = FALSE]

# scale the raw library sizes dgelist
dge <- calcNormFactors(dge_keep)

# Construct Design Matrices
design.mat <- model.matrix(~ 0 + dge$samples$group)

# fit the GLM
Common_GLM <- estimateGLMCommonDisp(dge, design.mat)
Tagwise_GLM <- estimateGLMTagwiseDisp(Common_GLM, design.mat)

fit <- glmFit(Tagwise_GLM, design.mat)

# conducts likelihood ratio tests
lrt_A_B <- glmLRT(fit, contrast = c(1, -1, 0))
lrt_A_C <- glmLRT(fit, contrast = c(1, 0, -1))

# extant the main information: logFC, PValue, FDR
# data <- topTags(lrt, n = nrow(dgelist$counts))
# data <- as.data.frame(data)

# Identify the differentially OTUs though defalut parms
de_A_B <- decideTestsDGE(lrt_A_B, adjust.method = "fdr", 
                                p.value = 0.5)
de_A_C <- decideTestsDGE(lrt_A_C, adjust.method = "fdr", 
                                p.value = 0.5)

#  1 indicating each otu is classified as significantly significant up-regulated 
rich_A <- rownames(otu)[de_A_B == 1 & de_A_C == 1]

# construct enrich_A data
enrich_A <- data.frame(rownames = rich_A, 
                    enrich = rep("A", length(rich_A)))

# 重复上述步骤获得enrich_C, enrich_D
enrich_index <- rbind(enrich_A, enrich_C, enrich_D)

可视化

三元图 - ggtern

上一篇 下一篇

猜你喜欢

热点阅读