[R/phyloseq] Neutral community model(NCM) 그리고 해석하기

작성: 2023-10-16-월

 

| Sloan's neutral community model(NCM) 이란?

NCM은  환경 생태쪽에서 발전한 개념으로 모든 종이 환경에서 중립적이라는 전제를 가지고 있다. 이 모델의 주요 특징은 아래와 같다. 

 

1. 중립성(Neutrality): 메타커뮤니티 내에서는 모든 종이 생태적으로 중립적이다. 이것은 모든 종이 기능적으로 동일한 경쟁능력과 적응성을 가진다는 것을 의미한다

 

2. 확률적 프로세스(Stochastic Processes): 모델은 확률적 또는 무작위 프로세스(종화, 멸종, 이주 및 이탈)를 포함한다.

 

3. 제로-합 동태학(Zero-Sum Dynamics): 커뮤니티의 개체 수는 상대적으로 일정하게 유지된다. 

 

즉 모델이 데이터와 일치할 수록 , OTUs들은 중립적이고 이주 가능성이 낮다는 것을 의미한다. 

 

 

1. 그림의 성분

(Burns AR et al., 2016)

 

✅ X축과 Y축

- 그림의 x축은 log 스케일로 상대적인 분포도를 나타내고, y축은 빈도를 나타낸다. 이 값에 따라 중립 모델의 95% 신뢰구간에 있는 모델은 대부분 빈도가 증가하면서 상대적 분포가 같이 늘어나는 중립적인 패턴을 보인다. 

 

✅ plot의 구성성분 

- 실선: 예측된 빈도 

- 점선 : 95% 신뢰구간

 

- 점 : 각 OTUs(ASVs)

> 위 점이 뜻하는 것은?

예측된 모델모다 빈도가 높은 OTUs는 초록색으로, 빈도가 낮은 OTUs는 주황색으로 표시되어 있다. 

- 초록색 점들은 예측보다 빈도가 높고, 이주 능력이 높으며, 다른 장소로 분산될 수 있다는 것을 나타낸다. 즉, host로 인해 유지되고 선택 되었음을 의미한다.

- 주황색은 예측보다 빈도가 낮으먀, host에 의해 분산이 제한되었다는 것을 나타낸다.

 

✅ 통계값

- R²  값: 중립 모델에 대한 적합도(모델이 데이터를 얼마나 잘 설명하는지 나타낸다.) 

- m 값: 추정 이동률,    m값이 높으면 이주가 빈번하게 일어난 다는 뜻으로, 여러 지역으로 이동하여 고르게 분포될 가능성이 높다. 

 

 

 

2. 그림 해석하기

NCM은 관찰된 OTU의 빈도와 메타 커뮤니티에서의 그들의 풍부도를 나타낸다. NCM은 데이터가 중립적인지, niche-based process를 따르는지 평가한다. 이때 upper line (fit 모델에서 95% 신뢰구간 이상)에 많이 존재하였을 때(R²), 이는 중립이 아닌 niche-based  process를 따른다는 가설을 뒷받침한다.NCM은 확률적인 프로세스를 강조하며, niche-based process는 종 간 경쟁과 상호작용을 강조한다.  

 

(Kim HJ et al., 2018)

 

위 데이터는 중국의 일부 도시에서 피부 마이크로바이옴의 NCM모델의 예측 결과를 나타낸다. 대도시인 Beijing은 R²값이  0.598로 비 대도시인 Hohhot 의 0.749보다 낮은 값을 가진다. 이는 Beijing이  niche-based 모델에 더 적합하다는 것을 말한다. 

 

즉 "대도시의 피부 미생물은 중립성이 약하며, 이는 신속하게 회복함으로써 교란으로부터 회복하기 어렵다는 것은 나타내며, 이는 미생물 군집의 취약성을 예측하는 요인일 수 있다."라고 위 논문에서는 전개하고 있다. 

 

 

| R에서 그리는 법

1. 예제 데이터 

- QIIME2 tutorial인 moving-picture의 데이터입니다. 사람의 양 손바닥, 혀, 장의 마이크로바이옴 데이터를 담고 있습니다. 현재 저장된 형식은 phyloseq object입니다. phylsoeq형식은 biom format을 차용하여 R에서 편리하게 다루도록 개발된 것으로,  metadata와 otu table, taxonomy 그리고 계통정보를 한 번에 저장할 수 있습니다. 

