Bioconductor通讯

发布的瓦莱丽Obenchain, 2016年1月

Bioconductor通讯每季回顾核心基础设施的发展、社区项目和未来方向。我们的目标是普遍感兴趣的主题以及对软件影响最大的主题。

在本问题中,我们将跟踪的发展InteractionSet分类和审查管理包存储库的技巧。Mike Love谈论构建基因表达实验的设计矩阵和Jim MacDonald带我们参观Bioconductor注释包。

内容

设计差异基因表达矩阵

Mike Love是一名博士后拉斐尔·伊的实验室在达纳·法伯癌症研究所和哈佛大学陈曾熙公共卫生学院的生物统计和计算生物学系工作。他开发了基因组学和表观遗传学的定量方法,教授edX课程,偶尔微博关于生物统计学和R.许多人知道他是非常受欢迎的作者和主要支持者DESeq2RNASeq差异基因表达数据包。

参观支持网站显示了他每天回答的问题的数量不仅与DESeq2包但是关于一般的基因表达分析。

在支持网站上许多与deseq2相关的帖子中,创建一个适当的设计矩阵是一个常规的帖子,似乎会引起相当多的困惑。下面Mike分享了他对哪些关键概念会导致最多问题的一些观察和想法。

关于“实验设计”的一些背景知识

实验设计指样品之间的相互关系,包括生物和实验信息(克隆1、处理B等)和技术信息(样品制备和加工的批次)。得到正确的实验设计——这发生在实验发生之前——是非常重要的,因为错误的实验设计会导致无法解释的数据。

最好的做法是在表格中跟踪实验设计(包括制剂批次),可以是CSV或TSV文件,也可以是Excel电子表格,当进行生物信息分析时可以导出到CSV。这是向计划中的实验或已经发生的实验的人解释实验设计的最简单的方法。

设计矩阵模型矩阵矩阵,在统计学中通常用X,这将用于统计建模。设计矩阵中的每一行描述一个样品,每一列提供关于该样品的信息,例如,样品是否经过处理,样品是在第1、2或3批次中,等等。对于设计矩阵的每一列,模型都有一个匹配系数,通常表示为β,以描述样本之间的表达差异。这些系数是对数尺度上的相加差异,因此RNA-seq计数或微阵列表达值的乘法差异(倍数变化),因此它们是日志褶皱变化.然后利用实验数据估计这些系数。

一个单一的实验设计并不意味着一个单一的设计矩阵,但往往涉及许多选择。例如,可以包含技术批次的系数(通常是一个好主意),也可以不包含,这将给出一个带有加法列的设计矩阵。设计矩阵编码了研究人员想要对样本做出的假设,并与感兴趣的生物问题有关。在显著性检验中,生物问题被表述为零假设,一个或多个系数等于零,结果是ap值.只有假设是合理的,并且数据的模型很好地指定了,p值才是对概率的有意义的估计。

我应该说,我是从教科书(约翰·赖斯的)中学到这些话题的数理统计与数据分析“,桑福德,韦斯伯格的应用线性回归,Bioconductor书),以及从阅读大量的帖子Bioconductor沃尔夫冈·休伯、戈登·史密斯、西蒙·安德斯、詹姆斯·麦克唐纳、亚伦·伦等人的支持论坛。

情况与控制

简单的设计似乎不会造成太大问题。例如,对照和处理样本,或对照,处理1和处理2。使用它们可以很容易地建模R的内置公式而且model.matrix函数,然后输入到limma刨边机或其他Bioconductor包。DESeq2直接把公式表达式并在内部转换为设计矩阵。

混杂和批处理效应

有时,定量/计算问题以错误信息的形式出现,这表明实验设计中存在固有问题。其中之一就是兴趣的比较抱愧蒙羞配合技术因素,如样品制备批次。当对照和处理样品各自分批制备时,存在混杂的典型情况,但也常见实验设计不良的情况,如:

条件 批处理
控制 1
控制 1
治疗。一个 1
治疗。一个 1
治疗。B 2
治疗。B 2
治疗。C 2
治疗。C 2

