表达数量性状位点（eQTL）的概念及其相关分析

表达数量性状位点（expression quantitative trait locus, eQTL）是一类能够影响基因表达量的遗传位点（大部分都是单核苷酸多态性，SNP），具有一定的生物学意义。迄今为止最全的eQTL数据库是GTEx（https://www.gtexportal.org/home/），如今已更新到第八版了。

一般而言，eQTL主要分为两类：
（1）顺式eQTL（cis-eQTL）：它主要是指与所调控基因相距较近的eQTL，一般多位于所调控基因的上下游1Mb区域；
（2）反式eQTL（trans-eQTL）：与cis-eQTL恰恰相反，反式是指距离所调控基因位置比较远的eQTL，有时候距离甚至超过5Mb。

因此，对于eQTL分析而言，我们通常需要考虑两点，SNP和基因表达水平的关联度以及SNP与基因的距离。

由于大量eQTL数据库的开发，我们现在可以直接利用别人的结果寻找SNP调控的基因，最常用的就是GTEx数据库了，大家可以自行学习了解。接下来我将介绍如何利用自己的数据计算并确定相关eQTL。

利用原始数据做eQTL分析，我们至少需要三个文件，
第一个是样本信息文件，该文件包含样本的年龄，性别和人种等等；
第二个是基因表达量文件，它表示的是每个基因在每个样本中的表达含量；
第三个是基因型数据，也即每个样本的基因型数据。

接下来我们看看数据格式：

第一个是样本信息文件，除开第一列，其它列都代表不同的样本，每一行代表的是样本的表型信息。

第二个是基因表达量信息，同样地，除开第一列，其它列都代表不同的样本，每一行代表的是不同的基因（一般来说基因表达数据需要先进行标准化转换）。

第三个是基因型数据，同样地，除开第一列，其它列都代表不同的样本，每一行代表的是不同的基因型（SNP），一般基因型数据用0，1，2这三个数字编码，代表的是效应等位基因剂量。举个简单的例子，SNP1的等位基因分别是A和C，如果我们以A为效应等位基因，那么基因型AA的剂量便是2，AC为1，CC为0。

有了这些数据，我们就可以简单分析SNP和基因表达量的关系了

其数学模型如下：

gene1 ~ snp1 + sex + age + error_term

这里gene1（因变量）一般就是一个基因的表达量，snp1（自变量）就是一个SNP的基因型，两者拟合，矫正相关干扰项（如sex和age等），error_term是指回归模型的误差项。

如果想区分顺式还是反式eQTL，这时候就需要结合基因与SNP的位置信息了。

利用“MatrixEQTL”包进行eQTL实战分析

这里我们使用的是该包提供的内置数据集，代码如下：

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="") # 获取基因表达量文件位置
expression_file = data.table::fread(expression_file_name, header=T) # 读取基因表达量文件，可以在R中查看
covariates_file_name = paste(base.dir, "/data/Covariates.txt", sep="") # 读取协变量文件
covariates_file = data.table::fread(covariates_file_name, header=T) # 读取协变量文件，可在R中查看
output_file_name = tempfile() # 将输出文件设置为临时文件
pvOutputThreshold = 1e-2 # 定义gene-SNP associations的显著性P值
errorCovariance = numeric() # 定义误差项的协方差矩阵 (用的很少)
snps = SlicedData$new() # 创建SNP文件为S4对象（S4对象是该包的默认输入方式）
snps$fileDelimiter = "\t"      # 指定SNP文件的分隔符
snps$fileOmitCharacters = "NA" # 定义缺失值
snps$fileSkipRows = 1          # 跳过第一行（适用于第一行是列名的情况）
snps$fileSkipColumns = 1       # 跳过第一列（适用于第一列是SNP ID的情况
snps$fileSliceSize = 2000      # 每次读取2000条数据
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 )
# 文件的输入部分结束
me = Matrix_eQTL_engine( # 这是进行eQTL分析的主要函数
    snps = snps, # 指定SNP 文件
    gene = gene, # 指定基因表达量文件
    cvrt = cvrt, # 指定协变量文件
    output_file_name = output_file_name, # 指定输出文件
    pvOutputThreshold = pvOutputThreshold, # 指定显著性P值
    useModel = useModel, # 指定使用的计算模型
    errorCovariance = errorCovariance, # 指定误差项的协方差矩阵
    verbose = TRUE,
    pvalue.hist = TRUE,
    min.pv.by.genesnp = FALSE,
    noFDRsaveMemory = FALSE)
res <- me$all$eqtls # 把eQTL的显著结果存储到变量res里
res # 查看结果

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