前言

我們知道GWAS描述的是snp與表型之間的關(guān)系，即利用線性模型尋找與某種性狀顯著相關(guān)聯(lián)的位點(diǎn)。但是發(fā)揮生物學(xué)功能的往往是蛋白質(zhì)，蛋白質(zhì)是由轉(zhuǎn)錄本翻譯而來(lái)，那么建立基因表達(dá)量與表型的關(guān)聯(lián)將會(huì)使得分析更進(jìn)一步

這張圖講述了GWAS與轉(zhuǎn)錄組數(shù)據(jù)整合分析的流程

基于snp和轉(zhuǎn)錄組數(shù)據(jù)（PrediXcan）

可參考文章A gene-based association method for mapping traits using reference transcriptome data
首先我們先從原理上介紹下這款軟件的工作流程

第一步
我們對(duì)測(cè)序數(shù)據(jù) call 出來(lái)snp以后會(huì)得到一張表格：

顯然 C 突變?yōu)?T （C/T為雜合突變，T/T為純和突變）會(huì)促使表型值上升，與表型值成正相關(guān)

那么，判斷突變與表型值呈正相關(guān)還是負(fù)相關(guān)，我們可以設(shè)置一個(gè)對(duì)照（比方說(shuō)設(shè)置純和未突變的為對(duì)照）看看突變是否會(huì)引起表型值的上升

并且測(cè)序數(shù)據(jù) call 出來(lái)的snp信息可以和表型值相關(guān)聯(lián)

而我們RNA-seq的表格為：

第二步（可選，若用戶沒(méi)用基因表達(dá)矩陣執(zhí)行此步驟）
利用snp數(shù)據(jù)和轉(zhuǎn)錄組各個(gè)基因表達(dá)量的數(shù)據(jù)，建立它們之間的線性的權(quán)重關(guān)系

其中，w_k,g 代表第 k 個(gè)snp與第 g 個(gè)基因表達(dá)量之間的權(quán)重（權(quán)重 w_k,g 是作者利用機(jī)器學(xué)習(xí)的思想，利用GTEx Project, GEUVADIS 和 DGN數(shù)據(jù)庫(kù)中基因型數(shù)據(jù)和基因表達(dá)數(shù)據(jù)做訓(xùn)練集，即利用已知的基因表達(dá)量 T_g 和 snp 的基因型數(shù)據(jù) X_k 通過(guò) LASSO 和 elastic net 來(lái)計(jì)算權(quán)重 w_k,g）；T_g 代表第 g 個(gè)基因的表達(dá)量；X_k 代表第 k 個(gè)snp的基因型（因子型變量轉(zhuǎn)換為[數(shù)值型變量] 0，1，2）

因此如果你沒(méi)有現(xiàn)成的轉(zhuǎn)錄組數(shù)據(jù)，你可以利用作者已開(kāi)發(fā)好的模型（利用機(jī)器學(xué)習(xí)的方法開(kāi)發(fā)了部分模型）。作者已經(jīng)利用已發(fā)的數(shù)據(jù)做了模型訓(xùn)練了，如果你有對(duì)應(yīng)組織的snp數(shù)據(jù)，可以到PredictDB下載對(duì)應(yīng)的模型，用于預(yù)測(cè)基因的表達(dá)量信息。

第三步
關(guān)于權(quán)重 w_k,g 的計(jì)算可以利用LASSO和 elastic net 來(lái)計(jì)算，由此可知，對(duì)于某一個(gè) Tissue 來(lái)說(shuō)，第 g 個(gè)基因的表達(dá)量可以用snp來(lái)線性表示，那么結(jié)合trait的值（我們?cè)谶@里稱為trait的表達(dá)量），再次擬合一個(gè)線性模型：

其中 Tg 第 g 個(gè)基因的表達(dá)量；γ 代表回歸系數(shù)；Y 代表表型值（表型表達(dá)量））；Y_n = γ₁×T_1,n+ γ₂×T_2,n + ... + γ_m×T_m,n，n代表第n個(gè)id，m代表m個(gè)基因
這樣就可以將基因表達(dá)量與表型值聯(lián)系起來(lái)了，其本質(zhì)就是基因表達(dá)量和表型直接的關(guān)系

