N-of-1 基因异常值检测

用于RNA-seq 表达数据

此软件包通过构建一个从N-of-1样本获得信息的背景数据集的共识分布来识别基因表达数据的异常值。有关更多信息，请参阅模型解释，我们的JCO论文或我们的预印本。

该模型已进行了显著的性能改进，并且与论文中的模型不再相同，但产生的结果基本上是相同的。另一个优点是现在每个基因/数据集组合的先验超参数直接在模型中共享，而以前由于运行时间它们是被近似的。

此工作流程以基因表达数据作为输入，并输出以下内容

SAMPLE_UUID
├── _info
│   ├── _gelman-rubin.tsv
│   ├── _pearson_correlations.txt
│   ├── _pval_runs.tsv
│   └── _run_info.tsv
├── model.pkl
├── pvals.tsv
├── ranks.tsv
├── traceplot.png
├── weights.png
└── weights.tsv

• _info 子目录包含有关模型运行的二级信息
  • _gelman-rubin.tsv - 包含每个模型参数的 Gelman-Rubin诊断 的TSV，包括参数的中位数，应约为 ~1.0。
  • _pearson_correlations.txt - 基因p值之间皮尔逊相关性的单列列表（除非提供了 -d）
  • _pval_runs.tsv - 随着添加背景数据集的基因p值表
  • _run_info.tsv - 软件参数的TSV，用于可重复性和模型运行时间
• model.pkl — 包含PyMC3 model和trace的python pickle。模型必须使用--save-model标志运行，并且可以通过

import pickle
with open(pkl_path, 'rb') as buff:
    data = pickle.load(buff)
model, trace = data['model'], data['trace']

• pvals.tsv — 模型训练的所有基因的p值
• ranks.tsv — 所有组按成对欧几里得距离测量的中位数排名
• traceplot.png — 来自PyMC3线性模型系数和模型错误的迹图
• weights.png — 所有背景数据集模型权重的箱线图
• weights.tsv — 模型权重与背景数据集相关的平均值和标准差

快速入门

1. 安装

pip install --pre gene-outlier-detection

2. 下载先决条件 输入
3. 运行模型

outlier-detection --sample /data/tumor.hd5 \
        --background /data/gtex.hd5 \
        --name TCGA-OR-A5KV-01 \
        --gene-list /data/drug-genes.txt \
        --col-skip 5

依赖关系和安装

此工作流程已在 ubuntu 18.04 和 Mac OSX 上进行测试，但应也能在其他基于 Unix 的系统或 Anaconda 安装下运行。

1. Python 3.6
2. 如果使用 Docker 版本或 Toil 工作流程版本，请使用 Docker
3. HDF5 库（如果输入为 HDF5 格式）
   1. conda install -c anaconda hdf5
4. PyMC3 的 Theano 的 C++ / GCC 编译器
   1. 您有几个选项
      1. conda install theano
      2. apt-get update && apt-get install -y libhdf5-serial-dev build-essential gcc

您可能需要修改您的 ~/.theanorc 以支持更大的括号深度，以便为此模型。

[gcc]
cxxflags = -fbracket-depth=1024

模型说明

以下板图中的 G 代表基因，D 代表数据集（背景数据集）。

我们方法的核心是一个针对 N-of-1 样本基因表达的贝叶斯统计模型。该模型隐式假设样本的基因表达可以被近似为背景数据集基因表达的凸混合。这种混合的系数在基因之间是共享的，类似于每个数据点都是背景数据集中一个基因的表达向量的一线性模型。此外，我们还将每个基因来自每个背景数据集的表达建模为随机变量本身。这使我们能够将来自某些背景集小样本量的统计不确定性直接纳入模型，而不会通过合并淹没任何背景集的贡献

可以使用马尔可夫链蒙特卡罗（MCMC）方法探索模型，以获得近似其后验分布的 y（每个基因）样本。如果我们有一个感兴趣的基因的观察到的表达值（来自 N-of-1 癌症样本），我们可以将其与样本值进行比较。大于（或小于）观察值的样本比例是该表达值的后验预测 p 值的估计。后验预测 p 值可以看作是度量与比较集的期望相比，该表达值是多大程度的异常值。

