Network analyses using R 

組成データ/Compositional data

定数和制約 :

正の数かつ合計が定数のデータ(Aitchison 1986)

① 対数比解析 (logratio analysis) ... 変数変換し、組成データを実数へ変換

② 単体解析 (simplicial analysis) ... 組成データを距離空間として再定義する

(太田&新井 2006）

分析方法

：Sparse Cooccurrence Network Investigation for Compositional Data [相関(correlation/Cooccurrence)に基づく]）(Shaffer et al. 2020) 

*Single-Cell rEgulatory Network Inference and Clusteringとは別物

仮定条件

ref. : SCNIC解説（日本語）

Compositional dataの解説（日本語）: (Ohta and Arai 2006)


：Sparese InverseE Covariance estimation for Ecological Association and Statistical Inference [共分散(corvariance)に基づく] *SparCCを含む

以下を含むpileline
  • data transformationを含む（input=normalizeしない）
  • sparse inverse covariance estimation
  • model selection

λ決定の検討が必要

  • SpiecEasiでは、①最適λを決定し、②そのモデルをサンプルにfittingする（getRefit()）
  • parameters: lambda.min.ratio, nlambda ...R package huge()のヘルプに記述がある（GitHub-#57）
  • Optimum lambda①: （GitHub-#190）
    • lambda.min.ratioが"too low"の場合、全ネットワークは連結している
    • [res.spieceasi]$select$stars$summaryがthresh値（閾値=defaultで0.05）を境界とする値を探す。直線的に増加する
spiec.easiのパラメータ例）
      • lambda.min.ratio=0.01
      • nlambda=100
      • pulsar.params=list(thresh=0.05, rep.num=20,seed=888,ncores=9)
  • Optimum lambda②: （GitHub-#85）
    • 高次元（多OTUs）& 少数サンプルはSpiecEasiには向かないので、OTU数を減らしてmodelのパラメータを減らす必要がある。
    • rep.run...少なすぎると正確なパラメータ決定ができない
  • λ：小=複雑 ↔︎ 大=単純

method: 

pulsar package@R ref. Kurtz's web

Terms: 
  • 隣接行列：有限グラフを表す正方行列（0以外の要素で表す）
  • precision matrix=inverse covariance matrix：精度行列
  • neighborhood selection (Meinshausen and Bühlmann 2010)
  • sparse inverse covariance estimation (Banerjee et al 2008; Yuan and Lin 2007)
  • Graphical Lasso（GLASSO）(Friedman et al 2008) lasso=輪投げ
  • Quadratic approximation for sparse Inverse estimation (QUIC: Hsieh et al 2014)
  • glasso:
  • mb (meinshausen-buhlmann's neighborhoood selection)

グラフのモデル選択方法：regularization parameter="ℷ (lambda)": model complexityをcalibrationするので重要
  1. Information criteria (IC) ... AIC, Bayesian IC (BIC)（Yuan and Lin 2007; Foygel and Drton 2010）
☞ 低次元（少数OTUs）に向くが、高次元には向かない(Liu et al. 2010)
  1. Stability Approach to Regularization Selection (StARS) (Liu et al. 2010: 高次元で良い再現性)
☞ Networkの希薄化. 実証的. random subsamplingによる. stabilityに関するβを固定するので、正規化よりはstabilityに依存する調整となる☞ StARSの短所: ①計算コスト： N個のsubsample setでneighborhooodやinverse covariance selectionを実行（Nは指定可能）②βの最適化：edge stabilityは未知のgraph形態の影響を強く受けるので、不変的なβの値はない⬅︎ 改善：Müller et al 2016① N=2（Bounded StARS: 正規化の上限と下限の2点）に指定（一定以上に複雑なグラフの計算を省く）② "edge stability"をsubgraph (graphlet) stabilityへ一般化し、新たな変数としたgraphlet correlation distance (gcd: Yaveroglu et al 2014). edge安定性とgraphletの安定性により優れた計算となる
  1. 他の方法もあるが、regularization parameterを決定する必要がある(Liu et al. 2010)。

model selection R package: pulsar（Parallelized Utilities for Lambda Selection Along a Regularization path）
  • 高速化オプション：parallelizations, gcd, natural connectivity
  • モデル学習オプション：Glasso, neighborhood selection, QUIC, clime... etc...

3. rMAGMA package

: Microbial Association Graphical Model Analysis for Rref. GitHub微生物HTSデータのような、noisyな構造のcountデータを扱う：
  • an excess of 0 count
  • overdispersion of sequencing data
  • compositionality
  • possible covariate inclusion

ネットワークの解析


*直接的な相互作用のみを検出する場合は、潜在変数モデル（環境要因から環境選好性を推定）を検討する（Warton et al. 2015)

*Hubに注目（Toju et al. 2017, プレスリリース有り）

*生態学的地位の共有の効果を除外したネットワーク分析として時系列サンプリングがある


モジュラリティ(Q) 

："分割されたコミュニティ内の辺の数とコミュニティ外の辺の数の比較から、高密度のグループを抽出できているかを示す指標"（「ネットワーク分析」）

☞ d <- modularity(DAT)


コミュニティ 

：内部同士の結合が多く、外との結合が比較的少ないネットワークの部分集合（クラスター）

cf. locanな定義, 頂点の類似度に基づく定義(村田2009)

コミュニティ抽出法 ref. Takemoto's Web:

*速さ: fastgreedy.community < Girvan-Newman < 固有ベクトル < spinglass.communiy  *精度: spinglass.communiy > 固有ベクトル > fastgreedy. community > Girvan-Newman

【OUTPUT】


Non-Overlapping法：nodeが単一のcommunityのみに所属

モジュラリティ最大化法（Greedy algorithms）


: 媒介性が高いedgeを削除して分割

☞ d <- cluster_edge_betweenness(DAT)(ref. Workshop Web)

☞ d <- edge.betweenness.community(DAT) 
☞ d <- walktrap.community(DAT, modularity=T) 
:Q値が高くなるようにnodeをまとめる

☞ d <- cluster_fast_greedy(DAT)(ref. Workshop Web)

☞ d <- fastgreedy.community(DAT)
:グラフラプラシアンによりQ値が最大となるような分割を探す☞ d <- leading.eigenvector.community(DAT, options=list(maxiter=1000000,ncv=5)))

