R语言实现eQTL分析
今天给大家介绍一个GWAS分析过程中的一个重要的环节eQTL（表达数量性状位点）分析。eQTL指的是染色体上一些能特定调控mRNA和蛋白质表达水平的区域，其mRNA/蛋白质的表达水平量与数量性状成比例关系，通俗点讲就是把基因表达作为一种性状,研究遗传突变与基因表达的相关性。
接下来我们看下在R语言中如何实现这个过程。我们需要准备以下几个R包：
install.packages('SNPassoc')#获取SNPs原始数据
install.packages('synbreed')#将原始数据转化为0，1，2编码的矩阵。最重要的是这个包需要在高版本R语言安装最好是3.5以上版本。
install.packages('MatrixEQTL')#eQTL分析
然后我们看下流程分解：
首先我们需要载入SNPassoc中的原始数据：
library(SNPassoc)
data(SNPs)
SNPs[1:8,1:8]
SNPs数据的提取：
myDat<- SNPs[,-(2:5)]
row.names(myDat) <- myDat$id;
myDat <- myDat[,-1]
myDat[1:5,1:5]
myDat <- as.matrix(myDat)#将数据格式进行转化并且填补空值
数据准备完毕，接下来是数据的预处理，利用synbreed通过等位基因频率对基因型进行转化，将频率高的纯合体为0，杂合为1，频率低的纯合体为2。
library(synbreed)
cp <- create.gpData(geno = myDat)
cp.dat <- codeGeno(gpData =cp,label.heter = "alleleCoding", maf = 0.01, nmiss = 0.1,
impute = TRUE, impute.type ="random", verbose = TRUE)
gee <- cp.dat$geno#获取结果数据
最后就是eQTL分析，利用MatrixEQTL进行分析。我们需要先整理好预先用的数据。
我们在此以其自身所带的示例数据为基础展开说明，主要有三部分：
1. 数据
前期数据文件以及相关参数的设置：
#SNP数据
SNP_file_name = paste0(base.dir,"/data/SNP.txt");
snps_location_file_name = paste0(base.dir,"/data/snpsloc.txt");
# 基因表达数据
expression_file_name = paste0(base.dir,"/data/GE.txt");
gene_location_file_name = paste0(base.dir,"/data/geneloc.txt");
# 表型数据
# Set to character() for no covariates
covariates_file_name = paste0(base.dir,"/data/Covariates.txt");
#输出文件
output_file_name_cis = tempfile();
output_file_name_tra = tempfile();
#关联水平阈值
pvOutputThreshold_cis = 2e-2;
pvOutputThreshold_tra = 1e-2;
#误差协方差数据存储
# Set to numeric() for identity.
errorCovariance = numeric();
# errorCovariance =read.table("Sample_Data/errorCovariance.txt");
# 基因-SNP距离阈值
cisDist = 1e6;
#模型的选择modelANOVA,modelLINEAR, or modelLINEAR_CROSS。modelANOVA顾名思义是基于F-检验进行的分析；modelLINEAR是基于T-检验进行的分析；modelLINEAR_CROSS主要是在线性回归的基础上加了交叉项。最终还是基于T-检验。
useModel = modelLINEAR;
## SNP数据结构构建
snps = SlicedData$new();
snps$fileDelimiter = "\t"; # the TAB character
snps$fileOmitCharacters = "NA"; #denote missing values;
snps$fileSkipRows = 1; # one row of column labels
snps$fileSkipColumns = 1; # one column of row labels
snps$fileSliceSize = 2000; # read file in slices of 2,000 rows
snps$LoadFile(SNP_file_name);
## 基因表达数据构建
gene = SlicedData$new();
gene$fileDelimiter = "\t"; # the TAB character
gene$fileOmitCharacters = "NA"; #denote missing values;
gene$fileSkipRows = 1; # one row of column labels
gene$fileSkipColumns = 1; # one column of row labels
gene$fileSliceSize = 2000; # read file in slices of 2,000 rows
gene$LoadFile(expression_file_name);
## 表型数据
cvrt = SlicedData$new();
cvrt$fileDelimiter = "\t"; # the TAB character
cvrt$fileOmitCharacters = "NA"; #denote missing values;
cvrt$fileSkipRows = 1; # one row of column labels
cvrt$fileSkipColumns = 1; # one column of row labels
if(length(covariates_file_name)>0) {
cvrt$LoadFile(covariates_file_name);
}
## 基因和SNP位置信息
snpspos =read.table(snps_location_file_name, header = TRUE, stringsAsFactors = FALSE);
genepos =read.table(gene_location_file_name, header = TRUE, stringsAsFactors = FALSE);
2. 数据分析
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 = TRUE,
min.pv.by.genesnp = FALSE,
noFDRsaveMemory= FALSE);
接下来数对结果的描述：
结果中我们可以得到顺式作用（cis-）与反式作用(trans-)的位点信息。所谓顺式作用即位点位于基因内部；反式作用也就是不再基因内部的SNPs。
顺式作用结果：
show(me$cis$eqtls)
反式作用结果：
show(me$trans$eqtls)
其中的beta是指标准回归系数，是用来比较各个系数之间的绝对作用或者贡献的大小；stattistic展示相关系数r的大小。
当然我们也可以将结果可视化，比如查看所有的SNP 的一个P值的分布情况：
plot(me)
