eQTL教程，Bioc 2014

VJ Carey Stvjc在Channing Dot Harvard Dot Edu

基本数据

我们将使用通过RNA-seq获得的转录谱，应用于HapMap细胞系(Montgomery, Sammeth, Gutierrez-Arcelus, Lach, Ingle, Nisbett, Guigo和Dermitzakis(2010))。报告的基因水平计数数据已在重新计算项目(Frazee, Langmead, and Leek(2011))中使用基因的ENSEMBL命名法表示为ExpressionSet (montpick_set . rdata)。对这些计数进行过滤，以排除所有样本中具有一致零值的基因，进行rlog转化(DESeq2)，然后过滤，只包括在基因MADs总体分布的第40百分位上有MAD的基因。

在此处给出了转换转换的特征代码。您不需要执行此操作，因为转换数据可用于加载，如下所述。（exs2se（）是此小插图源中可用的一系列自定义代码。）

库（BioBase）加载（“/ data / eqtl2014 / montpick_eset.rdata”）mp = montpick.eset注释（mp）=“org.hs.eg.db”tmp = select（org.hs.eg.db，keytype =“ensembl”，columns =“EntrezID”，键= FeatureNames（MP））EZS = split（TMP [，2]，TMP [，1]）EZS = Sapply（EZS，“[”，1）删除=哪个（重复（EZS）|是。）]）mpse = exs2se（mp，annodb = org.hs.eg.db，probekeytype =“EntrezID”）

转换概述的表达数据

我们获得如下变换计数数据的基本形象。

load（“/ data / eqtl2014 / rmpdepf.rda”）＃rlogged montgomery pickrell deseqdataset，一些正，过滤的rmpdepf

##类：摘要分析## DIM：7054 129 ## EXPTDATA（1）：MIAME ##测定（1）：''## Rownames（7054）：8813 57147 ... 56099 56104 ## Rowdata元数据列名称（6）：底座Basevar ... Dispfit ## Rlogintercept ## Colnames（129）：Na06985 Na06986 ... Na19239 Na19257 ## Coldata名称（5）：sample.id num.tech.reps人口研究## solductor

Amazon S3桶中的基因型数据

1000个基因组项目中生成的所有数据都可以在S3（简单存储服务）容器中提供。我们有一个URL，可用于访问VCF格式的CHR22的基因型。

加载（“/ data / eqtl2014 / s3url.rda”）s3url

## [1]“http://1000genomes.s3.amazonaws.com/release/20110521/all.chr22.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz”

图书馆（VariantAnnotation）

head = scanvcfheader（s3url）str（head，max.level = 4）

##正式类'VCFHEADER'[包“VariantAnnotation”]，3个插槽## .. @参考：CHR（0）## .. @样本：CHR [1：1092]“HG00096”“HG00097”“HG00099”“hg00100“... ## .. @ header：正式的'simpledataframelist'[包”讽刺“]用4个插槽## .. .. @ ElementType：Chr”DataFrame“## .. .. @ElementMetadata：null ## .. .. .. @ metadata：list（）## .. .. .. @ listdata：4

一个天真的分析

GGTools包包括功能，Cisassoc（），可以在陷阱索引的VCF和概括的表达数据中进行CIS中的CIS中的添加剂遗传关联。该测试是广义线性模型的分数测试适合具有连续响应的独立样品，响应方差大致独立于响应均值。这些假设通常不适用于来自RNA-SEQ实验的计数数据，但它们可能在RLOG转化后约束地保持。在执行一些测试后，我们将检查违规行为的严重程度。

这是计算的代码，在CHR17上的200个基因，关联检验统计量及其在基因型表达置换下的版本。不要执行这些计算，它们将花费教程太长时间。教程数据提供了答案，请参见下面的内容。

图书馆（GGTools）

## Loading required package: GGBase ## Loading required package: snpStats ## Loading required package: survival ## Loading required package: splines ## Loading required package: Matrix ## ## attach package: 'Matrix' ## ##以下对象被'package:VariantAnnotation'屏蔽:## ## expand ## ##以下对象被“package:IRanges”屏蔽:## ## expand ## ##加载所需的package: data。表# #数据。table 1.9.2 For help type: help("data.table") ## KernSmooth 2.23 loaded ##版权M. P. Wand 1997-2009 ## ## attach package: 'GGtools' ## ##以下对象被'package:stats'屏蔽:## ## getCall . cfg . cfg

