笔记内容

• 以iris(鸢尾花)数据集为例，tmap各步骤input, output,结果解读，参数调整和选取。包括以下脚本：
  envfit_analysis.py
  Network_generator.py
  SAFE_analysis.py
  SAFE_visualization.py
• 不涉及具体的算法解读，仅为走一遍数据分析流程做参考
• 每个可执行的脚本(.py)都说明其input, output及目的
• 小结及注意

tmap documents
tmap github

envfit_analysis.py

目的：

1. 用到rpy2，使用R的vegan包中的方法做envfit。tmap的算法并不涉及envfit，该结果仅用于与tmap结果比较。如果不需要可以用--dont_analysis跳过。
2. 处理input的metadata.
   默认如果metadata中的某个变量(column)缺失值占比超过60%，则删掉该变量；对于类别变量(category variable)做one-hot处理;对连续型变量用中位数填充缺失值。
   另外处理metadata的函数为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'  # 这里设置为iris,则Ouput文件名均以iris开头
                    --dont_analysis # 设置则不做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

目的:

根据inputdata(这里为iris.csv)生成网络图. 下图为tmap生成网络图的简要步骤，参数设置影响着网络图对刻画和提取高维数据特征的精细程度。

参数选择的文档参考： https://tmap.readthedocs.io/en/latest/param.html
先Network_generator.py -h，配合上述文档看看参数说明。

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得到的是一个.pickle的文件，通过quick_vs.py迅速看一下graph长什么样
quick_vis.py -G iris.graph 
             -O iris.html 
             -M iris_meta.csv  # 可以选择性的传个metadata, 看看某个变量在网络上的分布如何
             -col 'type'       # 指定某个变量
             -t 'categorical'  # 指定的变量是数值型(numerical)还是字符型(categorical)

input:

iris.csv (同上)

output:

iris.graph 及其log output，包含了graph的参数信息，生成了多少node和edge等。可以保存下来留做参考或debug.

Graph
Contains 91 nodes and 133 samples  
During constructing graph, 17 (88.67%) samples lost # 这里是指有88.67%的samples被留了下来，17个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]

对网络图的着色：颜色代表的是原始数值，而不是网络富集的统计学显著性(SAFE score)。对类别变量则为1和0两种着色，对连续型变量则将其原始数值map到每个点上。

SAFE_analysis.py

目的

计算每个node的SAFE (spatial analysis of functional enrichment) score

                 ↓↓↓注意这里要指定calculating mode
SAFE_analysis.py both -G iris.graph 
                      -M temp.envfit.metadata temp.envfit.data  # 可以输入多个metadata文件，包括构建.graph的data文件
                      -P SAFE    # 输出文件的开头字符
                      -i 1000    # permutation的次数
                      -p 0.05    # 计算显著性的cutoff value
                      --raw      # 最好保留中间文件
                      -v

input:

iris.graph： 即上一步得到的.graph文件
temp.envfit.metadata, temp.envfit.data : -M 可以输入多个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个node, 于是该dataframe有91行，对应每个node,每个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及其显著性开展的后续分析。主要实现三个分析目的：

三个分析目的实现示意图

# ranking: 当然你也可以直接就用 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  # 这里也可以再加一个iris_envfit.csv, 和envfit做比较。不加则如下图的ranking output 
                              -O ranking.html  

# stratification: 显示所有的显著富集feature
SAFE_visualization.py stratification -G iris.graph
                                     -S1 SAFE_raw_enrich  # 每个node的SAFE socre
                                     -O strtification.html

# stratification: 指定某一个显著富集的feature
SAFE_visualization.py stratification -G iris.graph
                                     -S1 SAFE_raw_enrich  
                                     --col 'petal length(cm)'
                                     # --allnodes 则显示所有的点对该feature的富集情况
                                     -O strtification.html

# stratification指定两个或多个feature，则--col 'XXX' 'XXXX'，不需要加--allnodes

# ordination: -S2 可以放metadata和data两个文件
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的个数，legend中括号内数字为enriched的点个数

stratification_petal_length.html 右图为加了--allnodes的情况

ordination.html

小结及注意

1. Network_generator.py构建了整个tmap分析中最基础的网络结构。在后续的分析中不会再改变这个基本结构。如果画出了很奇怪的network，建议根据数据本身的情况调整参数-r, -f等。网络图的构建参数的选择具有经验性，多多尝试 (-_-)
2. 构建好了网络结构，网络上每一个node都是一簇基于microbiome profile(如果你的input是一个OTU table之类的丰度表)相似的样本，node之间有edge相连，则表明两个node之间有共有样本。
3. 对于网络结构上的每一个node，都可以计算某一个feature的SAFE score。这个SAFE score量化了一种“富集程度的统计学显著性”：对该feature而言，该node周围的node，有该feature值偏高的富集趋势。通过随机混洗(permutation)打乱node在network上的位置，如果发现这种富集趋势显著不同于随机结果，则说明这种富集是有统计学意义的。
4. SAFE_analysis.py有both/enrich/decline三种SAFE score计算模式可以选择。enrich即第三条中说的富集程度，decline即“该node周围的node, 有该feature值偏低的富集趋势”。both则将两种模式都做了。对于微生物组研究，可以只选择enrich的方式。因为某个物种的丰度升高可能更加有意义，更能够解读生物学意义。
5. 另外.graph中有茫茫多的属性，用pickle读入.graph之后，graph.__dir__()查看所有属性。这里例举几个常用的：

目前代码仍在不断完善中，可能存在各种小bug，如有问题请大家及时同我们联系，帮助我们改进代码 orz