前言
我们知道GWAS描述的是snp与表型之间的关系，即利用线性模型寻找与某种性状显著相关联的位点。但是发挥生物学功能的往往是蛋白质，蛋白质是由转录本翻译而来，那么建立基因表达量与表型的关联将会使得分析更进一步
这张图讲述了GWAS与转录组数据整合分析的流程
基于snp和转录组数据（PrediXcan）
可参考文章A gene-based association method for mapping traits using reference transcriptome data
首先我们先从原理上介绍下这款软件的工作流程
第一步
我们对测序数据 call 出来snp以后会得到一张表格：
其中 id 代表每一个个体sample，rs 代表基因组上被call出来的snp（每一列的 rs 代表同一位置的snp），用连续型变量表示（假设二倍体参考基因组上某位点为 C/C，那么C/C = 0，C/T = 1，T/T = 2，这样转换的目的是将因子型变量（不同的基因型）转换成数值型（连续型）变量，方便利用线性模型建模）；基于上表的结果，每个snp都有三种基因型（用0，1，2表示），而每个snp的每一种基因型又对应着不同的表型值，因此我们就可以建模看一下每一个snp不同的基因型会对表型值有什么样的影响
显然 C 突变为 T （C/T为杂合突变，T/T为纯和突变）会促使表型值上升，与表型值成正相关
那么，判断突变与表型值呈正相关还是负相关，我们可以设置一个对照（比方说设置纯和未突变的为对照）看看突变是否会引起表型值的上升
并且测序数据 call 出来的snp信息可以和表型值相关联
而我们RNA-seq的表格为：
其中 id 代表每一个个体sample，Tissue 代表不同的组织，每一个 Tissue 对于一张表达矩阵，g 代表每一个基因
第二步（可选，若用户没用基因表达矩阵执行此步骤）
利用snp数据和转录组各个基因表达量的数据，建立它们之间的线性的权重关系
建立关系公式如下：
其中，wk,g 代表第 k 个snp与第 g 个基因表达量之间的权重（权重 wk,g 是作者利用机器学习的思想，利用GTEx Project, GEUVADIS 和 DGN数据库中基因型数据和基因表达数据做训练集，即利用已知的基因表达量 Tg 和 snp 的基因型数据 Xk 通过 LASSO 和 elastic net 来计算权重 wk,g）；Tg 代表第 g 个基因的表达量；Xk 代表第 k 个snp的基因型（因子型变量转换为[数值型变量] 0，1，2）
因此如果你没有现成的转录组数据，你可以利用作者已开发好的模型（利用机器学习的方法开发了部分模型）。作者已经利用已发的数据做了模型训练了，如果你有对应组织的snp数据，可以到PredictDB下载对应的模型，用于预测基因的表达量信息。
第三步
关于权重 wk,g 的计算可以利用LASSO和 elastic net 来计算，由此可知，对于某一个 Tissue 来说，第 g 个基因的表达量可以用snp来线性表示，那么结合trait的值（我们在这里称为trait的表达量），再次拟合一个线性模型：
其中 Tg 第 g 个基因的表达量；γ 代表回归系数；Y 代表表型值（表型表达量））；Yn = γ1×T1,n+ γ2×T2,n + ... + γm×Tm,n，n代表第n个id，m代表m个基因
这样就可以将基因表达量与表型值联系起来了，其本质就是基因表达量和表型直接的关系
对于该模型，我们可以这样理解，对于每一个基因 g 来说，在各个sample（id）中的表达量不同，而每一个sample（id）的表型值也不同，因此可以建立基因 g 在不同sample（id）中的表达量与在sample（id）的表型值之间的线性关系（如上图）
其中，每个点对应不同的sample（id）；每个点对应的横坐标为基因 g 在不同sample中的表达量；每个点对应的纵坐标为不同sample对应的表型值
该表代表基因与表型的关系，那么回归系数 γ 的为正，那么代表基因表达量越高，则对性状的影响成正相关；反之为负，则代表基因表达量越高，则对性状的影响成负相关，后面的pval为回归系数的显著性
注: PrediXcan的使用
运行PrediXcan需要输入三个文件：转录组表达矩阵，基因型文件和样本信息文件：
- 基因型文件：该文件每一行表示一个SNP，包含的信息分别为：chromosome rsid position allele1 allele2 MAF，后面的每一列的内容是每一个样本在该SNP allele2的dosage，最好是每一条染色体分开制作文件。
- 样本信息文件：直接将PLINK的fam文件导入即可。
- 基因表达矩阵
可选，利用已有的模型预测一个基因表达矩阵
./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）
这种方法不依赖于call出来的snp和转录组数据，而是直接利用GWAS-summary数据来建立基因表达与表型值之间的关系
我们知道，GWAS-summary描述的是不同的snp（基因型）和表型值之间的关系
那么模型基于已经训练好的snp与基因 g 表达量之间的权重关系，推测该权重关系是基于PredictDB训练好的权重进行计算的
其中，wIg 代表第 I 个snp对gene g表达量产生影响的权重（该权重即snp I 是否突变对gene g表达量产生变化的回归系数）；βl 代表第 l 个snp 对表型值的影响的回归系数（也称为效应值）；se(βl) 代表所有回归系数（效应值）的标准误；σl 代表所有回归系数（效应值）的标准差；σg 代表基因 g 在各个样本中表达量的标准差；Zg 即为基因 g 对表型值的回归系数（效应值）
其中，wIg 是已经训练好的snp与基因 g 表达量之间的权重
因此 Zg 为基因 g 表达量和表型值之间线性模型的回归系数，越大说明对表型的正向影响越大；反之越小代表对表型的负向影响越大
