tmap数据分析基本步骤和结果解析
2020-01-16 本文已影响0人
GPZ_Lab
笔记内容
- 以iris(鸢尾花)数据集为例，tmap各步骤input, output,结果解读，参数调整和选取。包括以下脚本：
envfit_analysis.py
Network_generator.py
SAFE_analysis.py
SAFE_visualization.py- 不涉及具体的算法解读，仅为走一遍数据分析流程做参考
- 每个可执行的脚本(.py)都说明其input, output及目的
- 小结及注意
envfit_analysis.py
目的：
- 用到rpy2，使用R的vegan包中的方法做envfit。tmap的算法并不涉及envfit，该结果仅用于与tmap结果比较。如果不需要可以用
--dont_analysis跳过。 - 处理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:
| file | what's this |
|---|---|
| iris_envfit.csv | 如果设置了--dont_analysis则没有输出R的envfit分析结果 |
| iris.envfit.data | iris.csv文件 |
| iris.envfit.dis | 样本*样本的距离矩阵，可用-m设置哪种距离矩阵，默认为braycurtis
|
| iris.envfit.metadata | 处理后的metadata文件，字符变量one-hot处理，数值变量清洗并填充空值 |
| iris.envfit.raw.metadata | 原metadata |
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:
| file | what's this |
|---|---|
| SAFE_raw_coldict | 保存了iris.envfit.data, iris.envfti.metadata的column信息。import pickle后，用pickle.load(open('SAFE_raw_coldict', 'rb'))在Python中打开看看 |
| SAFE_raw_decline | 保存了decline mode的SAFE SCORE及其参数。SAFE SCORE以dataframe形式储存，每个node的每个feature对应着一个SAFE SCORE。参考下面的code. |
| SAFE_raw_enrich | 保存了enrich mode的SAFE SCORE及其参数。同上。 |
| SAFE_iris.envfit.data_enrich.csv | 在enrich mode下，row均为iris.envfit.data的column name, 有6个column, 为根据SAFE score做的一系列summary, 详见下面的code |
| SAFE_iris.envfit.data_decline.csv | decline mode版本的 SAFE score summary，同上 |
| SAFE_iris.envfit.metadata_enrich.csv | row均为iris.envfit.metadata的column name,同上 |
| SAFE_iris.envfit.metadata_decline.csv | row均为iris.envfit.metadata的column name,同上 |
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及其显著性开展的后续分析。主要实现三个分析目的：
三个分析目的实现示意图
| 分析目的 | how | 意义 |
|---|---|---|
| ranking | 即上表中SAFE enriched score: 对每个feature，求其所有具有显著性node SAFE score之和后排序 | 与linear regression中的effective size类似，为feature的重要性排序。ranking中SAFE score之和越大，则其对微生物组变化的贡献越大 |
| stratification | stratification是挑选出每个node SAFE score最大的feature，且通过permutation验证统计学显著性，将每个Node的enriched feature着色 | ranking是从数据整体出发，看最重要的feature, stratification则是从每个node出发，找到每个node显著富集的feature。 一簇Node如果相互连接，且enriched feature一致，则可以通过node找到该feature富集的samples (subgroup). |
| ordination | 对SAFE score做PCoA | PCoA图中的点为feature, 互相靠近的feature则：在基于微生物profile的情况下，富集趋势趋于一致 |
# 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
小结及注意
-
Network_generator.py构建了整个tmap分析中最基础的网络结构。在后续的分析中不会再改变这个基本结构。如果画出了很奇怪的network，建议根据数据本身的情况调整参数-r,-f等。网络图的构建参数的选择具有经验性，多多尝试 (-_-) - 构建好了网络结构，网络上每一个node都是一簇基于microbiome profile(如果你的input是一个OTU table之类的丰度表)相似的样本，node之间有edge相连，则表明两个node之间有共有样本。
- 对于网络结构上的每一个node，都可以计算某一个feature的SAFE score。这个SAFE score量化了一种“富集程度的统计学显著性”：对该feature而言，该node周围的node，有该feature值偏高的富集趋势。通过随机混洗(permutation)打乱node在network上的位置，如果发现这种富集趋势显著不同于随机结果，则说明这种富集是有统计学意义的。
-
SAFE_analysis.py有both/enrich/decline三种SAFE score计算模式可以选择。enrich即第三条中说的富集程度，decline即“该node周围的node, 有该feature值偏低的富集趋势”。both则将两种模式都做了。对于微生物组研究，可以只选择enrich的方式。因为某个物种的丰度升高可能更加有意义，更能够解读生物学意义。 - 另外
.graph中有茫茫多的属性，用pickle读入.graph之后，graph.__dir__()查看所有属性。这里例举几个常用的：
code(读入的graph用G表示) |
output |
|---|---|
| G.info | Network_generator.py的输出，快速了解graph的基本情况 |
| G.nodes | 所有node的ID，是networkx的一个类, 用G.nodes.data()看一看 |
| G.edges | 所有的edges, 用G.edges.data()看一看 |
| G.sample_names | 所有node中的samples |
| G.data | 所有ndoe的坐标 |
| G.size | ouput为一个dictionary,key为nodeID，value为里面包含的samples个数 |
| G.node2sample(nodeID) | 需要输入一个参数nodeID, output为该node中包含的样本ID |
| G.sample2node(sampleID) | 需要输入一个参数sampleID, output为该sample出现的nodeID |
| G.show() | 将网络图简单plot |
目前代码仍在不断完善中，可能存在各种小bug，如有问题请大家及时同我们联系，帮助我们改进代码 orz
