筆記內(nèi)容

• 以iris(鳶尾花)數(shù)據(jù)集為例，tmap各步驟input, output,結(jié)果解讀便锨，參數(shù)調(diào)整和選取澜倦。包括以下腳本：
  envfit_analysis.py
  Network_generator.py
  SAFE_analysis.py
  SAFE_visualization.py
• 不涉及具體的算法解讀，僅為走一遍數(shù)據(jù)分析流程做參考
• 每個(gè)可執(zhí)行的腳本(.py)都說明其input, output及目的
• 小結(jié)及注意

tmap documents
tmap github

envfit_analysis.py

目的：

1. 用到rpy2嘹叫，使用R的vegan包中的方法做envfit嫩絮。tmap的算法并不涉及envfit丛肢，該結(jié)果僅用于與tmap結(jié)果比較。如果不需要可以用--dont_analysis跳過剿干。
2. 處理input的metadata.
   默認(rèn)如果metadata中的某個(gè)變量(column)缺失值占比超過60%蜂怎，則刪掉該變量；對于類別變量(category variable)做one-hot處理;對連續(xù)型變量用中位數(shù)填充缺失值置尔。
   另外處理metadata的函數(shù)為tmap.api.general的process_metadata_beta
   inputdata(即-I后輸入的文件)在這里沒有處理派敷，只是讀入，并且確保其與metadata所含有的樣本(row)一致。

envfit_analysis.py  -I iris.csv 
                    -M iris_meta.csv 
                    -O iris_envfit.csv 
                    -tn 'iris'  # 這里設(shè)置為iris,則Ouput文件名均以iris開頭
                    --dont_analysis # 設(shè)置則不做envfit分析篮愉，僅處理metadata和data
                    --keep  # 保留output中處理后的.data, .dis及.metadata文件
                    -v

input:

# load并head一下iris.csv: 
sepal length (cm)   sepal width (cm)    petal length (cm)   petal width (cm)
0     5.1           3.5                 1.4                 0.2
1     4.9           3.0                 1.4                 0.2
2     4.7           3.2                 1.3                 0.2
3     4.6           3.1                 1.5                 0.2
4     5.0           3.6                 1.4                 0.2

# load并head一下iris_meta.csv
    type
0   setosa
1   setosa
2   setosa
3   setosa
4   setosa

output:

Network_generator.py

目的:

根據(jù)inputdata(這里為iris.csv)生成網(wǎng)絡(luò)圖. 下圖為tmap生成網(wǎng)絡(luò)圖的簡要步驟，參數(shù)設(shè)置影響著網(wǎng)絡(luò)圖對刻畫和提取高維數(shù)據(jù)特征的精細(xì)程度颠蕴。

參數(shù)選擇的文檔參考： https://tmap.readthedocs.io/en/latest/param.html
先Network_generator.py -h泣刹，配合上述文檔看看參數(shù)說明。

Network_generator.py -I iris.csv 
                     -O iris.graph  # 生成Python3 pickle格式的graph
                     -f 'PCOA'      # 選擇降維方式犀被，default為PCOA
                     -ol 0.75       # 具體見下圖
                     -eps 95
                     -ms 3
                     -v

# Network_generator.py得到的是一個(gè).pickle的文件椅您，通過quick_vs.py迅速看一下graph長什么樣
quick_vis.py -G iris.graph 
             -O iris.html 
             -M iris_meta.csv  # 可以選擇性的傳個(gè)metadata, 看看某個(gè)變量在網(wǎng)絡(luò)上的分布如何
             -col 'type'       # 指定某個(gè)變量
             -t 'categorical'  # 指定的變量是數(shù)值型(numerical)還是字符型(categorical)

input:

iris.csv (同上)

output:

iris.graph 及其log output，包含了graph的參數(shù)信息寡键，生成了多少node和edge等掀泳。可以保存下來留做參考或debug.

Graph
Contains 91 nodes and 133 samples  
During constructing graph, 17 (88.67%) samples lost # 這里是指有88.67%的samples被留了下來西轩，17個(gè)samples被去掉了员舵。
Used params:
cluster params
n_jobs: 1
leaf_size: 30
metric: euclidean
metric_params: None
algorithm: auto
min_samples: 3
p: None
eps: 0.475659777979
=================
cover params
r: 20
overlap: 0.75
=================
lens params
lens_0:
metric: precomputed
components: [0, 1]

對網(wǎng)絡(luò)圖的著色：顏色代表的是原始數(shù)值，而不是網(wǎng)絡(luò)富集的統(tǒng)計(jì)學(xué)顯著性(SAFE score)藕畔。對類別變量則為1和0兩種著色马僻，對連續(xù)型變量則將其原始數(shù)值map到每個(gè)點(diǎn)上。

SAFE_analysis.py

目的

計(jì)算每個(gè)node的SAFE (spatial analysis of functional enrichment) score

                 ↓↓↓注意這里要指定calculating mode
SAFE_analysis.py both -G iris.graph 
                      -M temp.envfit.metadata temp.envfit.data  # 可以輸入多個(gè)metadata文件注服，包括構(gòu)建.graph的data文件
                      -P SAFE    # 輸出文件的開頭字符
                      -i 1000    # permutation的次數(shù)
                      -p 0.05    # 計(jì)算顯著性的cutoff value
                      --raw      # 最好保留中間文件
                      -v

input:

iris.graph： 即上一步得到的.graph文件
temp.envfit.metadata, temp.envfit.data : -M 可以輸入多個(gè)metadata文件

output:

import pickle

# 看看SAFE_raw_coldict
coldict = pickle.load(open('SAFE_raw_coldict', 'rb'))
print(coldict.keys())  # dict_keys('iris.envfit.metadata', 'iris.envfit.data')
print(coldict['iris.envfit.metadata']) # ['type_setosa', 'type_versicolor', 'type_virginica']

