Matrix eQTL分析

目前针对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.