《微信小程序进阶实战之分答应用开发（中级项目）》（完整版）

找到Matrix eQTL这个包，看下文章Matrix eQTL: ultra fast eQTL analysis via large matrix operations（https://doi.org/10.1093/bioinformatics/bts163
）
eQTL（表达数量性状位点）计算transcript-SNP 的关系，即分析SNP与基因的表达是否相关。由于计算数量巨大，很多人都用较小的数据来做。因此该作者开发了Matrix eQTL，用于处理大数据，支持additive linear and ANOVA models with covariates，并且可以将cis- and trans-eQTLs分开计算。
Matrix eQTL相较于其他软件如FastMap — 18.4 min, Merlin — 12.3 min, Plink — 9.0 min, Matrix eQTL — 5.7 min and snpMatrix — 3.3 min要快，它设置一个阈值，只有超过这个阈值的p值才会被计算。
采用的是线型回归模型，g为基因表达情况，s为SNP分型结果。

说明文档http://www.bios.unc.edu/research/genomic_software/Matrix_eQTL/runit.html
示例数据：http://www.bios.unc.edu/research/genomic_software/Matrix_eQTL/R.html
分析过程很简单，首先设置好要分析文件的路径和名称：

install.packages("MatrixEQTL")

# 设置数据目录，示例数据放在包的安装目录下了。
base.dir = find.package("MatrixEQTL")

#设置分析的模型
useModel = modelLINEAR; # modelANOVA or modelLINEAR or modelLINEAR_CROSS
#设置SNP文件的名称
SNP_file_name = paste(base.dir, "/data/SNP.txt", sep="");
# 设置表达数据文件的名称
expression_file_name = paste(base.dir, "/data/GE.txt", sep="");
# 设置协变量文件的名称
# 无协变量设置为character()
covariates_file_name = paste(base.dir, "/data/Covariates.txt", sep="");

output_file_name = tempfile();
提供了三种分析模型供选择
(1) modelLINEAR
Model: useModel = modelLINEAR
Equation: expression = α + ∑k βk⋅covariatek + γ⋅genotype_additive
Testing for significance of: γ
Test statistic: t-statistic
(2) modelANOVA

Model: useModel = modelANOVA
Equation: expression = α + ∑k βk⋅covariatek + γ1⋅genotype_additive + γ2⋅genotype_dominant
Testing for significance of: (γ1,γ2) pair
Test statistic: F-statistic
(3) modelLINEAR_CROSS

Model: useModel = modelLINEAR_CROSS
Equation:
expression = α + ∑k βk⋅covariatek + γ⋅genotype_additive + δ⋅genotype_additive⋅covariateK
Testing for significance of: δ
Test statistic: t-statistic
注意这里要设置一个p值的阈值，一般越大的数据量阈值设的越小，之前说过它会按这个阈值来计算结果，如果设的过大，分析耗时并且输出很多结果。输出的结果都储存在output_file_name里
pvOutputThreshold = 1e-2
# 设置协变量矩阵为 numeric()，很少用，默认
errorCovariance = numeric()
# 这里建立了一个SlicedData的新对象，用于存放martix的数据，并设置存放数据的格式
snps = SlicedData$new();
# 设置数据分隔符为tab
snps$fileDelimiter = "\t"; # the TAB character
# 设置缺失值为NA
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 pieces of 2,000 rows
snps$LoadFile( SNP_file_name );

## Load gene expression data

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);

## Load covariates

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
看下文件格式，snp文件用0,1,2表示，基因文件是表达量，cvrt是covariates：

image.png

image.png
设置好文件后可以用 Matrix_eQTL_engine主函数进行eQTL分析了，参数snps设置SNP文件，gene设置表达量文件，cvrt设置协变量。然后将每行的SNP和gene放到一块进行线性回归的分析。

me = Matrix_eQTL_engine(
snps = snps,
gene = gene,
cvrt = cvrt,
output_file_name = output_file_name,
pvOutputThreshold = pvOutputThreshold,
useModel = useModel,
errorCovariance = errorCovariance,
verbose = TRUE,
pvalue.hist = TRUE,
min.pv.by.genesnp = FALSE,
noFDRsaveMemory = FALSE);
运行完后得到的me对象是一个list：

image.png
输出文件的每行eqtl为：SNP name, a transcript name, estimate of the effect size, t- or F-statistic, p-value, and FDR。

