작성 시작 : 2023-05-16 ~ 2023-05-17
Regression 은 주로 연속형 데이터의 예측에 사용된다. 아래는 Microbiome 분야에서 사용되는 regression plot의 종류를 알아보자
사실 원하는 두 데이터만 있으면 linear regression은 뭐든 그릴 수 있다.
| 목표
마이크로바이옴 데이터로 여러 종류의 regression plot을 그리고 이를 함수를 만들어보자.
| 예제 데이터
- qiime2 moving pictures Tutorial에 나오는 데이터로, 사람의 4 부위에 해당하는 마이크로바이옴 데이터를 담고 있다.
- 이 데이터는 phyloseq 데이터로, biom형식으로 구성되어 있다. [참고]
library(phyloseq) # 사용할 데이터를 다루는 패키지
library(ggplot2) # 시각화 패키지
library(dplyr) # 데이터를 다룰 때 사용
library(ggpmisc) # 통계 결과값을 나타낼 때 사용
# import data 파일 읽어오기
ps <- readRDS("./ps.rds")
ps
# phyloseq-class experiment-level object
# otu_table() OTU Table: [ 770 taxa and 34 samples ]
# sample_data() Sample Data: [ 34 samples by 8 sample variables ]
# tax_table() Taxonomy Table: [ 770 taxa by 7 taxonomic ranks ]
# phy_tree() Phylogenetic Tree: [ 770 tips and 768 internal nodes ]
# Show metadata 메타데이터 살펴보기
sample_data(ps)
# read count to relative abundance 상대적 분포로 바꾸기 (각 read수를 샘플의 총 read 수로 나누어 줌)
ps.rel <- transform_sample_counts(ps, function(x) x/sum(x) )
| 시각화
1. Read 수에 따른 Alpha diversity 변화
우리는 phyloseq을 사용하니, 이를 사용해서 alpha diversity(shannon index)를 계산해 보자.
alpha_div <- estimate_richness(ps.rel, measures = "Shannon")그다음 phyloseq에서 otu table 가져오고, 각 샘플의 read 수의 총합을 담고 있는 column을 추가(dplyr::mutate())해보자.
otu <- otu_table(ps) %>% t() %>% as.data.frame() %>%
mutate(Sample_read = rowSums(.))위 함수에서 t()를 해주는 이유는 우리가 사용하는 데이터의 otu_table은 샘플의 이름이 열의 이름이며, OTU 가 각 행에 위치하기 때문이다. 만약 본인의 데이터가 그 반대의 경우라면 t()를 사용하지 않아도 된다.
alpha diversity와 otu테이블의 데이터를 샘플이름(rownames)을 기준으로 합쳐주자.
alpha_read <- merge(otu, alpha_div, by = "row.names")
위 데이터에서 regression model을 만들어 보자.
model <- lm(Sample_read~Shannon, alpha_read)
model
# Call:
# lm(formula = Sample_read ~ Shannon, data = alpha_read)
#
# Coefficients:
# (Intercept) Shannon
# 1459 1035
이 데이터를 사용해서 그림을 그리고, regression model을 표시해 보자. 기본 값으로 stat_smooth를 추가했을 때, confidence interval(회색영역)은 95%를 기준으로 그려진다.
ggplot(alpha_read, aes(Sample_read, Shannon)) +
geom_point() +
stat_smooth(color = "blue", method="lm")
여기서 우리는 그림 안에 R^2과 p-value값을 추가하고 싶다. 이럴 때는 수동으로 추가할 수 도 있지만, 출처[3]에 나오는 함수를 사용해 보자.
ggplot(alpha_read, aes(Sample_read, Shannon)) +
geom_point() +
stat_smooth(color = "blue", method="lm") +
stat_poly_eq(formula = as.formula(y ~ x),
label.x = "centre",
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., sep = "~~~")),
color= "red",
parse = TRUE)+
stat_fit_glance(method = 'lm',
method.args = list(formula = y ~ x),
label.x = label_pos,
aes(label = paste("~italic(p) ==", round(..p.value.., digits = 3),
"~italic(R)^2 ==", round(..r.squared.., digits = 2),
sep = "~")),
color= "red",
parse=TRUE)+
theme_minimal() +
labs(title = "Regression plot ", y = "shannon", x = "Read Number")
위 그래프에 나와있는 p-value과 R^2은 아래 regression model에서 구한 값과 같다.
summary(model)$coefficients[2,4] # p-value
summary(model)$r.squared # R squared
자 그렇다면 저 복잡한 시각화코드를 함수로 만들어 보자. 함수로 만드는 이유는 간단한 함수로 긴 코드를 구현 가능하게 해 주기 때문이다. ggplot2를 사용해서 함수화 할 때는 aes() 부분을 약간 수정할 필요가 있다.
reg_plot <- function(df, x = x, y = y, Ylab = NULL, Xlab = NULL, Title = NULL, label_pos = "right") {
ggplot(df, aes_string(x, y)) + # aes -> aes_string
geom_point() +
stat_smooth(color = "blue", method="lm") +
stat_poly_eq(formula = as.formula(y ~ x),
label.x = "centre",
eq.with.lhs = "italic(hat(y))~`=`~",
aes(label = paste(..eq.label.., sep = "~~~")),
color= "red",
parse = TRUE)+
stat_fit_glance(method = 'lm',
method.args = list(formula = y ~ x),
label.x = label_pos,
aes(label = paste("~italic(p) ==", round(..p.value.., digits = 3),
"~italic(R)^2 ==", round(..r.squared.., digits = 2),
sep = "~")),
color= "red",
parse=TRUE)+
theme_minimal() +
labs(title = Title, y = Ylab, x = Xlab)
}
이제 만든 함수를 사용해서 데이터를 시각화해 보자.
reg_plot(alpha_read, "Sample_read", "Shannon")동일한 결과가 출력된다. 이후의 그림도 이 reg_plot을 사용해 보자.
2. 메타 데이터에 대한 Alpha diversity 변화
alpha_div <- estimate_richness(ps.rel, measures = "Shannon")
alpha_meta <- merge(sample_data(ps.rel), alpha_div, by = "row.names")
reg_plot(alpha_meta, "days.since.experiment.start", "Shannon",
Ylab = "Alpha diversity", Xlab = "Time", Title = "Regression plot 02")
그림만 보면 metadata의 days.since.experiment.start 와 alpha diversity가 연관이 있어 보인다. 그러나 p-value값은 유의하지만 R^2 값은 낮다. 이 데이터는 사람의 4개의 부위에서 채취된 데이터로 각 데이터의 diversity가 다르다. 그러므로 여기서 유의미한 결과를 해석하기는 어렵다.
만약 항생제를 투여하고 이를 일정 기간 동안 주기적으로 관찰한 연구의 경우, regression plot을 통해 diversity의 결과를 유의미하게 볼 수 있을 것이다.
3. 특정 taxa와 메타데이터와 비교
otu <- as.data.frame(otu_table(ps.rel))
otu <- otu %>% mutate(Sum = rowSums(.))
otu %>% arrange(-Sum)여기서 풍부도가 높은 d29fe3c70564fc0f69f2c03e0d1e5561를 사용해서 비교해 보자.
otu_table <- as.data.frame(t(otu_table(ps.rel)))
metadata <- as.data.frame(sample_data(ps.rel))
# tax <- as.data.frame(tax_table(ps.rel))
# tax["d29fe3c70564fc0f69f2c03e0d1e5561", "Species"]
otu_meta <- merge(otu_table, metadata, by = "row.names")
reg_plot(otu_meta, "days.since.experiment.start", "d29fe3c70564fc0f69f2c03e0d1e5561",
Ylab = "OTU", Xlab = "Time", Title = "Regression plot 03")
4. 메타데이터 간 회귀분석
meta2 <- sample_data(ps.rel) %>% as.data.frame()
meta2$Random <- sample(1:300, 34)
reg_plot(meta2, "days.since.experiment.start", "Random",
Ylab = "Random data", Xlab = "Time", Title = "Regression plot 04")마땅히 비교할 데이터가 없으니 sample()함수를 통해 랜덤한 데이터를 만들어 준다.
결과적으로 정말 아무 의미 없어 보이는 plot이 완성됐다. BMI와 혈당지수 등 의학적으로 연관이 큰 데이터를 이용하는 것이 의미 있겠다.
5. 각 taxa 간 비교
otu <- as.data.frame(t(otu_table(ps.rel)))
x <- "3c9c437f27aca05f8db167cd080ff1ec" # "Prevotella melaninogenica"
y <- "1d2e5f3444ca750c85302ceee2473331" # Haemophilus parainfluenzae
reg_plot(otu,
x = "`3c9c437f27aca05f8db167cd080ff1ec`", # "Prevotella melaninogenica
y = "`1d2e5f3444ca750c85302ceee2473331`", # Haemophilus parainfluenzae
Ylab = "OTU1", Xlab = "OTU2", Title = "Regression plot 05")변수의 이름이 숫자로 시작할 때 종종 에러가 나기 때문에 ` `(~키)을 사용해서 변수를 감싸주자.
생각보다 그럴듯한 plot이 완성되었다!
아마도 저 두 개의 taxa가 풍부도가 높은 친구들이라 p-value값에서 유의성이 나타난 것 같다.
| 참고
1. chatGTP