ps.rds

0.08MB

 

 

 

2. R코드

-  참고문헌 중에 "코드 참고  2"에 나오는 코드를 사용하겠습니다. 

(1의 코드가 더 짧지만, 2가 ggplot2를 사용하기 더 용이합니다)

 

데이터 불러오기

library(dplyr)
library(phyloseq)
library(devtools)
library(tibble)
library(ape)
library(vegan)
library(FSA)
library(eulerr)

# plot 
library(grid)
library(gridExtra)

# model
library(Hmisc)
library(minpack.lm)
library(stats4)

ps <- readRDS("./ps.rds")  # 770 taxa and 34 samples
ps.g <- subset_samples(ps, body.site %in% c("gut"))
ps.t <- subset_samples(ps, body.site %in% c("tongue"))

 

Gut microbiome 관찰

- 먼저 예제 샘플에서 우리가 추출한 데이터를 살펴봅시다. 

phyloseq = ps.g # 770 taxa and 8 samples

# 아래 ggplot2를 위해 metadata에 열 추가
meta <- sample_data(phyloseq)
meta$X.Sample <- rownames(meta)
sample_data(phyloseq) <- meta

## 각 샘플에 속한 read수 확인 
read_count = data.frame("count" = colSums(otu_table(phyloseq))) %>%
  rownames_to_column(var="X.Sample") %>%
  inner_join(data.frame(sample_data(phyloseq)), by="X.Sample") %>%
  arrange(-count) %>%
  mutate(X.Sample=factor(X.Sample, levels=X.Sample))

# 각 샘플에 대한 read수를 bar plot으로 그려보자. 
# 여기서는 2,000과10,000 read threshold를 나타낸다. 
ggplot(data=read_count, aes(x=X.Sample, y=log10(count), fill=body.site)) +
  geom_bar(stat="identity") +
  labs(x="Sample", y="Log10(Read count)") +
  geom_hline(yintercept=log10(10000)) +
  geom_hline(yintercept=log10(2000)) +
  theme(text = element_text(size=16),
        axis.text.x = element_blank())

 

분석에 들어가기 전에, NCM분석은 rarefy된 샘플을 사용합니다. 아래 코드에 제시되어있지만

1. 이를 위해선 각 샘플의 depth를 확인하고

2. 적절한 depth의 수를 기준으로 rarefy를 진행해야 합니다. 

 

위 그림을 보니, 모든 샘플이 맨 윗선(depth = 10,000)을 넘지 못하지만, 적어도 6000read정도 수준의 depth를 가집니다.  

우리는 9개 샘플 중 최소 depth인 6388 값을 기준으로 rarefy를 진행해야 합니다. 

(만약 최소 depth가 2000~3000이라면, 일정 depth(10,000)를 기준으로 샘플을 drop하여야 결과의 신뢰도를 높일 수 있습니다. 특정 depth는 딱 정해진 것은 없지만 적어도 만 리드가 넘어야 한다고 생각합니다.)

 

# Rarefy to an even depth
set.seed(42)  # setting seed for reproducibility
phyloseq.rare = rarefy_even_depth(phyloseq)

# Normalize read counts (this gives relative abundance)
phyloseq.norm = transform_sample_counts(phyloseq.rare, function(x) x/sum(x))

 

이제 모델을 계산하자. 

# phyloseq object에서 OTU table을 추출한다
OTU.table = t(otu_table(phyloseq.rare)) # communities/samples as rows and taxa as columns

# 각 OTU table에서 평균 abundance read수를 계산한다 
# rarefy이후이니 모두 6388 
N <- mean(apply(OTU.table, 1, sum)) # 6388

# 각 taxa에 대한 평균 abudance를 계산한다
p.m <- apply(OTU.table, 2, mean)
p.m <- p.m[p.m != 0]
p <- p.m/N
p.df = data.frame(p) %>%
  rownames_to_column(var="OTU")

# 각 taxa에 대한 occurrence frequency를 계산한다
OTU.table.bi <- 1*(OTU.table>0)
freq.table <- apply(OTU.table.bi, 2, mean)
freq.table <- freq.table[freq.table != 0]
freq.df = data.frame(OTU=names(freq.table), freq=freq.table)

# 위에서 계산한 데이터를 합친다
C <- inner_join(p.df,freq.df, by="OTU") %>%
  arrange(p)
