[R/SpiecEasi] OTUtable, phyloseq 데이터로 Network plot을 그려보자

[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) taxa bar 2) diversity plot_PCoA, bax-plot 등 3) network 분석 이 많이 보이는데 

네트워크 분석은 count table을 기반으로 각각 종의 상관관계를 선으로 연결하여 표시한 것으로 환경데이터에서 많이 사용된다

Network 분석을 위한 R package를 학습 해보자

 

 

• 예제 원문 
• 관련 논문 
• 예제 데이터 : amgut1.filt (American gut Progect Data from SpeicEasi package)
  • 289 sample(row) * 127 OTUs(column)

 

---

SpeicEasi Pipeline

즉 OTU table을 기반으로 여러가지를 한다

 

SpeicEasi 설치

library(devtools) # github에서 데이터 다운 받을 수 있게 하는 패키지
install_github("zdk123/SpiecEasi")
library(SpiecEasi)

 

기본 사용법

기본 단계는 아래와 같다

1. load the data and normalize counts to to common scale (min depth)
2. fit count margins to the model
3. generate a synthetic network
4. generate some synthetic data
5. clr transformation
6. inverse covariance estimation along a lambda (sparsity) path
7. stability selection using the StARS criterion
8. evaluate performance

실제 데이터는 1-4단계를 건너뛰어라

 

그럼 일단 예제 데이터에서 1-4단계를 맛보자 

분석 데이터 살펴보기 

# 데이터 불러오기 
data(amgut1.filt)

- SpiedEasi분석의 input 파일은 이와 같은 feature table(count table) 이다

- qiime에서 Dada2나 Deblur로 만든 feature-table.biom을 txt로 변환후 행렬 변환만 하면 SpeicEasi의 input파일을 만들 수 있다

amgut1.filt.tsv

0.10MB

 

 

 

 

 

Analysis of American Gut data

자 1-4단계 없이 바로 데이터에 적용시켜 보자 normalization이 spiec.easi 함수 내부에서 시행된다는걸 잊지 말자.

또한 spiec.easi 에서 rep.num(stars repetitions)옵션에서 큰 수를 사용해야 한다

 

 네트워크 분석에 사용되는 총 3가지 방법(md/ glasso/ sparcc)을 비교해보자 (15분 걸림)

se.mb.amgut <- spiec.easi(amgut1.filt, method='mb', lambda.min.ratio=1e-2,
                          nlambda=20, pulsar.params=list(rep.num=50))

se.gl.amgut <- spiec.easi(amgut1.filt, method='glasso', lambda.min.ratio=1e-2,
                          nlambda=20, pulsar.params=list(rep.num=50))                          

sparcc.amgut <- sparcc(amgut1.filt)
sparcc.graph <- abs(sparcc.amgut$Cor) >= 0.3 # abs : 절댓값 씌워줌 
diag(sparcc.graph) <- 0
library(Matrix)
sparcc.graph <- Matrix(sparcc.graph, sparse=TRUE) ## Create igraph objects


se.mb.amgut
# Pulsar-selected refit of sparseiCov 
# Path length: 20 
# Graph dim:   127 
# Criterion:
#   stars... sparsity 0.0179

se.gl.amgut
# Pulsar-selected refit of sparseiCov 
# Path length: 20 
# Graph dim:   127 
# Criterion:
#   stars... sparsity 0.0221

str(sparcc.amgut)
# List of 2
#  $ Cov: num [1:127, 1:127] 5.1095 -0.0176 -0.0268 0.2589 -1.2604 ...
#  $ Cor: num [1:127, 1:127] 1 -0.00469 -0.00476 0.06594 -0.21757 ...

 

바로 plot으로 나타낼 수는 없음으로 변환단계를 거쳐주자

#시각화를 위한 객체 생성 완료 
ig.mb     <- adj2igraph(getRefit(se.mb.amgut))

ig.gl     <- adj2igraph(getRefit(se.gl.amgut))

ig.sparcc <- adj2igraph(sparcc.graph)