Matrix eQTL可以区分顺式（cis，local）和反式（trans，distant）eQTL，主要用Matrix_eQTL_main函数来分析。其包括以下几个参数：

* `pvOutputThreshold.cis` – p-value threshold for cis-eQTLs.
* `output_file_name.cis` – detected cis-eQTLs are saved in this file.
* `cisDist` – maximum distance at which gene-SNP pair is considered local.
* `snpspos` – data frame with information about SNP locations, must have 3 columns - SNP name, chromosome, and position. See [sample SNP location file](http://www.bios.unc.edu/research/genomic_software/Matrix_eQTL/Sample_Data/snpsloc.txt).
* `genepos` – data frame with information about gene locations, must have 4 columns - the name, chromosome, and positions of the left and right ends. See [sample gene location file](http://www.bios.unc.edu/research/genomic_software/Matrix_eQTL/Sample_Data/geneloc.txt).

下面来看具体代码：

# source("Matrix_eQTL_R/Matrix_eQTL_engine.r");
library(MatrixEQTL)

## Location of the package with the data files.
base.dir = find.package('MatrixEQTL');
# base.dir = '.';

## Settings

# Linear model to use, modelANOVA, modelLINEAR, or modelLINEAR_CROSS
useModel = modelLINEAR; # modelANOVA, modelLINEAR, or modelLINEAR_CROSS

# Genotype file name
SNP_file_name = paste(base.dir, "/data/SNP.txt", sep="");
snps_location_file_name = paste(base.dir, "/data/snpsloc.txt", sep="");

# Gene expression file name
expression_file_name = paste(base.dir, "/data/GE.txt", sep="");
gene_location_file_name = paste(base.dir, "/data/geneloc.txt", sep="");

# Covariates file name
# Set to character() for no covariates
covariates_file_name = paste(base.dir, "/data/Covariates.txt", sep="");

# Output file name
output_file_name_cis = tempfile();
output_file_name_tra = tempfile();

# Only associations significant at this level will be saved
pvOutputThreshold_cis = 2e-2;
pvOutputThreshold_tra = 1e-2;

# Error covariance matrix
# Set to numeric() for identity.
errorCovariance = numeric();
# errorCovariance = read.table("Sample_Data/errorCovariance.txt");

# Distance for local gene-SNP pairs
cisDist = 1e6;

## Load genotype data

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);

## Load gene expression data

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);

## Load covariates

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);
}

## Run the analysis
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);

## Results:

cat('Analysis done in: ', me$time.in.sec, ' seconds', '\n');
cat('Detected local eQTLs:', '\n');
show(me$cis$eqtls)
cat('Detected distant eQTLs:', '\n');
show(me$trans$eqtls)

## Plot the Q-Q plot of local and distant p-values

plot(me)
因此，分析自己的数据需要准备
genotype
expression
covariates
gene location
SNP location
这五个文件，前三个需要每列的样本名对应且顺序一致。
作者也提供了生成模拟数据的代码：

# Create an artificial dataset and plot the histogram and Q-Q plot of all p-values
library('MatrixEQTL')

# Number of samples
n = 100;

# Number of variables
ngs = 2000;

# Common signal in all variables (population stratification)
pop = 0.2 * rnorm(n);

# data matrices
snps.mat = matrix(rnorm(n*ngs), ncol = ngs) + pop;
gene.mat = matrix(rnorm(n*ngs), ncol = ngs) + pop + snps.mat*((1:ngs)/ngs)^9/2;

# data objects for Matrix eQTL engine
snps1 = SlicedData$new( t( snps.mat ) );
gene1 = SlicedData$new( t( gene.mat ) );
cvrt1 = SlicedData$new( );
rm(snps.mat, gene.mat)

# Slice data in blocks of 500 variables
snps1$ResliceCombined(500);
gene1$ResliceCombined(500);

# name of temporary output file
filename = tempfile();

# Perform analysis recording information for
# a histogram
meh = Matrix_eQTL_engine(
snps = snps1,
gene = gene1,
cvrt = cvrt1,
output_file_name = filename,
pvOutputThreshold = 1e-100,
useModel = modelLINEAR,
errorCovariance = numeric(),
verbose = TRUE,
pvalue.hist = 100);
unlink( filename );
# png(filename = "histogram.png", width = 650, height = 650)
plot(meh, col="grey")
# dev.off();