虽然处理A可以与对照进行比较,处理C可以与处理B进行比较,但不能在批次之间进行比较。有些比较无法进行的原因是由于基因表达的差异,例如处理B与对照的效果,无法从不同样品制备批次之间可能产生的差异中分离出来。这里最有效的解决方案是使用块设计,其中每个批包含所有可能的条件。每个批次至少应该包括对照样品,以便利用这些样品来估计批处理效果。

这两个环节解释了为什么批处理效应对高通量实验(或任何实验)构成了一个大问题:

阻塞、交互和嵌套设计

块实验设计和其他的实验设计,例如测试条件之间交互作用的重要性,或嵌套交互作用,可以在limma用户指南的章节中阅读单通道实验设计

该指南详细描述了如何用不同的方式表述设计矩阵,以回答相同的问题,并解释了不同的参数化如何影响对结果的解释。limma作者推荐的方法可以应用到其他方面Bioconductor包。

先进的设计

还有一些非常复杂的设计,其中包含许多技术和生物因素,研究人员需要进行很多比较,但却不知道如何进行比较。在这种情况下,对于那些发现自己不知道使用什么设计或如何解释系数的人,我强烈建议他们考虑与当地统计学家或具有线性建模或定量分析背景的人合作。

解释定量分析是一件很困难的事情Bioconductor在很大程度上简化了高通量分析的分析,期望复杂的结果可以由没有定量背景的人编译或解释不一定合理。我认为找一个能够提供帮助的合作者更安全、更合理,如果他们从项目一开始就被纳入其中,这样的合作者可以帮助识别实验设计中的问题。

回到顶部

开始使用Bioconductor注释包

吉姆·麦克唐纳(Jim MacDonald)是华盛顿大学环境与职业健康科学系的生物统计学家。他分析了从表达(微阵列,RNA-Seq)到基因组(SNP阵列,DNA-Seq, ChIP-Seq,甲基化阵列,BS-Seq)和其他“组学”数据的HTS数据域。他一直在大力参与指导Bioconductor项目,并贡献和维护大量的注释包。

在2015年10月发布的版本中,由于布法罗的搬迁,我们失去了一些员工,导致人手短缺。吉姆介入并负责所有的内部建设Bioconductor注释包。Jim对注释世界的全面理解可以从他在支持网站.在本节中,我们联合起来(90% Jim, 10% Val)概述了一些关键包,以及如何使用它们来回答一些常见的分析问题。

主要包

本节重点介绍使用最频繁的部分Bioconductor注释包。

值得注意的是,一些注释包与特定的基因组构建相关联,而另一些则不是。的TxDb家族包含基因/转录本/外显子等的位置。基于给定的构建。BSgenome而且SNPlocs是特定于构建的包的其他示例。因为基因组组装需要将整个基因组的结构拼凑在一起,所以新的构建每隔几年才会发布一次。特定于构建的注释包中的数据可以相当稳定,并且可以持续多年。

其他软件包与基因在哪里被发现无关,因此与基因组构建无关,例如OrgDb家庭。这些软件包可以被认为是我们所掌握的关于某一特定日期某一特定有机体基因的所有信息的封装,因为我们知道,这些信息可能在下一周就会过时,至少部分会过时。它们包含RefSeq、GenBank或UniGene id等信息,表示临时转录本。这些工作正在进行中,正在根据公众的意见不断更新和修改。Bioconductor在发布时每6个月更新一次这些包。的AnnotationForge包提供了用于构建自己的函数OrgDb(或其他包)如果你想要更新潮的东西。

OrganismDb包包含特定于构建的包和非特定于构建的包的组合,这可能不适合您的用例。但是,很容易切换TxDb方法获取更合适的版本TxDb < -函数。

常见的任务

在考虑特定的任务之前,我们应该首先讨论如何确定输入()和输出()可用于特定的注释包。的keytypes函数将返回可用作输入的所有类型的注释。作为一个例子,让我们使用org.Hs.eg.db包中。

