表達(dá)數(shù)量性狀位點(diǎn)（expression quantitative trait locus, eQTL）是一類能夠影響基因表達(dá)量的遺傳位點(diǎn)（大部分都是單核苷酸多態(tài)性，SNP），具有一定的生物學(xué)意義。迄今為止最全的eQTL數(shù)據(jù)庫是GTEx（https://www.gtexportal.org/home/），如今已更新到第八版了。

一般而言，eQTL主要分為兩類：
（1）順式eQTL（cis-eQTL）：它主要是指與所調(diào)控基因相距較近的eQTL，一般多位于所調(diào)控基因的上下游1Mb區(qū)域；
（2）反式eQTL（trans-eQTL）：與cis-eQTL恰恰相反，反式是指距離所調(diào)控基因位置比較遠(yuǎn)的eQTL，有時候距離甚至超過5Mb。

因此，對于eQTL分析而言，我們通常需要考慮兩點(diǎn)，SNP和基因表達(dá)水平的關(guān)聯(lián)度以及SNP與基因的距離。

由于大量eQTL數(shù)據(jù)庫的開發(fā)，我們現(xiàn)在可以直接利用別人的結(jié)果尋找SNP調(diào)控的基因，最常用的就是GTEx數(shù)據(jù)庫了，大家可以自行學(xué)習(xí)了解。接下來我將介紹如何利用自己的數(shù)據(jù)計算并確定相關(guān)eQTL。

利用原始數(shù)據(jù)做eQTL分析，我們至少需要三個文件，
第一個是樣本信息文件，該文件包含樣本的年齡，性別和人種等等；
第二個是基因表達(dá)量文件，它表示的是每個基因在每個樣本中的表達(dá)含量；
第三個是基因型數(shù)據(jù)，也即每個樣本的基因型數(shù)據(jù)。

接下來我們看看數(shù)據(jù)格式：

第一個是樣本信息文件，除開第一列，其它列都代表不同的樣本，每一行代表的是樣本的表型信息。

第二個是基因表達(dá)量信息，同樣地，除開第一列，其它列都代表不同的樣本，每一行代表的是不同的基因（一般來說基因表達(dá)數(shù)據(jù)需要先進(jìn)行標(biāo)準(zhǔn)化轉(zhuǎn)換）。

第三個是基因型數(shù)據(jù)，同樣地，除開第一列，其它列都代表不同的樣本，每一行代表的是不同的基因型（SNP），一般基因型數(shù)據(jù)用0，1，2這三個數(shù)字編碼，代表的是效應(yīng)等位基因劑量。舉個簡單的例子，SNP1的等位基因分別是A和C，如果我們以A為效應(yīng)等位基因，那么基因型AA的劑量便是2，AC為1，CC為0。

有了這些數(shù)據(jù)，我們就可以簡單分析SNP和基因表達(dá)量的關(guān)系了

其數(shù)學(xué)模型如下：

gene1 ~ snp1 + sex + age + error_term

這里gene1（因變量）一般就是一個基因的表達(dá)量，snp1（自變量）就是一個SNP的基因型，兩者擬合，矯正相關(guān)干擾項(xiàng)（如sex和age等），error_term是指回歸模型的誤差項(xiàng)。

如果想?yún)^(qū)分順式還是反式eQTL，這時候就需要結(jié)合基因與SNP的位置信息了。

利用“MatrixEQTL”包進(jìn)行eQTL實(shí)戰(zhàn)分析

這里我們使用的是該包提供的內(nèi)置數(shù)據(jù)集，代碼如下：

install.packages("MatrixEQTL") # 安裝R包
library("MatrixEQTL") # 加載R包
base.dir = find.package("MatrixEQTL") # 獲取該包的路徑信息
useModel = modelLINEAR # 三種模型可選(modelANOVA or modelLINEAR or modelLINEAR_CROSS)
SNP_file_name = paste(base.dir, "/data/SNP.txt", sep="") # 獲取SNP文件位置
SNP_file = data.table::fread(SNP_file_name, header=T) # 讀取SNP文件，可以在R中查看
expression_file_name = paste(base.dir, "/data/GE.txt", sep="") # 獲取基因表達(dá)量文件位置
expression_file = data.table::fread(expression_file_name, header=T) # 讀取基因表達(dá)量文件，可以在R中查看
covariates_file_name = paste(base.dir, "/data/Covariates.txt", sep="") # 讀取協(xié)變量文件
covariates_file = data.table::fread(covariates_file_name, header=T) # 讀取協(xié)變量文件，可在R中查看
output_file_name = tempfile() # 將輸出文件設(shè)置為臨時文件
pvOutputThreshold = 1e-2 # 定義gene-SNP associations的顯著性P值
errorCovariance = numeric() # 定義誤差項(xiàng)的協(xié)方差矩陣 (用的很少)
snps = SlicedData$new() # 創(chuàng)建SNP文件為S4對象（S4對象是該包的默認(rèn)輸入方式）
snps$fileDelimiter = "\t"      # 指定SNP文件的分隔符
snps$fileOmitCharacters = "NA" # 定義缺失值
snps$fileSkipRows = 1          # 跳過第一行（適用于第一行是列名的情況）
snps$fileSkipColumns = 1       # 跳過第一列（適用于第一列是SNP ID的情況
snps$fileSliceSize = 2000      # 每次讀取2000條數(shù)據(jù)
snps$LoadFile( SNP_file_name ) # 載入SNP文件
# 下面的代碼同snps的那部分一致
gene = SlicedData$new()
gene$fileDelimiter = "\t"
gene$fileOmitCharacters = "NA"
gene$fileSkipRows = 1
gene$fileSkipColumns = 1
gene$fileSliceSize = 2000
gene$LoadFile( expression_file_name )
cvrt = SlicedData$new()
cvrt$fileDelimiter = "\t"
cvrt$fileOmitCharacters = "NA"
cvrt$fileSkipRows = 1
cvrt$fileSkipColumns = 1
cvrt$fileSliceSize = 2000
cvrt$LoadFile( covariates_file_name )
# 文件的輸入部分結(jié)束
me = Matrix_eQTL_engine( # 這是進(jìn)行eQTL分析的主要函數(shù)
    snps = snps, # 指定SNP 文件
    gene = gene, # 指定基因表達(dá)量文件
    cvrt = cvrt, # 指定協(xié)變量文件
    output_file_name = output_file_name, # 指定輸出文件
    pvOutputThreshold = pvOutputThreshold, # 指定顯著性P值
    useModel = useModel, # 指定使用的計算模型
    errorCovariance = errorCovariance, # 指定誤差項(xiàng)的協(xié)方差矩陣
    verbose = TRUE,
    pvalue.hist = TRUE,
    min.pv.by.genesnp = FALSE,
    noFDRsaveMemory = FALSE)
res <- me$all$eqtls # 把eQTL的顯著結(jié)果存儲到變量res里
res # 查看結(jié)果

參考：
https://zhuanlan.zhihu.com/p/378403055
https://zhuanlan.zhihu.com/p/378403493
https://blog.csdn.net/weixin_43569478/article/details/108079878