# Perform the same analysis recording information for
# a Q-Q plot
meq = Matrix_eQTL_engine(
snps = snps1,
gene = gene1,
cvrt = cvrt1,
output_file_name = filename,
pvOutputThreshold = 1e-6,
useModel = modelLINEAR,
errorCovariance = numeric(),
verbose = TRUE,
pvalue.hist = "qqplot");
unlink( filename );
# png(filename = "QQplot.png", width = 650, height = 650)
plot(meq, pch = 16, cex = 0.7)
# dev.off();
阅读推荐：

生信技能树公益视频合辑：学习顺序是linux，r，软件安装，geo，小技巧，ngs组学！

B站链接：https://m.bilibili.com/space/338686099

YouTube链接：https://m.youtube.com/channel/UC67sImqK7V8tSWHMG8azIVA/playlists

生信工程师入门最佳指南：https://mp.weixin.qq.com/s/vaX4ttaLIa19MefD86WfUA

学徒培养：https://mp.weixin.qq.com/s/3jw3_PgZXYd7FomxEMxFmw

---------------------
作者：weixin_34295316
来源：CSDN
原文：https://blog.csdn.net/weixin_34295316/article/details/87552069
版权声明：本文为博主原创文章，转载请附上博文链接！

《微信小程序进阶实战之分答应用开发（中级项目）》（完整版）相关推荐

《学做程序经理》完整版
文/Joel Spolsky 译/罗小平指派一名优秀的程序经理,是团队产出优秀软件的重要前提之一.你的团队里可能没有这样的人,其实绝大多数团队都没有. Charles Simonyi,这位曾与 ...
coji小机器人_让孩子学做程序员 Coji编程机器人体验
原标题:让孩子学做程序员 Coji编程机器人体验玩具机器人可能每个孩子小时候的梦想,不过对于90后小编来说,小时候见过的所谓"机器人"好像最多也就是那种可遥控移动或放音乐的电子玩 ...
微软ping程序源代码完整版
微软ping程序源代码完整版编写自己的一个ping程序,可以说是许多人迈出网络编程的第一步吧!!这个ping程序的源代码经过我的修改和调试,基本上可以取代windows中自带的ping程序. 各个模 ...
php退款,php实现小程序退款完整版
本文主要和大家分享php实现小程序退款完整版,功能前提:1. 使用 wx php sdk (小程序支付完整版) , 2. 配置证书时使用绝对路径希望能帮助到大家. 1. 上代码:/** * 退款 * ...
微软ping程序源代码完整版(附详细的注释)
作者:侯志江单位:天津大学软件学院 E-mail :tjuhzjemail@yahoo.com.cn 日期:2005年1月1日内容简介: 编写自己的一个ping程序 ...
零基础学python-Python入门教程完整版（懂中文就能学会）
提取码:sjfo 目录大纲: 本套教程15天学前环境搭建 1-3 天内容为Linux基础命令 4-13 天内容为Python基础教程 14-15 天内容为飞机大战项目演练视频概括: 第一阶段(1 ...
杉德支付php代码实现_php实现小程序支付完整版
本文实例为大家分享了php实现小程序支付的具体代码,供大家参考,具体内容如下环境: tp3.2 + 小程序微信支付功能开通 Step1:下载PHP 支付SDK(下载地址) 放到Library\ ...
下载python教程-Python基础教程下载【黑马程序员完整版】
课程介绍目录大纲: 1-3 天内容为Linux基础命令 4-13 天内容为Python基础教程 14-15 天内容为飞机大战项目演练视频概括: 第一阶段(1-3天): 该阶段首先通过介绍不同领域 ...
用java制作一个软件控制小车_Android手机控制智能小车的手机端程序（完整版）...
[实例简介] 本程序是我写的Android手机控制智能小车的手机端的全部的源程序,下载后直接就能用. [实例截图] [核心代码] 624ba65e-a75e-4ba0-8e72-6dbc0823fcb ...
iso 软件测试控制程序,最新ISO17025：2017一整套程序文件完整版(共38个程序)
XXXX检测实验室 HHRBDC/CX01-01~38-2018 (3) 标题:保护客户机密信息和所有权的程序 (8) 标题:保证实验室诚信度的政策和程序 (10) 标题:质量手册的管理 (11) 标 ...

《微信小程序进阶实战之分答应用开发（中级项目）》（完整版）

《微信小程序进阶实战之分答应用开发（中级项目）》（完整版）相关推荐

最新文章

热门文章