>库(org.Hs.eg.db) >键类型(org. hs .db) [1] "ACCNUM" "ALIAS" " integrbl " " integrblprot " " integrbltrans " [6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME" [11] "GO" GOALL" IPI" MAP" "OMIM" [16] "ONTOLOGY" ONTOLOGYALL" PATH" PFAM" PMID" [21] "PROSITE" "REFSEQ" SYMBOL" "UCSCKG" "UNIGENE" [26] "UNIPROT"

我们可以列出所有可用的对于一个给定的keytype使用函数。

> head(keys(org. hs . e.g. .db)) [1] "1" "2" "3" "9" "10" "11" > head(keys(org. hs . e.g. .db, "ENSEMBLPROT")) [1] "ENSP00000263100" "ENSP00000470909" "ENSP00000323929" "ENSP00000438599" [5] "ENSP00000445717" "ENSP00000385710"

我们可以得到所有可用的,或注释数据,我们可以映射出现。

>列(org.Hs.eg.db) [1] "ACCNUM" "ALIAS" " integrbl " " integrblprot " " integrbltrans " [6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME" [11] "GO" GOALL" IPI" MAP" "OMIM" [16] "ONTOLOGY" ONTOLOGYALL" PATH" PFAM" PMID" [21] "PROSITE" REFSEQ" SYMBOL" "UCSCKG" "UNIGENE" [26] "UNIPROT"

将制造商id映射到基因符号

一个常见的任务是通过将制造商的ID映射到更一般的东西来注释微阵列实验,例如HUGO基因符号,或NCBI (gene, GenBank, RefSeq, UniGene)或ensemble (ensemble基因,ensemble转录本)ID。例如,我们可以将人类基因1.0 ST数组中的Affymetrix ID映射到相应的HUGO符号。

>库(hugene10sttranscriptcluster.db) > hugene <- hugene10sttranscriptcluster.db ##最小化typing > select(hugene, "8012257", "SYMBOL")'select()'返回键和列之间的1:1映射

这是一个非常简单的示例,除了作为一个快速交互查询之外,可能没有什么用处。注意,我们没有指定keytype论点。默认的keytype理由选择是正在使用的包的中心ID。在本例中,它是PROBEID,由于我们使用的是探测集id,没有必要指定keytype

的一个更常见的用例是注释向量,返回一个或多个输出id(或).作为示例,我们将只使用5个从hugene10sttranscriptcluster.db包中,查询HUGO符号和Entrez基因id。

> ids <- keys(hugene)[15000:15005] > ids [1] "8005171" "8005191" "8005200" "8005202" "8005204" > annot <- c("SYMBOL","ENTREZID") > select(hugene, ids, annot)PROBEID SYMBOL ENTREZID 1 8005171 TRPV2 51393 2 8005191 LRRC75A-AS1 125144 3 8005191 SNORD49A 26800 4 8005191 SNORD49B 6920106 6 8005200 LRRC75A-AS1 125144 7 8005200 SNORD49A 26800 8 8005200 SNORD49B 692106 9 8005200 SNORD49A 26800 8 8005200 SNORD49B 692106 10 8005202 LRRC75A-AS1 125144 12 8005202 SNORD49B 692087 13 8005202 SNORD65 692106 14 8005204 CCDC144A 9720 15 8005204 CCDC144CP 348254 16 8005204 CCDC144B284047 17 8005204 ccdc144nl 339184 18 8005204 loc101929141 101929141 19 8005221 < na > < na >

关于上述结果,请注意三点。首先,返回的PROBEID列data.frame与输入id的顺序相同。第二,一些Affymetrix id与多个基因对应。返回所有的映射,并返回一个消息,其中一些映射为1:many.由于1:多个映射,返回的维数data.frame不匹配我们想要注释的数据的维度(例如,我们想要五个id的信息,却返回了19行数据)。第三,如果其中之一没有注释(最后一个),anNA返回值。

如果我们想要保证返回的数据是相同的顺序而且和输入长度相同吗我们可以用向量mapIds代替。然而,mapIds只能做一个keytype,并返回一个向量列表而不是一个data.frame.不像选择,它对键类型有默认值,mapIds需要第四个参数,指定keytype我们正在使用。

> mapIds(hugene, ids, "SYMBOL", "PROBEID") 8005171 8005191 8005200 8005202 8005204 "TRPV2" "LRRC75A-AS1" "LRRC75A-AS1" "SNORD49A" "CCDC144A" 8005221 NA

我们可以很容易地将其包装在一个小脚本中以返回一个data.frame每个只有一行关键

> d.f <- as.data.frame(lapply(annot, function(x) mapIds(hugene, ids, x, "PROBEID")) > names(d.f) <- annot > d.f SYMBOL ENTREZID 8005171 TRPV2 51393 8005191 LRRC75A-AS1 125144 8005200 LRRC75A-AS1 125144 8005202 SNORD49A 26800 8005204 CCDC144A 9720 8005221  

的默认值mapIds是取任意1:多个映射的第一个实例。这对于某些用例(例如RefSeq ID)是没问题的,但是对于其他情况(例如GO ID)就不那么有用了,因为我们希望返回所有的值。我们可以用multiVals参数来控制返回的内容。请注意,这个参数后面有一个省略号(...)参数,因此不能使用位置参数,而必须指定multiVals直接的论点。

> mapIds(hugene, ids, "SYMBOL", "PROBEID", multiVals = "list") $ ' 8005171 ' [1] "TRPV2" $ ' 8005191 ' [1] "LRRC75A-AS1" "SNORD49A" "SNORD49B" "SNORD65" $ ' 8005200 ' [1] "LRRC75A-AS1" "SNORD49A" "SNORD49B" "SNORD65" $ ' 8005204 ' [1] "CCDC144A" "CCDC144CP" "CCDC144B" "CCDC144NL" "LOC101929141" $ ' 8005221 ' [1] NA

如果我们希望为注释使用矩形格式,即保留所有的1:many映射,同时确保每一行直接映射到表达式值数组,则可以使用DataFrame相反,告诉mapIds返回一个CharacterList

> lst <- lapply(annot, function(x) mapIds(hugene, ids, x, "PROBEID", multiVals = "CharacterList") > d.f <- as(lst, "DataFrame") > names(d.f) <- annot > d.f DataFrame with 6行2列SYMBOL ENTREZID   8005171 TRPV2 51393 8005191 LRRC75A-AS1,SNORD49A,SNORD49B,…125144, 26800, 692087,……8005200 LRRC75A-AS1, SNORD49A, SNORD49B,…125144, 26800, 692087,……8005202 SNORD49A, LRRC75A-AS1, SNORD49B,…26800, 125144, 692087,……8005204 CCDC144A, CCDC144CP, CCDC144B,…9720, 348254, 284047,……8005221 NA NA

将Entrez基因ID映射到TRPV2染色体位置

根据以上数据,也许我们对TRPV2很感兴趣,想知道它的染色体位置。我们可以用Homo.sapiens包来获取该信息。虽然可以使用该基因的HUGO符号来获得位置,但使用Entrez基因ID更好,因为它更可能是唯一的。

>选择(Homo。智人,“51393”,c(“TXCHROM”、“TXSTART”,“TXEND”),“象征”)'select()'返回键和列之间的1:1映射ENTREZID TXCHROM TXSTART TXEND 1 51393 chr17 16318856 16340317

这只是告诉我们文本的开始和停止位置。如果我们想要外显子位置,我们也能得到。

>选择(Homo。智人," TRPV2", c("EXONCHROM","EXONSTART","EXONEND"), "SYMBOL")SYMBOL EXONCHROM EXONSTART EXONEND 1 TRPV2 chr17 16318856 16319147 2 TRPV2 chr17 16320876 16321182 3 TRPV2 chr17 16323429 16323562 4 TRPV2 chr17 16325913 163235203 5 TRPV2 chr17 16326783 16327081 6 TRPV2 chr17 16329413 16330583 7 TRPV2 chr17 16330763 16330861 9 TRPV2 chr17 163301631 163322901 10 TRPV2 chr17 16332131 16332296 11 TRPV2 chr17 16335098 16335164 12 TRPV2 chr17 16335280 16335614 13 TRPV2 chr1716336888 16337012 14 TRPV2 chr17 16338204 16338283 15 TRPV2 chr17 16340103 16340317 16 TRPV2 chr17 16325913 16326207 17 TRPV2 chr17 16330177 16330191 18 TRPV2 chr17 16330766 16330861

虽然这对单个基因是有用的,但对于大量的基因来说,它可能会变得笨拙。我们可以用transcriptsByexonsBy功能与TxDb.Hsapiens.UCSC.hg19.knownGene软件包,一次获取所有基因的信息,子集到我们关心的那些。

> trscpts <- transcriptsBy(TxDb.Hsapiens.UCSC.hg19)。> trscpts[["51393"]]具有2个范围和2个元数据列的GRanges对象:seqnames范围strand | tx_id tx_name    |   [1] chr17 [16318856, 16340317] + | 60527 uc002gpy。3 [2] chr17 [16318856, 16340317] + | 60528 uc002gpz。4 > exns <- exonsBy(TxDb.Hsapiens.UCSC.hg19)。> exns[["51393"]]具有18个范围和2个元数据列的GRanges对象:seqnames范围strand | exon_id exon_name    |  [1] chr17 [16318856, 16319147] + | 217024  [2] chr17 [16320876,16321182] + | 217025  [3] chr17 [16323429,16323562] + | 217026  [4] chr17 [16325913,16326207] + | 217027  ... ... ... ... ... ... ...[14] chr17 [16335098, 16335164] + | 217037  [15] chr17 [16335280,16335614] + | 217038  [16] chr17 [16336888, 16337012] + | 217039  [17] chr17 [16338204,16338283] + | 217040  [18] chr17 [16340103,16340317] + | 217041  ------- seqinfo: 93个序列(1个循环)来自hg19基因组

使用*范围对象超出了本通讯的范围,因此我们不再进一步探讨它们。有关更多信息,请参阅IRanges维格内特,还有农庄小插曲。

回到顶部

可再生的研究

管理包版本biocLite ()

Bioconductor遵循一年两次的时间表,一个在春天,一个在秋天。R每年只发布一个主要版本,通常在春季。因为每个Bioconductor发布与版本绑定R这种不对称的时间表造成了一些混乱。

当发布版本与Spring一致时,两者的开发分支R而且Bioconductor成为发布分支。在接下来的6个月里,包在两个Bioconductor' devel '和release分支是基于的' release '版本构建的R

在秋天,Bioconductor有一个版本R没有。的Bioconductor' devel '分支成为当前的' release ',并使用的发布版本R.新Bioconductor' devel '分支使用的' devel '版本R.建筑的目的Bioconductor“重击”对Rdevel是为了允许秋季的平稳过渡,具体来说,它允许Bioconductor' release '分支始终与R“释放”分支。

BiocInstallerPackage有几个功能来帮助管理干净的“发布”和“开发”包存储库。下面是一些常见的安装和版本不匹配问题的故障排除提示。

关于保持版本同步的更多信息可以在2021年欧洲杯比分预测 部分的网站。

回到顶部

基础设施

InteractionSet

艾伦·伦,莉兹·英-西蒙斯和马尔科姆·佩里正在拍摄一部InteractionSet用于存储和操作cia - pet和Hi-C实验数据。

ChIA-PET是染色质相互作用分析与配对末端标签。这些实验探索了与一些感兴趣的蛋白质有关的全基因组相互作用。该技术区别于Hi-C的关键步骤是抗体驱动的免疫沉淀步骤,以富集与特定蛋白质结合的染色质。染色质的相互作用只能被确定为基因组中与感兴趣的蛋白质有结合位点的部分。相互作用网络可以解释转录因子,绝缘体蛋白或转录机制。一个cia - pet实验提供了蛋白质在构建3D基因组组织中的潜在作用的信息。

Hi-C方法通过在全基因组尺度上识别长范围染色质相互作用,提供了关于三维基因组结构的信息。这些数据通常用于研究基因组结构的各个方面,如染色体区域、拓扑结构域、开放/封闭室室和染色质结构。

来自这两种技术的数据使对基因组区域之间的物理相互作用的研究成为可能。的InteractionSetPackage提供类来表示这些交互并存储相关的实验数据。其目的是为包开发人员提供稳定的类定义,可以通过大量的方法进行2021欧洲杯体育投注开户操作。

这个包定义了以下类:

类有排序和重复检测的方法;用于执行一维或二维重叠;用于计算相互作用的轨迹之间的距离;以及计算相互作用组的最小边界框。方法也可用于在类之间转换,或转换为标准Bioconductor像一个对象RangedSummarizedExperimentGRangesList

回到顶部

新功能R/Bioconductor

新增功能R(3.3)和Bioconductor本季度(3.3):

回到顶部

项目活动

认识到社区的贡献

欧洲Bioconductor2021欧洲杯体育投注开户开发人员会议去年12月,一个奖项被颁发给了对……做出贡献的个人Bioconductor支持网站论坛。

该奖项由F1000研究资助,该研究最近推出了一项专门的Bioconductor通道.的条款“对支持网站做出最大贡献的人”和“参加欧洲开发者大会的人”获得了该奖项。

恭喜获奖者Aaron Lun和Michael Love!每位获奖者都因发表在《F1000》杂志上的一篇文章而获得了出版费减免的奖励Bioconductor通道。其他有重要贡献的贡献者有Jim MacDonald, Gordon Smyth, Ryan Thompson和Steve Lianoglou。感谢所有花时间回答问题并在支持网站上分享他们经验的人。

感谢Mark Dunning和Laurent Gatto建议设立这个奖项(并组织了这次会议!),感谢Thomas Ingraham和F1000 Research的赞助。

回到顶部

10月发布

Bioconductor3.2于10月14日发布,包含1104个软件包,257个实验数据包,917个注释包。有80个新的软件包。

请注意这是最后的版本Bioconductor在雪豹上支持。雪豹用户应该计划在2016年春季下一个版本发布之前迁移到Mavericks或更新版本。

在这个版本中包含了80个新软件包。软件包摘要和官方发布时间表可以在网站

回到顶部

网站流量

以下是2015年第四季度(11月1日- 12月28日)的会话数和新用户数与2014年第四季度的比较。

2015年第四季度的会话数(18.05%)、用户数(11.27%)和页面浏览量(12.11%)均较2014年有所增长。平均浏览页面数(-5.39%)、平均会话持续时间(-0.54%)、跳转率(即单页访问;-3.30%)和新课程百分比(-8.86%)。

大多数用户仍然使用台式电脑(和笔记本电脑),但移动用户的数量正在稳步增长。在2015年第四季度,移动设备的新用户数量增长了39.8%,而桌面设备的新用户数量增长了6.7%。

网站流量:按设备划分的新用户
桌面(包括笔记本电脑)
2015年11月1日- 2015年12月27日 71111例(92.86%)
2014年11月1日- 2014年12月27日 66622例(93.96%)
变化百分比 6.74%
移动
2015年11月1日- 2015年12月27日 4230例(5.52%)
2014年11月1日- 2014年12月27日 3026例(4.27%)
变化百分比 39.79%
平板电脑
2015年11月1日- 2015年12月27日 1241例(1.62%)
2014年11月1日- 2014年12月27日 1260例(1.78%)
变化百分比 -1.51%


统计数字由谷歌分析

回到顶部

包下载和新提交

2015年10月、11月和12月软件包的唯一IP下载量分别为40085次、41499次和34216次。2014年同期的数字分别为44593、32728和30622。数据必须按月比较(与季度比较),因为有些ip在月份之间是相同的。的完整摘要,请参阅网站下载数据

2015年第四季度共增加17个软件,1个注释,5个实验数据包。

Bioconductor版本 软件 注释 实验数据
3.2 1104 895 257
3.3 1121 896 262

回到顶部

即将来临的事件

看到活动页面查看所有课程和会议的列表。

确认

感谢吉姆麦克唐纳和迈克爱的贡献部分,亚伦伦打样InteractionSet节和Bioconductor核心团队负责编辑审查。

回到顶部

向瓦莱丽发送评论或问题valerie.obenchain@roswellpark.org