# Remove rows with any zero (absent in either source pool or local communities). You already did this, but just to make sure we will do it again.
C.no0 <- C %>%
  filter(freq != 0, p != 0)

#Calculate the limit of detection
d <- 1/N

## 모델 생산
## Fit model parameter m (or Nm) using Non-linear least squares (NLS)
p.list <- C.no0$p
freq.list <- C.no0$freq
m.fit <- nlsLM(freq.list ~ pbeta(d, N*m*p.list, N*m*(1-p.list), lower.tail=FALSE), start=list(m=0.1))
m.ci <- confint(m.fit, 'm', level=0.95)
m.sum <- summary(m.fit)
m.coef = coef(m.fit)

freq.pred <- pbeta(d, N*coef(m.fit)*p.list, N*coef(m.fit)*(1-p.list), lower.tail=FALSE)
Rsqr <- 1 - (sum((freq.list - freq.pred)^2))/(sum((freq.list - mean(freq.list))^2))

# Get table of model fit stats
fitstats <- data.frame(m=m.coef, m.low.ci=m.ci[1], m.up.ci=m.ci[2], 
                       Rsqr=Rsqr, p.value=m.sum$parameters[4], N=N, 
                       Samples=nrow(OTU.table), Richness=length(p.list), 
                       Detect=d)

# Get confidence interval for predictions
freq.pred.ci <- binconf(freq.pred*nrow(OTU.table), nrow(OTU.table), alpha=0.05, method="wilson", return.df=TRUE)

# Get table of predictions
pred.df <- data.frame(metacomm_RA=p.list, frequency=freq.pred, 
                      frequency_lowerCI=freq.pred.ci[,2], 
                      frequency_upperCI=freq.pred.ci[,3]) %>%
  unique()

# Get table of observed occupancy and abundance
obs.df = C.no0 %>%
  rename(metacomm_RA = p, frequency=freq)

# 모델의 95% CI를 기준으로 어디에 속하는지 범주를 표시해 준다
obs.df <-  obs.df %>%
  left_join(pred.df, by = "metacomm_RA") %>%
  mutate(New_Column = case_when(
    frequency.x >= frequency_lowerCI & frequency.x <= frequency_upperCI ~ "fit",
    frequency.x > frequency_upperCI ~ "upper",
    frequency.x < frequency_lowerCI ~ "lower",
    TRUE ~ "NA"
  )) %>%
  select(-frequency.y, -frequency_lowerCI, -frequency_upperCI)

colnames(obs.df)[3] <- "frequency"

 

위 결과를 토대로 그림을 그려봅니다. 

 

fulldata.model.plot = 
  ggplot(data=obs.df) +
    geom_point(data=obs.df, aes(x=log10(metacomm_RA), y=frequency, color = New_Column), 
               alpha=.2, size=2) +
    geom_line(data=pred.df, aes(x=log10(metacomm_RA), y=frequency), color="black") + 
    geom_line(data=pred.df, aes(x=log10(metacomm_RA), y=frequency_lowerCI), linetype=2, color="black") + 
    geom_line(data=pred.df, aes(x=log10(metacomm_RA), y=frequency_upperCI), linetype=2, color="black") + 
    geom_text(data=fitstats, aes(label = paste("R^2 == ", round(Rsqr, 3))), 
              x=-3.8, y=0.95, size=4, parse=TRUE) +
    geom_text(data=fitstats, aes(label = paste("italic(m) ==", round(m, 3))), 
              x=-3.8, y=0.85, size=4, parse=TRUE) + 
    labs(x="Log10 abundance in\nmetacommunity", y="Frequency detected") +
    theme_bw() +
    theme(axis.line = element_line(color="black"),
          legend.position = "none",
          axis.title = element_text(size=14),
          axis.text = element_text(size=12))

fulldata.model.plot

 

 

Tongue microbiome 관찰

- 먼저 관찰하기 전에, 위에 적었던 과정을 모두 함수화 하자. input 파일은 phyloseq하나면 충분하다. 

########################
##    리드 수 확인    ##
########################