# 看看SAFE_raw_enrich
enrich = pickle.load(open('SAFE_raw_enrich','rb'))
print(enrich.keys())   # dict_keys(['params', 'data'])
enrich['data'].head() 
#  一共91個(gè)node, 于是該dataframe有91行，對應(yīng)每個(gè)node,每個(gè)feature的SAFE score
   type_setosa  type_versicolor  type_virginica  sepal length (cm)  sepal width (cm)  petal length (cm)  petal width (cm)
0     0.829397             -0.0            -0.0               -0.0          -0.00000               -0.0              -0.0
1     0.829397             -0.0            -0.0               -0.0          -0.00000               -0.0              -0.0
2     0.829397             -0.0            -0.0               -0.0           0.11586               -0.0              -0.0
3     0.829397             -0.0            -0.0               -0.0           0.30695               -0.0              -0.0
4     0.829397             -0.0            -0.0               -0.0           0.30695               -0.0              -0.0

下圖為SAFE_iris.envfit.data_enrich.csv及其解讀：

SAFE_iris.envfit.data_enrich.csv

SAFE_visualization.py

目的:

基于SAFE_analysis.py得到的SAFE score及其顯著性開展的后續(xù)分析袜茧。主要實(shí)現(xiàn)三個(gè)分析目的：

三個(gè)分析目的實(shí)現(xiàn)示意圖

# ranking: 當(dāng)然你也可以直接就用 SAFE_analysis.py中得到的SAFE summary表格自行畫圖瞎颗，比方說使用SAFE_iris.envfit.data_enrich.csv
SAFE_visualization.py ranking -G iris.graph 
                              -S2 SAFE_iris.envfit.metadata_enrich.csv  # 這里也可以再加一個(gè)iris_envfit.csv, 和envfit做比較件甥。不加則如下圖的ranking output 
                              -O ranking.html  

# stratification: 顯示所有的顯著富集feature
SAFE_visualization.py stratification -G iris.graph
                                     -S1 SAFE_raw_enrich  # 每個(gè)node的SAFE socre
                                     -O strtification.html

# stratification: 指定某一個(gè)顯著富集的feature
SAFE_visualization.py stratification -G iris.graph
                                     -S1 SAFE_raw_enrich  
                                     --col 'petal length(cm)'
                                     # --allnodes 則顯示所有的點(diǎn)對該feature的富集情況
                                     -O strtification.html

# stratification指定兩個(gè)或多個(gè)feature，則--col 'XXX' 'XXXX'哼拔，不需要加--allnodes

# ordination: -S2 可以放metadata和data兩個(gè)文件
SAFE_visualization.py ordination -S1 iris_raw_enrich 
                                 -S2 SAFE_iris.envfit.data_enrich.csv SAFE_iris.envfit.metadata_enrich.csv 
                                 -O ordination.html

input:

如代碼中所示引有。

output:

ranking.html

stratification.html node的大小表示其含有Samples的個(gè)數(shù)，legend中括號(hào)內(nèi)數(shù)字為enriched的點(diǎn)個(gè)數(shù)

stratification_petal_length.html 右圖為加了--allnodes的情況

ordination.html

小結(jié)及注意

1. Network_generator.py構(gòu)建了整個(gè)tmap分析中最基礎(chǔ)的網(wǎng)絡(luò)結(jié)構(gòu)倦逐。在后續(xù)的分析中不會(huì)再改變這個(gè)基本結(jié)構(gòu)譬正。如果畫出了很奇怪的network，建議根據(jù)數(shù)據(jù)本身的情況調(diào)整參數(shù)-r, -f等僻孝。網(wǎng)絡(luò)圖的構(gòu)建參數(shù)的選擇具有經(jīng)驗(yàn)性导帝，多多嘗試 (-_-)
2. 構(gòu)建好了網(wǎng)絡(luò)結(jié)構(gòu)，網(wǎng)絡(luò)上每一個(gè)node都是一簇基于microbiome profile(如果你的input是一個(gè)OTU table之類的豐度表)相似的樣本穿铆，node之間有edge相連您单，則表明兩個(gè)node之間有共有樣本。
3. 對于網(wǎng)絡(luò)結(jié)構(gòu)上的每一個(gè)node荞雏，都可以計(jì)算某一個(gè)feature的SAFE score虐秦。這個(gè)SAFE score量化了一種“富集程度的統(tǒng)計(jì)學(xué)顯著性”：對該feature而言，該node周圍的node凤优，有該feature值偏高的富集趨勢悦陋。通過隨機(jī)混洗(permutation)打亂node在network上的位置，如果發(fā)現(xiàn)這種富集趨勢顯著不同于隨機(jī)結(jié)果筑辨，則說明這種富集是有統(tǒng)計(jì)學(xué)意義的俺驶。
4. SAFE_analysis.py有both/enrich/decline三種SAFE score計(jì)算模式可以選擇。enrich即第三條中說的富集程度棍辕，decline即“該node周圍的node, 有該feature值偏低的富集趨勢”暮现。both則將兩種模式都做了。對于微生物組研究楚昭，可以只選擇enrich的方式栖袋。因?yàn)槟硞€(gè)物種的豐度升高可能更加有意義，更能夠解讀生物學(xué)意義抚太。
5. 另外.graph中有茫茫多的屬性塘幅，用pickle讀入.graph之后昔案，graph.__dir__()查看所有屬性。這里例舉幾個(gè)常用的：

目前代碼仍在不斷完善中匾乓，可能存在各種小bug呼伸，如有問題請大家及時(shí)同我們聯(lián)系，幫助我們改進(jìn)代碼 orz