☞ cluster_louvain(DAT) (ref. Workshop Web)

☞ multilevel.community(DAT)
:焼きなまし法☞ spinglass.community(DAT) *unconnected graphは不可

Overlapping法：nodeが複数のcommunityに属することを許す

①☞ getLinkCommunities(DAT)

*DAT=as_edgelist(matrix or data frame)ref.  TJO's Web 
gl.lc <- getLinkCommunities(as_edgelist(DAT))plot(gl.lc, type="graph")      # ネットワークplot(gl.lc, type="members")    # 複数コミュニティに属するノード

②☞ getOCG.clusters(DAT)    

*DAT=as_edgelist(matrix or data frame)ref. MSato's web....

モジュラリティ(Q) 

："分割されたコミュニティ内の辺の数とコミュニティ外の辺の数の比較から、高密度のグループを抽出できているかを示す指標"（「ネットワーク分析」）

☞ Q = max(d$modularity)☞ #Q <- max(d$membership)RES <- NULL      #Q値の保存num_try <- 100   #試行回数 for (i in 1:num_try){memb <- cut_at(d,i)    #d=コミュニティ検出結果 RES[i] <- modularity(DAT,memb)}plot(RES, type="b"); which(RES==max(RES))
#-Q値が最大となる分割/community検出（ランダムウォーク法） res.com <– walktrap.community(DAT,modularity=T)memb <- community.to.membership(DAT, res.com$merges, steps=which.max(res.com$modularity)-1)#- コミュニティ毎に色を設定V(DAT)$color <- rainbow(length(memb$csize))[memb$membership+1]#- ネットワーク描画plot(DAT,layout=layout.fruchterman.reingold, edge.arrow.size=0.5)

中心性

: どのnodeが重要かを測る尺度（複数の考え方がある）

Ref.: reruru's web・Workshop Web・Takemoto's Web
:そのnodeの次数（他と繋がるnode）

☞ degree(DAT)/(vcount(DAT)-1)

0~1の値にnormalize次数分布 ☞ plot(degree.distribution(DAT))
:そのnodeを通る測地線の数に基づく（他を仲介するnode）

☞ betweenness(DAT)

:各nodeの平均測地適距離の逆数☞ closeness(DAT) *unconnected graphは不可
：重要な(次数が高い)nodeに繋がっているかを考慮した中心性

☞ evcent(DAT)$vector

☞ center_eigen(DAT, directed=T, normalized=T)$vector [SYN] 

Spiec-Easi including SparCC

ref. GitHub; Kurtz et al. 2015

Installation

事前インストールが必要: 

@MacOS %xcode-select --install@R 
  • pulsar: モデル選択（SparCCに追加。高速モデルあり）
  • huge, devtools etc...

@R% library(devtools)@R% install_github("zdk123/SpiecEasi", ref='dev') # ref.GitHub#36@R% library(SpiecEasi)

Usage

@R:SpiecEasi

data=amgut1.filt（127 OTUs × 289 samples）
Using "Phylotag" package @R:SpiecEasi

Errors in SpiecEasi

@R:SpiecEasi

data=amgut1.filt（127 OTUs × 289 samples）
Using "Phylotag" package @R:SpiecEasi

MBPro intelでは途中で止まる。。。おそらくメモリ不足