Check_read <- function(phyloseq){
  meta <- sample_data(phyloseq)
  meta$SampleID <- rownames(meta)
  sample_data(phyloseq) <- meta
  
  ## Check how many reads you have in each of the samples. This will tell you if you need to re-do anything
  # Get read counts and make a new dataframe with this data
  read_count = data.frame("count" = colSums(otu_table(phyloseq))) %>%
    rownames_to_column(var="SampleID") %>%
    inner_join(data.frame(sample_data(phyloseq)), by="SampleID") %>%
    arrange(-count) %>%
    mutate(SampleID=factor(SampleID, levels=SampleID))
  
  print(read_count)
  # Now plot read count for each sample. The horizontal line represents a 2000 read threshold
  ggplot(data=read_count, aes(x=SampleID, y=log10(count), fill=body.site)) +
    geom_bar(stat="identity") +
    labs(x="Sample", y="Log10(Read count)") +
    geom_hline(yintercept=log10(10000)) +
    geom_hline(yintercept=log10(2000)) +
    theme(text = element_text(size=16),
          axis.text.x = element_blank())
  

  
}
  
  
  
  
########################
##      모델 함수     ##
########################


NCM_fit <- function(phyloseq, rarefy = "min"){
  
  if (rarefy == "min") {
    # Rarefy to an even depth
    set.seed(42)  # setting seed for reproducibility
    phyloseq.rare = rarefy_even_depth(phyloseq)
    } else {
      phyloseq.rare = phyloseq
      }
  # Normalize read counts (this gives relative abundance)
  phyloseq.norm = transform_sample_counts(phyloseq.rare, function(x) x/sum(x))
  
  
  # Extract the OTU table from pthe phyloseq object
  OTU.table = t(otu_table(phyloseq.rare))
  
  # Calculate the number of individuals in the meta community (Average read depth)
  N <- mean(apply(OTU.table, 1, sum))
  
  # Calculate the average relative abundance of each taxa across communities
  p.m <- apply(OTU.table, 2, mean)
  p.m <- p.m[p.m != 0]
  p <- p.m/N
  p.df = data.frame(p) %>%
    rownames_to_column(var="OTU")
  
  # Calculate the occurrence frequency of each taxa
  OTU.table.bi <- 1*(OTU.table>0)
  freq.table <- apply(OTU.table.bi, 2, mean)
  freq.table <- freq.table[freq.table != 0]
  freq.df = data.frame(OTU=names(freq.table), freq=freq.table)
  
  #Combine
  C <- inner_join(p.df,freq.df, by="OTU") %>%
    arrange(p)
  # Remove rows with any zero (absent in either source pool or local communities). You already did this, but just to make sure we will do it again.
  C.no0 <- C %>%
    filter(freq != 0, p != 0)
  
  #Calculate the limit of detection
  d <- 1/N
  
  ##Fit model parameter m (or Nm) using Non-linear least squares (NLS)
  p.list <- C.no0$p
  freq.list <- C.no0$freq
  m.fit <- nlsLM(freq.list ~ pbeta(d, N*m*p.list, N*m*(1-p.list), lower.tail=FALSE), start=list(m=0.1))
  m.ci <- confint(m.fit, 'm', level=0.95)
  m.sum <- summary(m.fit)
  m.coef = coef(m.fit)
  
  freq.pred <- pbeta(d, N*coef(m.fit)*p.list, N*coef(m.fit)*(1-p.list), lower.tail=FALSE)
  Rsqr <- 1 - (sum((freq.list - freq.pred)^2))/(sum((freq.list - mean(freq.list))^2))
  
  # Get table of model fit stats
  fitstats <- data.frame(m=m.coef, m.low.ci=m.ci[1], m.up.ci=m.ci[2], 
                         Rsqr=Rsqr, p.value=m.sum$parameters[4], N=N, 
                         Samples=nrow(OTU.table), Richness=length(p.list), 
                         Detect=d)
  
  # Get confidence interval for predictions
  freq.pred.ci <- binconf(freq.pred*nrow(OTU.table), nrow(OTU.table), alpha=0.05, method="wilson", return.df=TRUE)
  
  # Get table of predictions
  pred.df <- data.frame(metacomm_RA=p.list, frequency=freq.pred, 
                        frequency_lowerCI=freq.pred.ci[,2], 
                        frequency_upperCI=freq.pred.ci[,3]) %>%
    unique()
  
  # Get table of observed occupancy and abundance
  obs.df = C.no0 %>%
    rename(metacomm_RA = p, frequency=freq)
  
  obs.df <-  obs.df %>%
    left_join(pred.df, by = "metacomm_RA") %>%
    mutate(New_Column = case_when(
      frequency.x >= frequency_lowerCI & frequency.x <= frequency_upperCI ~ "fit",
      frequency.x > frequency_upperCI ~ "upper",
      frequency.x < frequency_lowerCI ~ "lower",
      TRUE ~ "NA"
    )) %>%
    select(-frequency.y, -frequency_lowerCI, -frequency_upperCI)
  
  colnames(obs.df)[3] <- "frequency"
  
  
  
  fulldata.model.plot = 
    ggplot(data=obs.df) +
      geom_point(data=obs.df, aes(x=log10(metacomm_RA), y=frequency, color = New_Column), size=2, alpha = 0.8) +
      geom_line(data=pred.df, aes(x=log10(metacomm_RA), y=frequency), color="black") + 
      geom_line(data=pred.df, aes(x=log10(metacomm_RA), y=frequency_lowerCI), linetype=2, color="black") + 
      geom_line(data=pred.df, aes(x=log10(metacomm_RA), y=frequency_upperCI), linetype=2, color="black") + 
      geom_text(data=fitstats, aes(label = paste("R^2 == ", round(Rsqr, 3))), 
                x=-3.8, y=0.95, size=4, parse=TRUE) +
      geom_text(data=fitstats, aes(label = paste("italic(m) ==", round(m, 3))), 
                x=-3.8, y=0.85, size=4, parse=TRUE) + 
      labs(x="Log10 abundance in\nmetacommunity", y="Frequency detected") +
      theme_bw() +
      theme(axis.line = element_line(color="black"),
            legend.position = "none",
            axis.title = element_text(size=14),
            axis.text = element_text(size=12))
  
  return(fulldata.model.plot)

}

 