库（DoparPristial）RegisterDoseQ（）R17 = RMPDEPF [哪个（SEQNAMES（RMPDEPF）==“CHR17”），] SET.SEED（1234）S3URLB = SUB（“22”，“17”，S3URL）C17= Cisassoc（R17，TabixFile（S3URLB），CISRADIUS = 100000，LBMAF = .025）C17D = AddMindist（C17）

这是输出中的提取物。我们手动添加了一个柱子“Mindist”，其中核苷酸的数量来自SNP到最近的基因体末端。“Mindist”为SNP躺在基因体内为零。

c17d [1:3, 1:5)

##带有3个范围和5个元数据专栏的Granges：## SEQNAMES范围股票|paramRangeid ref ##    |<因子>  ## [1] 17 [36926181,66926181] * |3927 T ## [2] 17 [36926303,36926303] * |3927 C ## [3] 17 [36926384,36926384] * |3927 C ## Alt ChiSQ PermScore_1 ##    ## [1] 2.9869 0.51171 ## [2] T 4.3823 1.01236＃3] T 0.4235 0.05441 ## --- ## SEQLENGTH：## 17 ## na

名称（MCols（C17D））

## [1]“paramRangeid”“ref”“alt”“alt”“chisq”## [5]“permscore_1”“permscore_2”“permscore_3”“snp”## [9]“maf”“probeId”“得分”“fdr“## [13]”Mindist“

计算插件fdr

下面的算法来自统计学习的要素由Hastie，Tibshirani和Friedman（Springer 2001）。

我们可以使用GGTools中的PIFDR（）函数来获得上面测试的Gene-SNP对的估计FDR。

FDR = PIFDR（C17D $ CHISQ，C（C17D $ PERMSCORE_1，C17D $ PERMSCORE_2，C17D $ PERMSCORE_3））总和（FDR == 0）

## [1] 430

总和（fdr <= .05）

## [1] 5292

敏感性分析，以及“算什么？”

以下计算产生了灵敏度分析显示，面对显着相关的SNP-基因对的枚举。

库(data.table)库(foreach)

## foreach：旋转分析的简单，可扩展并行编程##使用革命性r来可扩展性，容错等。## http://www.revolutionAnalytics.com.

registerDoseQ（）dt1 = data.table（as.data.frame（mcols（c17d）））

s1 = eqsens_dt（dt1）

Plotsens（S1，Ylab =“给定FDR的重要SNP-基因对计数”）

##加载所需的包:ggplot2

块DOPL

如果代替计数SNP-基因对，我们希望计数显示与某些基因型相关联的基因，该插入FDR程序可用于每个基因的最大关联评分，用于观察和置换的测试。

s2 = eqsens_dt（dt1，by =“probes”）

Plotsens（S2，Ylab =“鉴于FDR的eqtl证据的基因计数”）

块PL2的情节

这里有一个指示，如果我们在MAF和距离比我们在初始运行中更急剧过滤，我们会获得更大的产量。

练习1。

我们如何计算FDR用于基因明假设H0G：基因G的平均表达与所有SNP CIS到G的SNP含量无关，在CIS半径小于50kb和MAF大于5％的情况下，大于5％？

可视化

以下代码为所选SNP-基因对创建基因型表达关联的可视化。

plotone = function（summex，vcf.tf，suppex，thath，genome =“hg19”，......）{sampidsinsumm = colnames（summex）sampidsinvcf = sampsinvcf（vcf.tf）Oksamp =相交（Sampidsinsumm，SampidsInvcf）StopIfnot（长度（OKSAMP）> 0）SUMMEX = SUMMEX [，OKSAMP] VP = SCANVCFPARAM（固定=“ALT”，INFO = NA，GENO =“GT”，SAMPLES = OKSAMP，= OKSAMP）VDATA = READVCF（VCF.TF，基因组= genome，param = vp）rdd = rowdata（vdata）vdata = gdatools ::: snvsonly（vdata）gtdata = fing（geno（vdata）[[1]]图（As.numeric（As.numeric（Marmatays（Summex）[[1] [1]]] [airsa，]）〜gtdata，...）} r17 = rmpdepf [哪个（as.character（seqnames（rmpdepf））==“chr17”），] s3urlb = sub（“22”，“17”，S3URL）Plotone（R17，TabixFile（S3URLB），1，C17D [1]）

