前言

我们知道GWAS描述的是snp与表型之间的关系，即利用线性模型寻找与某种性状显著相关联的位点。但是发挥生物学功能的往往是蛋白质，蛋白质是由转录本翻译而来，那么建立基因表达量与表型的关联将会使得分析更进一步

这张图讲述了GWAS与转录组数据整合分析的流程

基于snp和转录组数据（PrediXcan）

可参考文章A gene-based association method for mapping traits using reference transcriptome data
首先我们先从原理上介绍下这款软件的工作流程

第一步
我们对测序数据 call 出来snp以后会得到一张表格：

显然 C 突变为 T （C/T为杂合突变，T/T为纯和突变）会促使表型值上升，与表型值成正相关

那么，判断突变与表型值呈正相关还是负相关，我们可以设置一个对照（比方说设置纯和未突变的为对照）看看突变是否会引起表型值的上升

并且测序数据 call 出来的snp信息可以和表型值相关联

而我们RNA-seq的表格为：

第二步（可选，若用户没用基因表达矩阵执行此步骤）
利用snp数据和转录组各个基因表达量的数据，建立它们之间的线性的权重关系

其中，w_k,g 代表第 k 个snp与第 g 个基因表达量之间的权重（权重 w_k,g 是作者利用机器学习的思想，利用GTEx Project, GEUVADIS 和 DGN数据库中基因型数据和基因表达数据做训练集，即利用已知的基因表达量 T_g 和 snp 的基因型数据 X_k 通过 LASSO 和 elastic net 来计算权重 w_k,g）；T_g 代表第 g 个基因的表达量；X_k 代表第 k 个snp的基因型（因子型变量转换为[数值型变量] 0，1，2）

因此如果你没有现成的转录组数据，你可以利用作者已开发好的模型（利用机器学习的方法开发了部分模型）。作者已经利用已发的数据做了模型训练了，如果你有对应组织的snp数据，可以到PredictDB下载对应的模型，用于预测基因的表达量信息。

第三步
关于权重 w_k,g 的计算可以利用LASSO和 elastic net 来计算，由此可知，对于某一个 Tissue 来说，第 g 个基因的表达量可以用snp来线性表示，那么结合trait的值（我们在这里称为trait的表达量），再次拟合一个线性模型：

其中 Tg 第 g 个基因的表达量；γ 代表回归系数；Y 代表表型值（表型表达量））；Y_n = γ₁×T_1,n+ γ₂×T_2,n + ... + γ_m×T_m,n，n代表第n个id，m代表m个基因
这样就可以将基因表达量与表型值联系起来了，其本质就是基因表达量和表型直接的关系

对于该模型，我们可以这样理解，对于每一个基因 g 来说，在各个sample（id）中的表达量不同，而每一个sample（id）的表型值也不同，因此可以建立基因 g 在不同sample（id）中的表达量与在sample（id）的表型值之间的线性关系（如上图）
其中，每个点对应不同的sample（id）；每个点对应的横坐标为基因 g 在不同sample中的表达量；每个点对应的纵坐标为不同sample对应的表型值

该表代表基因与表型的关系，那么回归系数 γ 的为正，那么代表基因表达量越高，则对性状的影响成正相关；反之为负，则代表基因表达量越高，则对性状的影响成负相关，后面的pval为回归系数的显著性

注: PrediXcan的使用

运行PrediXcan需要输入三个文件：转录组表达矩阵，基因型文件和样本信息文件：

1. 基因型文件：该文件每一行表示一个SNP，包含的信息分别为：chromosome rsid position allele1 allele2 MAF，后面的每一列的内容是每一个样本在该SNP allele2的dosage，最好是每一条染色体分开制作文件。
2. 样本信息文件：直接将PLINK的fam文件导入即可。
3. 基因表达矩阵

可选，利用已有的模型预测一个基因表达矩阵

./PrediXcan.py 
--predict 
--dosages genotype/ 
--dosages_prefix chr 
--samples samples.txt 
--weights model/DGN-HapMap-2015/DGN-WB_0.5.db # 已存在的模型
--output_prefix results/DGN-HapMap

建立基因表达与表型的关系

./PrediXcan.py 
--assoc 
--pheno My_pheno.txt 
--mpheno 1 
--pred_exp results/TW_Brain_Frontal_predicted_expression.txt #预测的基因表达矩阵，作者也可自行提供
--logistic 
--output_prefix results/DGN-HapMap

基于GWAS-summary数据（Summary-PrediXcan）

可参考文章Exploring the phenotypic consequences of tissue specific gene expression variation inferred from GWAS summary statistics

这种方法不依赖于call出来的snp和转录组数据，而是直接利用GWAS-summary数据来建立基因表达与表型值之间的关系
我们知道，GWAS-summary描述的是不同的snp（基因型）和表型值之间的关系

其中，wIg 代表第 I 个snp对gene g表达量产生影响的权重（该权重即snp I 是否突变对gene g表达量产生变化的回归系数）；βl 代表第 l 个snp 对表型值的影响的回归系数（也称为效应值）；se(βl) 代表所有回归系数（效应值）的标准误；σl 代表所有回归系数（效应值）的标准差；σg 代表基因 g 在各个样本中表达量的标准差；Zg 即为基因 g 对表型值的回归系数（效应值）

其中，wIg 是已经训练好的snp与基因 g 表达量之间的权重

因此 Zg 为基因 g 表达量和表型值之间线性模型的回归系数，越大说明对表型的正向影响越大；反之越小代表对表型的负向影响越大