目前针对QTL分析有很多软件，比如Plink、 Merlin、R/qtl、FastMap等等,今天给大家分享一下如何使用Matrix eQTL进行分析。首先我们先了解一下什么是eQTL以及eQTL的分类。
表达数量性状基因座（expression Quantitative Trait Loci, eQTL）：指染色体上一些能特定调控mRNA和蛋白质表达水平的区域，其mRNA/蛋白质的表达水平量与数量性状成比例关系。其中eQTL又分为cis-eQTL和trans-eQTL。
顺式作用eQTL（cis-eQTL）：指某个基因的eQTL定位到该基因所在的基因组区域，表明可能是该基因本身的差别引起的mRNA水平变化。
**反式作用eQTL(trans-eQTL): ** 指某个基因的eQTL定位到其他基因组区域，表明其他基因的差别控制该基因mRNA水平的差异。
Matrix eQTL是一款针对大型矩阵可以超快运行进行eQTL分析的软件(http://www.bios.unc.edu/research/genomic_software/Matrix_eQTL/)。
运行这款软件，需要提前准备5个文件文件，基因型文件SNP.txt, 表达文件GE.txt, 协变量文件Covariates.txt, 基因位置文件geneloc.txt和SNP位置文件snploc.txt。我们先看看这几个文件是什么格式。
file 1：SNP.txt
SNP.txt.jpg
行代表SNP，列代表个体
file 2: snpsloc.txt
snploc.txt.jpg
包含3列，第一列SNP id, 第二列染色体， 第三列SNP位置
file 3: GE.txt
GE.txt.jpg
行代表基因表达量，列代表个体
file 4: geneloc.txt
geneloc.txt.jpg
包含3列，第一列基因id, 第二列染色体， 第三列基因起始位置，第三列基因终止位置。
file 5: Covariates.txt
Covariates.txt.jpg
协变量文件包含个体的性别和年龄。
接下来我们一起学习一下如何使用Matrix eQTL进行分析。
BiocManager::install("MatrixEQTL") #安装R包
测试所有基因-SNP对并绘制所有p值的直方图
base.dir = find.package("MatrixEQTL") #获取R包MatrixEQTL的路径
useModel = modelLINEAR #三种模型可选(modelANOVA or modelLINEAR or modelLINEAR_CROSS),这里我们选择modelLINEAR
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文件
#加载基因表达文件
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)
#结果
cat('Analysis done in: ', me$time.in.sec, ' seconds', '\n'); #查看分析所用时间
cat('Detected eQTLs:', '\n');
show(me$all$eqtls) #查看超过阈值的eQTL
#画所有p-value直方图
plot(me)
Histogram for all p value.jpeg
Test local and distand gene-SNP pairs separately and plot Q-Q plots of local and distant p-values
Matrix eQTL可以区分顺式（cis，local）和反式（trans，distant）eQTL，主要用Matrix_eQTL_main函数来分析。
library(MatrixEQTL) #加载R包
base.dir = find.package('MatrixEQTL'); #找到示例数据所在路径
###设置
#使用modelLINEAR模型
useModel = modelLINEAR; # modelANOVA, modelLINEAR, or modelLINEAR_CROSS
#基因型文件名
SNP_file_name = paste(base.dir, "/data/SNP.txt", sep="");
snps_location_file_name = paste(base.dir, "/data/snpsloc.txt", sep="");
#基因表达文件名
expression_file_name = paste(base.dir, "/data/GE.txt", sep="");
gene_location_file_name = paste(base.dir, "/data/geneloc.txt", sep="");
#协变量文件名
covariates_file_name = paste(base.dir, "/data/Covariates.txt", sep=""); #如果没有协变量，设置为character()
#输出文件名，分为顺式和反式
output_file_name_cis = tempfile(); #顺式
output_file_name_tra = tempfile(); #反式
#只有在高于阈值重要的gene-SNP关联对才会被保存
pvOutputThreshold_cis = 2e-2;
pvOutputThreshold_tra = 1e-2;
#误差协方差矩阵
#设置为numeric()
errorCovariance = numeric();
#errorCovariance = read.table("Sample_Data/errorCovariance.txt");
#Distance for local gene-SNP pairs
cisDist = 1e6;
#加载基因型数据
snps = SlicedData$new();
snps$fileDelimiter = "\t"; # 指定SNP文件的分隔符Tab键
snps$fileOmitCharacters = "NA"; #定义缺失值;
snps$fileSkipRows = 1; #跳过第一行（适用于第一行是列名的情况）
snps$fileSkipColumns = 1; #跳过第一列（适用于第一列是SNP ID的情况）snps$fileSliceSize = 2000; #每次读取2000条数据snps$LoadFile(SNP_file_name); #载入SNP文件
#加载基因表达量数据
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;
if(length(covariates_file_name)>0) {
cvrt$LoadFile(covariates_file_name);
}
#读取基因型和基因位置文件
snpspos = read.table(snps_location_file_name, header = TRUE, stringsAsFactors = FALSE);
genepos = read.table(gene_location_file_name, header = TRUE, stringsAsFactors = FALSE);
#运行分析
me = Matrix_eQTL_main(
snps = snps,
gene = gene,
cvrt = cvrt,
output_file_name = output_file_name_tra,
pvOutputThreshold = pvOutputThreshold_tra,
useModel = useModel,
errorCovariance = errorCovariance,
verbose = TRUE,
output_file_name.cis = output_file_name_cis,
pvOutputThreshold.cis = pvOutputThreshold_cis,
snpspos = snpspos,
genepos = genepos,
cisDist = cisDist,
pvalue.hist = "qqplot",
min.pv.by.genesnp = FALSE,
noFDRsaveMemory = FALSE);
unlink(output_file_name_tra);
unlink(output_file_name_cis);
#结果:
cat('Analysis done in: ', me$time.in.sec, ' seconds', '\n');
cat('Detected local eQTLs:', '\n');
show(me$cis$eqtls) #检测到的顺式eQTLs
snps gene statistic pvalue FDR beta
1 Snp_05 Gene_03 38.812160 5.515519e-14 5.515519e-12 0.4101317
2 Snp_04 Gene_10 3.201666 7.608981e-03 3.804491e-01 0.2321123
cat('Detected distant eQTLs:', '\n');
show(me$trans$eqtls) #检测到的反式eQTLs
snps gene statistic pvalue FDR beta
1 Snp_13 Gene_09 -3.914403 0.002055817 0.1027908 -0.2978847
2 Snp_11 Gene_06 -3.221962 0.007327756 0.1619451 -0.2332470
3 Snp_14 Gene_01 3.070005 0.009716705 0.1619451 0.2147077
#Plot the Q-Q plot of local and distant p-values
plot(me)
QQ-plot.jpeg
参考文献：
Shabalin A A. Matrix eQTL: ultra fast eQTL analysis via large matrix operations[J]. Bioinformatics, 2012, 28(10): 1353-1358.