➡︎ MBPro M1はいける

--- ERROR!!---> ig2.mb <- adj2igraph(getRefit(se.mb.amgut2),vertex.attr=list(name=taxa_names(amgut2.filt.phy)))
 *** caught segfault ***address 0x1, cause 'memory not mapped'
Traceback: 1: t(object) 2: t(object) 3: .class1(object) 4: as(t(object), "sparseVector") 5: .local(object, ...) 6: Matrix::isSymmetric(adjmatrix) 7: Matrix::isSymmetric(adjmatrix) 8: graph.adjacency.sparse(adjmatrix, mode = mode, weighted = weighted,     diag = diag) 9: igraph::graph.adjacency(Adj, mode = "undirected", weighted = TRUE,     diag = diag)10: adj2igraph(getRefit(se.mb.amgut2), vertex.attr = list(name = taxa_names(amgut2.filt.phy)))
Possible actions:1: abort (with core dump, if enabled)2: normal R exit3: exit R without saving workspace4: exit R saving workspaceSelection: 

Phyloseq

ref. GitHub-Joey711 (Tutrial)

ref. Ushio's blog (基本的使用方法)

ref. 疫学とR（phyloseqのデータ作成）

ref. GitHub

ref. Phyloseq tutorial (Vaulot's GitHub)

ref. Rで学ぶデータサイエンス8（ネットワーク解析）

グラフのsave

ref. ぬいぐるみライフ？

library(ggplot2)ggsave(file="xxxxx.png", plot=plot.obj, dpi=100, width=6.4, height=4.8)*画像形式は拡張子で自動判別

どのように使う？

導入

ref. ：:日本地図を使った解説

igraph on R

Workshop（日本語）：Network analysis with R and igraph: NetSci X Tutorial （日本語版） 

Workshop（英語）： Network analysis with R and igraph: NetSci X Tutorial 

manualに沿った解説：igraphでネットワーク解析（Konishi's） 

Workshopの追加説明（詳細な解説）:  Takemoto's Labo 

「Rで学ぶデータサイエンス–ネットワーク分析–」に沿った実践例: TJO's ブログ 

コミュニティ抽出：Murata 2009

-----------------

.gmlファイルの読み込み

read_graph("./data/network/network__newid.gml", format="gml")


------

plot parameters (ref. Nemoto's web)

ノード

  • vertex.color ノードの色
  • vertex.frame.color ノードの縁の色
  • vertex.shape “none”, “circle”, “square”, “csquare”, “rectangle”, “crectangle”, “vrectangle”, “pie”, “raster”, “sphere”の中から1つ
  • vertex.size ノードの大きさ (デフォルトは15)
  • vertex.size2 ノードの2つ目の大きさ (ノードが長方形の時のもう一方の辺の長さなど)
  • vertex.label ノードにラベルをつけるために用いられる文字列ベクトル
☞ vertex.label=NA（ラベルなし）
  • vertex.label.family ラベルのフォント名 (“Times”, “Helvetica” など)
  • vertex.label.font フォント: 1 プレーン, 2 ボールド, 3, イタリック, 4 ボールドイタリック, 5 シンボル
  • vertex.label.cex フォントのサイズ (倍率)
  • vertex.label.dist ラベルと頂点の距離
  • vertex.label.degree 頂点に対するラベルの位置：0だと右側、“pi”だと左側、“pi/2”だと下、“-pi/2”だと上

エッジ

  • edge.color エッジの色
  • edge.width エッジの太さ (デフォルトは1)
  • edge.arrow.size 矢印のサイズ (デフォルトは1)
  • edge.arrow.width 矢印の太さ (デフォルトは1)
  • edge.lty 線の種類：0 か “blank”, 1か“solid”, 2か“dashed”, 3か“dotted”, 4か“dotdash”, 5か“longdash”, 6か“twodash”
  • edge.label エッジにラベルをつけるために用いられる文字列ベクトル
  • edge.label.family ラベルのフォント名 (ノードに同じ)
  • edge.label.font ノードに同じ
  • edge.label.cex エッジのラベルののフォントサイズ
  • edge.curved エッジの曲がり具合を0-1で表わす(FALSEだと0, TRUEだと0.5に設定される)
  • arrow.mode エッジに矢印があるかどうかを特定するベクトル(0だと矢印なし、1だと後ろ向き、2だと前向き、3だと両方向)

その他

  • margin プロット周囲にある空スペースのマージンで長さ4のベクトル
  • frame もしTRUEなら、プロットは枠付きになる
  • main 設定するとプロットにタイトルが付く
  • sub 設定するとプロットにサブタイトルが付く

微生物群集統計解析

Qiime2を用いた16S rRNAアンプリコン解析 

ANCOM (analysis of composition of microbiomes) (Mandal et al 2015)