在上一期内容中,米老鼠和大家介绍了eQTL的相关概念和分析原理,今天我就带大家用“MatrixEQTL”进行一下实战演练。
这里我们使用的是该包提供的内置数据集,代码如下:
install.packages("MatrixEQTL")#安装R包library("MatrixEQTL")#加载R包base.dir=find.package("MatrixEQTL")#获取该包的路径信息useModel=modelLINEAR#三种模型可选(modelANOVAormodelLINEARormodelLINEAR_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-SNPassociations的显著性P值errorCovariance=numeric()#定义误差项的协方差矩阵(用的很少)snps=SlicedData$new()#创建SNP文件为S4对象(S4对象是该包的默认输入方式)snps$fileDelimiter="\t"#指定SNP文件的分隔符snps$fileOmitCharacters="NA"#定义缺失值snps$fileSkipRows=1#跳过第一行(适用于第一行是列名的情况)snps$fileSkipColumns=1#跳过第一列(适用于第一列是SNPID的情况snps$fileSliceSize=#每次读取条数据snps$LoadFile(SNP_file_name)#载入SNP文件#下面的代码同snps的那部分一致gene=SlicedData$new()gene$fileDelimiter="\t"gene$fileOmitCharacters="NA"gene$fileSkipRows=1gene$fileSkipColumns=1gene$fileSliceSize=gene$LoadFile(expression_file_name)cvrt=SlicedData$new()cvrt$fileDelimiter="\t"cvrt$fileOmitCharacters="NA"cvrt$fileSkipRows=1cvrt$fileSkipColumns=1cvrt$fileSliceSize=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#查看结果
关于“MatrxiEQTL”包的用法就简单介绍到这里,感兴趣的朋友可以深入学习一下,掌握cis-eQTL和trans-eQTL的分析方法。
生信与临床