块的dooneplot

练习2。

这里描绘的SNP和基因的名称是什么？Replot，用此信息注释（改善x和y标签以获得最大的解释清晰度）。

练习3。

显示的关联的FDR是什么?

我们可以简单地表现出四个关系。

opar = par（nogreadonly = true）par（mfrow = c（2,2））for（i在2：5）plotone（r17，tabixfile（s3urlb），i，c17d [i]）

块块DOFFF

Par（opar）

练习4。

描绘四个基因-SNP对具有非常小的FDR的关联。四个基因应该是不同的。

暗示：

order(fdr)[1:10] c17d[7659] args(plotOne) plotOne(r17, TabixFile(s3urlb)， "114881"， c17d[7659])

练习5：主监管机构。

找到与多个不同基因相关联的FDR的SNP少于一个百分点。可视化关系。

与“gwas命中”的连接

我们可以按如下方式获取当前NHGRI GWAS目录的图像。

library(gwascat) curgw = makeCurrentGwascat() curgw

然而，截至2014年7月28日，当前目录使用的坐标是GRCh38。我们使用rtracklayer的liftOver和适当的链文件来获得hg19寻址方案，丢失了少量不可映射的轨迹。

加载（“/ data / eqtl2014 / curgwhg19.rda”）curgwhg19 [1：3,1：5]

# #农庄与3和5元数据列:# # seqnames范围链| Date.Added.to.Catalog # # < Rle > < IRanges > < Rle > | <人物> # # 1 chr19(7739177、7739177)* | 04/16/2014 # # 2 chr2(191113009、191113009)* | 07/26/2014 # # 3 chr6(31197514、31197514)* | 07/26/2014 # # PUBMEDID第一。Author Date Journal ##     ## 1 24123702 Chung CM 03/03/2014 Diabetes Metab Res Rev ## 2 24325914 Chu M 12/09/2013癌变## 3 24324648 Yatagai Y 12/04/2013 seqlength# # chr19 chr16 chr1 chr17 chr21 chr13…chr15 chr9 chr12 chr14 chrX ## NA NA NA NA NA NA…呐呐呐呐呐

练习6.

在FDR 5％或更低的SNP展示EQTL协会和GWAS击中时，SNP与GWAS频率如何频繁将此频率与SNP表现出在FDR 95％或更高的情况下进行比较。

进一步练习

练习6：原产地群体。

我们在这分析中忽视了原产地群体。使用可视化评估基因OSBPL7和SNP RS17774008之间的关联是否在来自CEU和YRI群体的受试者中具有类似的形式。

您可以使用S3URL获取访问CHR22的基因型数据。在CHR22基因和SNP上运行Cisassoc需要5分钟左右。单独获得CEU和YRI组的分析。关于FDR 5％左右的SNP-Gene对的重要性是什么协调？

练习7：用计数表达值有效分析。

/data/eqtl2014/montpick_eset.rdata具有原始计数。（注意，探针标识符位于EnsemBl词汇表中。）使用DESEQ2或EDGER获取所选SNP-基因对的FDR。写一个为任何给定选择而执行此操作的函数。原则上，需要为每个基因SNP对单独建模平均方差关系，因此获得常规搜索的统计最佳解决方案可能是费力的。在进行eqtl的高斯的分析之前，了解使用盲目和快速设置为true的使用rlog的成本很好，而不是对eqtl的配对特定的最佳建模的负二进制测试。

参考书目

A. C. Frazee, B. Langmead和J. T. Leek。重新计票:一个可供分析的RNA-seq基因计数数据集的多实验资源。在：BMC生物信息学12.1（2011年1月），p。449.Sisn：1471-2105。DOI：10.1186 / 1471-2105-12-449。。

[2] S. B. Montgomery，M. Sammeth，M.Gutierrez-Arcelus等。“在高加索人群中使用第二代测序的转录组遗传学。”在：自然464.7289（2010年4月），第773-7页。ISSN：1476-4687。DOI：10.1038 / Nature08903。。