通过利用 PyMC3 的向量化方法并替换有问题的学生-T 抽样为正态分布，该模型在速度上得到了大幅提升。另一个好处是，不再需要预先拟合学生-T 分布的计算“技巧”，因为背景中的表达值直接在两个 R.V. 之间共享，这两个 R.V. 建模每个基因/数据集组合的表达。下面是一个包含 125 个基因、3 个数据集和背景数据集中 48,500 个值的新的模型图

定义输入

模型需要两个 样本 by 基因 矩阵，一个包含 N-of-1 样本，另一个包含用作背景比较集的样本。它们 应 包含相同的基因集，并且背景数据集必须至少包含一个具有区分组的标签（例如组织、亚型、实验等）的元数据列。元数据列需要在数据框中的 最左侧 列，随后是所有基因。这允许一个简单的启发式方法（--col-skip）指定基因开始的位置，而无需用户提供额外的输入文件。

在这里，tissue 是用于区分样本的一个组。我们会用 --col-skip=5 运行，因为有 5 个元数据列在基因开始之前。

Docker 容器

可以使用以下方式执行包含程序的 Docker 容器

docker run --rm -v $(pwd):/data jvivian/gene-outlier-detection \
        outlier-model \
        --sample /data/inputs/tumor.hd5 \
        --background /data/inputs/gtex.hd5 \
        --name=TCGA-OR-A5KV-01 \
        --gene-list /data/inputs/drug-genes.txt \
        --out-dir /data/outputs/ \
        --col-skip=5

参数

运行程序时使用的参数的解释。

• --sample
  • 包含样本的矩阵（.tsv / .hd5）的路径。
• --background
  • 背景数据集的矩阵路径，至少包含一个用于标记不同背景组的标签向量的列。
• --name
  • 样本矩阵中对应要运行的样本的行名。
• --out-dir
  • 输出目录
• --gene-list
  • 用于模型训练和计算p值的基因的单列文件。在100-200个基因后，最好将基因分成批次并行运行。
• --col-skip
  • 在背景矩阵中跳过的元数据列数。此数字之后的列应该是带有表达值的基因。
• --group
  • 背景矩阵中区分不同背景数据集的分类列向量名称。
• --num-backgrounds
  • 模型训练中包含的最大背景分类组数。模型将从单个背景数据集开始运行，并迭代添加更多数据集，直到p值收敛或达到num-backgrounds。
• --max-genes
  • 运行的基因最大数量。即如果提供了基因列表，将使用ANOVA进行额外基因的分析。如果基因列表中没有足够组织特异性基因，则此操作有助于改进beta系数。建议将模型的max-genes设置为比--gene-list中存在的基因至少多10-20个。
• --pval-convergence-cutoff
  • 停止添加额外背景数据集的p值Pearson相关系数截止值。
• --num-training-genesgi
  • 如果基因列表为空，将使用ANOVA选择基因集。通常在测试之外不常用。
• --tune
  • 开始MCMC采样过程的调整步骤数。默认值为500，但在模型难以收敛的情况下，750或1000可能更有用。
• --disable-iter
  • 此标志禁用迭代运行，并使用--num-backgrounds运行单个模型。
• --save-model
  • 此标志将保存模型/跟踪的序列化版本。对于调试或检查所有模型参数非常有用。

Toil版本的工作流程

可在此处找到Toil版本的工作流程here。这允许在集群或云计算集群上以规模运行多个样本，但需要Python 2.7和pip install pandas toil==3.19.0

Toil现在兼容Python 3，因此建议在Toil代码更新之前不使用此方法运行程序。

引用

如果您在研究中利用了此方法，请引用我们的JCO论文。

Vivian, John, et al. "Bayesian Framework for Detecting 
Gene Expression Outliers in Individual Samples." 
JCO Clinical Cancer Informatics 4 (2020): 160-170.