什么是eQTL分析？
eQTL全称“Expression Quantitative Trait Locianalysis”，即“表达数量位置的基因座”，它指的是染色体上一些能特定调控mRNA和蛋白表达水平的区域，其mRNA/蛋白质的表达水平量与数量性状成比例关系。eQTL analysis是将基因表达水平的变化和基因型连接起来，研究遗传突变与基因表达的相关性。
eQTL可以分为顺式和反式两种，总的来说，顺式eQTL指的是那些和所调控基因距离非常近的SNP位点，一般多位于所调控基因的上下游1Mb区域，通常作用于局部基因的遗传变异，而反式eQTL则是与顺式eQTL相对的一种位点，通常距离所调控的基因距离较远（见下图）：
eQTL分析的本质是以全部的DNA变异位点为自变量，轮流以每种mRNA表达量为因变量，用大量的个体数据做样本进行线性回归，得到每一个SNP位点和每一个mRNA表达量间的关系。
使用MatrixEQTL进行eQTL分析
这里需要注意的是表达量文件的横坐标是基因，纵坐标是对应的样本，需要提前对表达量数据进行标准化处理
协变量文件：
这里每一行为一个协变量，每一列为一个样本，大家可以自行决定是否是用协变量对结果进行校正
基因的位置文件：
分别为“基因id”、“染色体”、“起始位置”
SNP的位置坐标文件：
分别为“snpid”、“染色体”、“位置”
下面是使用MatrixEQTL进行eQTL分析的脚本，把上面准备好的文件直接带入进相应的变量就可以了：
################PCA计算主成分，作为分群的协变量#加载R包library("MatrixEQTL")##读入数据#data.table::fread读取大型文件更高效SNP <- data.table::fread("LDfiltered.012.vcf",sep="\t")#Express <- data.table::fread("Express_data_Normalize.csv",sep="\t")#Cov <- data.table::fread("covariates.txt",sep="\t")###处理SNP数据#取需要的几列dim(SNP)new_SNP <- SNP[,-(1:2)]my_SNP <- new_SNP[,-(2:7)]#SNP_scaled <- scale(SNP)write.table(my_SNP,"SNP.txt",sep="\t",quote = FALSE,col.names = TRUE,row.names=FALSE)##设置模型、P值等参数useModel = modelLINEAR # 三种模型可选(modelANOVA or modelLINEAR or modelLINEAR_CROSS)cisDist=1e7#判断基因SNP属于cis的距离阈值output_file_name_tra="trans-eqtl.txt"pvOutputThreshold_tra=1e-2#trans-eqtl-pvalue阈值output_file_name_cis="cis-eqtl.txt"pvOutputThreshold_cis=5e-2#cis-eqtl-pvalue阈值#output_file_name = tempfile() # 将输出文件设置为临时文件#pvOutputThreshold = 1e-2 # 定义gene-SNP associations的显著性P值errorCovariance = numeric() # 定义误差项的协方差矩阵 (用的很少)##设置SNP文件参数SNP_file="SNP.txt"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 ) # 载入SNP文件##设置表达量文件参数expression_file="Express_data_Normalize.csv"expr = SlicedData$new()expr$fileDelimiter = "\t"expr$fileOmitCharacters = "NA"#gene$fileSkipRows = 1expr$fileSkipColumns = 1expr$fileSliceSize = 2000expr$LoadFile( expression_file )##设置协变量文件参数covariates_file="all_covariates.2.csv"cvrt = SlicedData$new()cvrt$fileDelimiter = "\t"cvrt$fileOmitCharacters = "NA"#cvrt$fileSkipRows = 1cvrt$fileSkipColumns = 1cvrt$fileSliceSize = 2000cvrt$LoadFile( covariates_file )#读取基因染色体位置文件gene_location_file_name="genesloc.txt"genepos=read.table(gene_location_file_name,header=TRUE,stringsAsFactors = FALSE)#读取SNP位置文件snps_location_file_name="snpsloc.txt"snpspos=read.table(snps_location_file_name,header=TRUE,stringsAsFactors = FALSE)##进行eQTL计算my_eQTL = Matrix_eQTL_main( # 这是进行eQTL分析的主要函数snps = snps, # 指定SNP 文件gene = expr, # 指定基因表达量文件cvrt = cvrt, # 指定协变量文件cisDist=cisDist,snpspos=snpspos,genepos=genepos,output_file_name = output_file_name_tra, # 指定输出文件pvOutputThreshold = pvOutputThreshold_tra, # 指定显著性P值output_file_name.cis = output_file_name_cis, # 指定输出文件pvOutputThreshold.cis = pvOutputThreshold_cis, # 指定显著性P值useModel = useModel, # 指定使用的计算模型errorCovariance = errorCovariance, # 指定误差项的协方差矩阵verbose = TRUE,pvalue.hist = "qqplot",min.pv.by.genesnp = FALSE,noFDRsaveMemory = FALSE)#out <- my_eQTL$all$eqtls # 把eQTL的显著结果存储到变量res里#out # 查看结果pdf("p_value_qqplots.pdf")plot(my_eQTL)dev.off()
参考资料：
eQTL分析相关原理：
https://xsliulab.github.io/Workshop/2021/week32/eqtl%E8%AE%A1%E7%AE%97%E6%96%B9%E6%B3%95.html