對(duì)于該模型，我們可以這樣理解，對(duì)于每一個(gè)基因 g 來(lái)說(shuō)，在各個(gè)sample（id）中的表達(dá)量不同，而每一個(gè)sample（id）的表型值也不同，因此可以建立基因 g 在不同sample（id）中的表達(dá)量與在sample（id）的表型值之間的線性關(guān)系（如上圖）
其中，每個(gè)點(diǎn)對(duì)應(yīng)不同的sample（id）；每個(gè)點(diǎn)對(duì)應(yīng)的橫坐標(biāo)為基因 g 在不同sample中的表達(dá)量；每個(gè)點(diǎn)對(duì)應(yīng)的縱坐標(biāo)為不同sample對(duì)應(yīng)的表型值

該表代表基因與表型的關(guān)系，那么回歸系數(shù) γ 的為正，那么代表基因表達(dá)量越高，則對(duì)性狀的影響成正相關(guān)；反之為負(fù)，則代表基因表達(dá)量越高，則對(duì)性狀的影響成負(fù)相關(guān)，后面的pval為回歸系數(shù)的顯著性

注: PrediXcan的使用

運(yùn)行PrediXcan需要輸入三個(gè)文件：轉(zhuǎn)錄組表達(dá)矩陣，基因型文件和樣本信息文件：

1. 基因型文件：該文件每一行表示一個(gè)SNP，包含的信息分別為：chromosome rsid position allele1 allele2 MAF，后面的每一列的內(nèi)容是每一個(gè)樣本在該SNP allele2的dosage，最好是每一條染色體分開(kāi)制作文件。
2. 樣本信息文件：直接將PLINK的fam文件導(dǎo)入即可。
3. 基因表達(dá)矩陣

可選，利用已有的模型預(yù)測(cè)一個(gè)基因表達(dá)矩陣

./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

建立基因表達(dá)與表型的關(guān)系

./PrediXcan.py 
--assoc 
--pheno My_pheno.txt 
--mpheno 1 
--pred_exp results/TW_Brain_Frontal_predicted_expression.txt #預(yù)測(cè)的基因表達(dá)矩陣，作者也可自行提供
--logistic 
--output_prefix results/DGN-HapMap

基于GWAS-summary數(shù)據(jù)（Summary-PrediXcan）

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

這種方法不依賴于call出來(lái)的snp和轉(zhuǎn)錄組數(shù)據(jù)，而是直接利用GWAS-summary數(shù)據(jù)來(lái)建立基因表達(dá)與表型值之間的關(guān)系
我們知道，GWAS-summary描述的是不同的snp（基因型）和表型值之間的關(guān)系

其中，wIg 代表第 I 個(gè)snp對(duì)gene g表達(dá)量產(chǎn)生影響的權(quán)重（該權(quán)重即snp I 是否突變對(duì)gene g表達(dá)量產(chǎn)生變化的回歸系數(shù)）；βl 代表第 l 個(gè)snp 對(duì)表型值的影響的回歸系數(shù)（也稱為效應(yīng)值）；se(βl) 代表所有回歸系數(shù)（效應(yīng)值）的標(biāo)準(zhǔn)誤；σl 代表所有回歸系數(shù)（效應(yīng)值）的標(biāo)準(zhǔn)差；σg 代表基因 g 在各個(gè)樣本中表達(dá)量的標(biāo)準(zhǔn)差；Zg 即為基因 g 對(duì)表型值的回歸系數(shù)（效應(yīng)值）

其中，wIg 是已經(jīng)訓(xùn)練好的snp與基因 g 表達(dá)量之間的權(quán)重

因此 Zg 為基因 g 表達(dá)量和表型值之間線性模型的回歸系數(shù)，越大說(shuō)明對(duì)表型的正向影響越大；反之越小代表對(duì)表型的負(fù)向影響越大