[Microbiome분석 단계]
1) Preprocessing
2) OTU clustering
3) Taxonomy classification
4) Diversity
- Alpha diversity
- Beta diversity
6) Differential abundance analysis
7) Functional analysis
8) Network analysis(Correlation)
+) Machine learning
🟦 목차
1. phyloseq basic function
2. Abundance 1% < 은 Other로 치환한 bar plot
3. Top10과 나머지는 Other로 처리한 bar plot
4. 다 귀찮고 그냥 Top n골라서 보고 싶다면?
🟦 Phyloseq 이란?
R 패티지 의 하나로 micobiome데이터를 분석 및 시각화하는데에 초점이 맞춰져 있다.
기본 구조는 Taxonomy파일과 OTU파일이며, 추가적으로 sample metadata와 phylogeny tree를 추가로 합쳐서 구성할 수도 있다.
일단 GlobalPatterns데이터를 활용해보자
library(phyloseq)
library(ggplot2)
library(dplyr)
data("GlobalPatterns")
GlobalPatterns
# phyloseq-class experiment-level object
# otu_table() OTU Table: [ 19216 taxa and 26 samples ]
# sample_data() Sample Data: [ 26 samples by 7 sample variables ]
# tax_table() Taxonomy Table: [ 19216 taxa by 7 taxonomic ranks ]
# phy_tree() Phylogenetic Tree: [ 19216 tips and 19215 internal nodes ]
🟦 Taxonomy bar plot 이란?
각 sample들에서 어떠한 미생물 종이 발견되었는지를 한눈에 볼 수 있게 그린 bar plot을 말한다.
각 샘플마다 시퀀싱된 read수가 모두 다르기 때문에, 각 샘플의 read생성수의 총 양으로 나눈 상대적 분포로 나타내게 된다.
🟦 Taxonomy bat plot을 그리는 방법
plot을 그리기 전에 각 샘플의 총 read수의 합으로 각 otu를 나누어 주도록 한다.
GlobalPatterns.rAB = transform_sample_counts(GlobalPatterns, function(x) x / sum(x) )
otu_table(GlobalPatterns.rAB)
#Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq' Also defined by ‘tidytree’ 는 무시해도 됨
이제 taxonomy bar plot을 그려볼 건데, 계산 속도를 위해 우리는 Phylum level로만 plot을 그려볼 것이다.
1. phyloseq basic function
- phyloseq의 기본 함수를 사용한 bar plot
phyloseq::plot_bar(physeq.rAB, fill = "Phylum") #Phylum taxa를 기준으로 보겠다
sample_data(GlobalPatterns.rAB)로 본 결과 SampleType열을 보니 각 샘플의 샘플링 지역을 볼 수 있다.
각 지역별로 나누어서 보고 싶다면? facet_wrap(~SampleType , scales = "free")를 추가하면 된다
phyloseq::plot_bar(GlobalPatterns.rAB, fill = "Phylum")+
geom_bar(aes(color = Phylum), stat = "identity", position = "stack") + # fill = Phylum
labs(x = "", y = "Relative Abundance\n") +
facet_wrap(~SampleType , scales = "free") +
theme(panel.background = element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())
2. Abundance 1% < 은 Other로 치환한 bar plot
- Major한 taxa를 집중 적으로 보고 싶을 때 사용
- 1%이하의 abundance를 가진 taxa를 Minor 하다고 보았을 때, 1% 미만을 Others로 처리한다
# 1) ggplot 을 사용하기 위해 데이터 해체하기
GlobalPatterns.rAB.melt <- tax_glom( GlobalPatterns.rAB, taxrank = 'Phylum')
# Error 메세지 무시해도 됨
GlobalPatterns.rAB.melt <- psmelt(GlobalPatterns.rAB.melt)
# 2) < 1% abundance 인 phylum 의 이름을 "< 1% Abundance"로 모두 바꾸기
GlobalPatterns.rAB.melt$Phylum[GlobalPatterns.rAB.melt$Abundance < 0.01 ] <- "< 1% Abundance"
# 3) 파래트를 지정하기 위해 남은 phylum의 종류 세기
Count = length(unique(GlobalPatterns.rAB.melt$Phylum))
Count # 20
# 4) 그냥 그리면 알파벳 순으로 그려짐 -> 빈도 순으로 지정하고 싶다면 아래 단계 적용
Order = tapply(GlobalPatterns.rAB.melt$Abundance, GlobalPatterns.rAB.melt$Phylum, sum) %>% sort(decreasing = T)%>% rownames() # 빈도 순으로 정렬
GlobalPatterns.rAB.melt$Phylum <- factor(GlobalPatterns.rAB.melt$Phylum, levels = c(Order[Order!="< 1% Abundance" ], "< 1% Abundance" ))
# 5) 색 지정하기
library(RColorBrewer)
qual_col_pals = brewer.pal.info[brewer.pal.info$category == 'qual',]
col_vector = unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))
col=sample(col_vector, Count)
# 6) stacked bar plot 그리기
p <- ggplot(data=GlobalPatterns.rAB.melt, aes(x=Sample, y=Abundance, fill=Phylum)) +
geom_bar(aes(fill = Phylum), stat="identity", position="fill") +
scale_fill_manual(values = col) +
theme(legend.position="bottom") +
theme(axis.text.x = element_text(angle = 90))+
theme(panel.background = element_blank())
p
3. Top10과 나머지는 Other로 처리한 bar plot
- Major한 Top 10개의 taxa만 보고 싶을 때 사용
- 보통 소수의 taxa가 대부분의 taxa를 차지할 때 사용
- Top 10을 빼고 모두 Others로 처리한다
GlobalPatterns.rAB.melt <- GlobalPatterns.rAB %>% tax_glom(taxrank = 'Phylum') %>%
psmelt()
Order = tapply(GlobalPatterns.rAB.melt$Abundance, GlobalPatterns.rAB.melt$Phylum, sum) %>%
sort(decreasing = T) # 빈도 순으로 정렬
Top10 <- names(Order[c(1:10)])
# NA 값 제거
table <-tax_table(GlobalPatterns)
table[is.na(table)] <- "_"
tax_table(GlobalPatterns) <- table
# melt
GlobalPatterns.rAB = transform_sample_counts(GlobalPatterns, function(x) x / sum(x) )
GlobalPatterns.rAB.melt <- GlobalPatterns.rAB %>% tax_glom(taxrank = 'Phylum') %>% psmelt()
# Top10만 고르기
Order = tapply(GlobalPatterns.rAB.melt$Abundance, GlobalPatterns.rAB.melt$Phylum, sum) %>% sort(decreasing = T) # 빈도 순으로 정렬
Top10 <- names(Order[c(1:10)])
# Top10을 제외하곤 Others로 바꾸기
GlobalPatterns.rAB.melt$Phylum[!(GlobalPatterns.rAB.melt$Phylum %in% Top10) ] <- "Others"
# Color 11개 색 지정하기
Count <- 11
qual_col_pals = brewer.pal.info[brewer.pal.info$category == 'qual',]
col_vector = unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))
col=sample(col_vector, Count)
# taxonomy bar plot 그리기
p <- ggplot(data=GlobalPatterns.rAB.melt, aes(x=Sample, y=Abundance, fill=Phylum)) +
geom_bar(aes(), position="fill", stat="identity") +
scale_fill_manual(values = col) +
theme(strip.background = element_blank(),
legend.position="bottom" ,
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank()) + # x axis rotation
theme_classic() +
theme(axis.text.x.bottom = element_text(angle = 45, vjust = 1, hjust=1),
axis.text = element_text(size = 7))+
theme(legend.position="right")
p
4. 다 귀찮고 그냥 Top n골라서 보고 싶다!
- 이때 "MicrobiotaProcess"라는 아주 편리한 package가 있으니 사용하면 된다.
- 기본 default값으로 돌리면 top30인 barplot을 그려준다.
- 지금은 topn=10으로 설정하여 top10과 나머지를 others로 처리한 plot을 그려보자
BiocManager::install("MicrobiotaProcess") # https://yulab-smu.top/MicrobiotaProcessWorkshop/articles/MicrobiotaProcessWorkshop.html
library(MicrobiotaProcess)
phylumtaxa <- get_taxadf(obj=GlobalPatterns, taxlevel=2) # taxlevel=2 <- phylum
pclas <- ggbartax(obj=phylumtaxa, topn = 10) + # Top 10만 보겠다
xlab(NULL) +
ylab("relative abundance (%)")
pclas
이렇게 총 4가지의 taxonomy barplot을 그리는 법을 알아보았다.
taxonomy의 분포를 보는건 bar plot뿐만이 아니라 pie, heatmap 등 여러 가지 시각화 자료가 존재한다.
나머지는 다른글에서 적어보겠다.