ig.mb 
# IGRAPH 42238df UNW- 127 144 -- 
# + attr: name (v/n), weight (e/n)
# + edges from 42238df (vertex names):
#   [1]  1-- 58  2-- 15  2-- 68  2--105  4-- 14  4-- 22  4-- 39  4-- 44  4-- 56  4-- 93  4--121  5-- 14  5-- 44  5-- 90  5--103
#  [16]  5--104  5--116  5--121  6-- 89  8-- 82 10-- 20 11-- 87 11--110 12-- 19 12-- 60 12-- 65 12-- 82 12--119 13-- 22 13-- 41
#  [31] 13-- 71 13-- 89 15-- 68 16-- 46 16-- 60 16-- 82 16-- 84 16--104 16--119 19-- 50 22-- 50 23-- 89 23--123 23--126 24-- 47
#  [46] 24-- 86 25-- 35 25-- 59 26-- 61 26-- 72 27-- 32 27--120 28-- 33 28-- 43 29-- 38 29-- 48 29-- 66 29--114 29--125 30-- 77
#  [61] 31-- 77 32-- 75 32-- 95 33-- 43 34-- 50 34-- 60 34-- 65 35-- 58 35-- 86 35--125 38-- 43 38-- 48 38-- 74 38--101 38--114
#  [76] 39-- 90 39--108 41-- 66 41-- 71 42-- 57 42--122 43-- 85 43--109 43--125 44-- 46 44-- 76 44--103 44--104 44--116 44--119
#  [91] 45-- 71 46-- 76 46--116 47-- 86 48--127 50-- 82 50--104 51-- 54 52-- 67 53-- 80 53--103 54-- 62 56-- 92 57--122 60-- 65
# [106] 61-- 72 62--105 62--112 62--120 63--107 64--122 66-- 77 66--125 67-- 75 68-- 75 68--105 69-- 78 69-- 94 71-- 89 74--114
# + ... omitted several edges

바로 plot으로 나타낼 수는 없음으로 변환단계를 거쳐주자

 

igraph package는 위의 1--3, 3--4 같은 형태의 포멧을 각자 연결해주어 plot을 그려주는 시각화 도구이다

library(igraph) # 네트워크를 시각화 하는 패키지
## set size of vertex proportional to clr-mean
vsize    <- rowMeans(clr(amgut1.filt, 1))+6 
am.coord <- layout.fruchterman.reingold(ig.mb)

par(mfrow=c(1,3)) # 3 화면으로 나누어서 보여주겠다
plot(ig.mb, layout=am.coord, vertex.size=vsize, vertex.label=NA, main="MB")
plot(ig.gl, layout=am.coord, vertex.size=vsize, vertex.label=NA, main="glasso")
plot(ig.sparcc, layout=am.coord, vertex.size=vsize, vertex.label=NA, main="sparcc")

잘은 몰라도 일단 선이 많이 생길 수록 모든 생물이 연결되어있다는 뜻이니 결과가 잘 나온것으로 볼 수 있음

 

 

Working with phyloseq

이미 Phyloseq format으로 만들어진 파일이 phyloseq package안에 존재 한다

이 기본 예제 데이터를 사용해보자

library(phyloseq)

## Load round 2 of American gut project
data('amgut2.filt.phy')

amgut2.filt.phy
saveRDS(se.mb.amgut2, "~/se.mb.amgut2.rds")
# R에서 읽을땐 다시  a <- readRDS("~/se.mb.amgut2.rds") 하면 된다

 

 

phyloseq 데이터를 바로 넣을 수 있다.

장점이 있는데 더 많은 옵션을 부과해서 시각화 할 수 있지만,

그만큼 네트워크를 구축하는데 시간이 많이 소요된다

se.mb.amgut2 <- spiec.easi(amgut2.filt.phy, method='mb', lambda.min.ratio=1e-2,
                           nlambda=20, pulsar.params=list(rep.num=50))
                                                 
ig2.mb <- adj2igraph(getRefit(se.mb.amgut2),  vertex.attr=list(name=taxa_names(amgut2.filt.phy)))

se.mb.amgut2.rds

5.81MB

 

plot_network() 에서 taxa를 기준으로 color는 Rank3(Class level)에 따라 관찰해보자

plot_network(ig2.mb, amgut2.filt.phy, type='taxa', color="Rank3")