함수가 잘 만들어졌는지 확인해 보자!

Check_read(ps.g)
p.gut <- NCM_fit(ps.g)

 

 

 

이제 혀의 마이크로바이옴을 관찰해 보자

Check_read(ps.t)
p.tng <- NCM_fit(ps.t)

 

장보다는 혀의 마이크로바이옴 데이터가 모델의 적합성이 더 높은 것을 알 수 있습니다. 

 

 

 

 

 

 

 

참고

1.  Sloan's neutral model

   Sloan WT, Woodcock S, Lunn M, Head IM, Curtis TP. Modeling taxa-abundance distributions in microbial communities using environmental sequence data. Microb Ecol. 2007 Apr;53(3):443-55. doi: 10.1007/s00248-006-9141-x. Epub 2006 Dec 13. PMID: 17165121.

2. neutral community model을 gut microbiome에 적용한 논문(초기 코드 제공)

   Burns, A. R., Stephens, W. Z., Stagaman, K., Wong, S., Rawls, J. F., Guillemin, K., & Bohannan, B. J. (2016). Contribution of neutral processes to the assembly of gut microbial communities in the zebrafish over host development. The ISME journal, 10(3), 655–664. https://doi.org/10.1038/ismej.2015.142

3. 코드 참고 1

  Chen W, Ren K, Isabwe A, Chen H, Liu M, Yang J. Stochastic processes shape microeukaryotic community assembly in a subtropical river across wet and dry seasons. Microbiome. 2019 Oct 22;7(1):138. doi: 10.1186/s40168-019-0749-8. Erratum in: Microbiome. 2019 Nov 14;7(1):148. PMID: 31640783; PMCID: PMC6806580.

  https://github.com/Weidong-Chen-Microbial-Ecology/Stochastic-assembly-of-river-microeukaryotes

4. 코드 참고 2

   Barnett SE, Youngblut ND, Buckley DH. Soil characteristics and land-use drive bacterial community assembly patterns. FEMS Microbiol Ecol. 2020 Jan 1;96(1):fiz194. doi: 10.1093/femsec/fiz194. PMID: 31834372.

   https://github.com/seb369/landuse_comm_assembly

5. 응용연구 

  Kim HJ, Kim H, Kim JJ, Myeong NR, Kim T, Park T, Kim E, Choi JY, Lee J, An S, Sul WJ. Fragile skin microbiomes in megacities are assembled by a predominantly niche-based process. Sci Adv. 2018 Mar 7;4(3):e1701581. doi: 10.1126/sciadv.1701581. PMID: 29532031; PMCID: PMC5842045. (Kim HJ et al., 2018)

6. ChatGPT (what is Sloan neutral model?)