跳至主要内容

Bioconductor 简介

学习如何对您的生物实验结果进行计算与统计分析。
· 8分钟 阅读

使用 AI 探索

在本教程中,您将熟悉 Bioconductor 生态。借此,您将能够对生物实验结果进行计算与统计分析,因为任何研究者都需要用数据证明其结论的显著性。换句话说,您将能够借助该项目判断所用数据是否与您的研究假设一致。此外,您还能以统计方式验证科学问题,而无需耗费过多时间与资金。

简而言之,本教程将厘清:

  • 什么是 Bioconductor?
  • Bioconductor 的 R 包代表什么?
  • 我们如何使用它?

简要概览

Bioconductor:高通量基因组数据的分析与理解

Bioconductor 是用于计算生物学与生物信息学的开放软件开发平台。它是一个面向计算生物学社区的开源、开放开发的软件项目。所谓开源,是指任何人都可以阅读和修改底层代码;所谓开放开发,是指任何人都可以参与代码的贡献与开发。它使用 R 语言——一门在数据科学领域广泛使用的语言,因其在数据分析方面灵活易用。R 还具备高质量图形能力,并且可以轻松与其他语言协作。Bioconductor 为基因组数据分析提供强大的统计与图形方法,同时便于在实验数据分析中整合 GenBank、GO、LocusLink 与 PubMed 等生物学元数据。

此外,它支持快速开发可扩展、可互操作且可伸缩的软件。借助 Bioconductor 还能推动高质量文档与可重复研究的发展,并为计算与统计方法的培训提供支持。简要地说,Bioconductor 可分为以下三个部分:

统计分析理解高通量
大规模数据生物学语境测序
技术伪影可视化微阵列
设计实验可重复性流式细胞术

计算生物学中的分析几乎总是需要将许多不同的工具组合起来以回答特定问题。通常这会形成可供其他研究者使用的包。因此,Bioconductor 还是一个带有一定规则与指导原则的 R 包软件仓库。3.7 版本拥有超过 1500 个包,全球有成千上万的研究者在使用。这些包通常配有优秀的文档。由于人人皆可贡献,针对同一问题往往会出现多种包可选,从而在包的开发者之间形成良性竞争。

官网提供如何安装、学习、使用以及回馈 Bioconductor 项目的资料。该项目始于 2001 年秋,目前在美国、欧洲与澳大利亚拥有 20 多位核心开发者。

bioconductor

概括而言,Bioconductor 专注于高通量基因组数据的统计分析与理解。之所以强调统计分析,是由数据的特性决定的。现代方法产生的海量数据必然需要统计性概括。作为研究者,我们必须在统计框架下开展聚焦、良好设计的实验。生成现代遗传数据的技术以及后续使用的实验室流程会引入技术伪影,需要通过统计技术加以处理。最后,紧凑的数据整合挑战也促使我们谨慎思考各步骤的统计后果。

Bioconductor R 包

Bioconductor 软件由 R 的扩展包组成。R 包是用于执行特定类型分析(如 affy、cluster、graph 包)的代码(R、C 或其他)、文档和/或数据的结构化集合。它们提供特定统计与图形方法的实现。

在最新的 Bioconductor 发行版中,已有 1500 余个包。它们大致可分为四类:

• 分析类软件包:用于分析高通量基因组数据的工具,例如 IMMAN、limma

• 注释类包:包含标识符映射、基因模型、通路等的静态数据库;例如 TxDb.Hsapiens.UCSC.hg19.knownGene

• 示例实验数据包:用于演示软件功能的数据集,例如 airway

• 工作流程包:描述生物信息学工作流程的文档,这些流程涉及多个 Bioconductor 包以阐释蛋白质组学等生物学领域,例如 highthroughputassays

主要的 Bioconductor 包

Bioconductor 的使用

安装

首次使用 Bioconductor,需先获取最新的 R 版本,然后启动 R 并输入以下命令安装最新的 Bioconductor:

## try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
biocLite()

通过上述命令,将安装该版本发行版中的所有核心包。要安装 Bioconductor 仓库中可用的特定包(如 IMMAN 和 GenomicRanges),可按各包落地页的说明执行,例如:

source("https://bioconductor.org/biocLite.R")
biocLite(c( "Biostrings", "GenomicRanges", "IMMAN"))

BiocInstaller 包中的 biocLite() 函数用于安装包,而非 R 常用的 install.packages()。原因在于 Bioconductor 使用独立于 CRAN 的仓库,且其发布节奏与 R 不同,这可能导致 R 与 Bioconductor 的发行时间不匹配。由 install.packages() 识别的版本往往并非最新的“release”。

使用示例

完成初始安装后,本部分通过一些简单示例让您熟悉 Bioconductor 的工作环境。这里汇集了该项目的不同组件,如高通量测序、统计分析与理解的两个实例。

首先,为熟悉高通量数据,您将接触 DNA 序列。为此,您需要使用 Biostring,这是处理大型生物序列的核心包。您可以使用 DNAStringSet() 函数来存储一组 DNAString(基于 IUPAC 扩展遗传字母表的 DNA 字符串)。它可将输入转换为 DNA 碱基类型的 XStringSet 对象。字母的编码方式有利于加速搜索算法。

library(Biostrings)

## Loading required package: BiocGenerics

## Loading required package: parallel

##
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:parallel':
##
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB

## The following objects are masked from 'package:stats':
##
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
##
##     anyDuplicated, append, as.data.frame, basename, cbind,
##     colMeans, colnames, colSums, dirname, do.call, duplicated,
##     eval, evalq, Filter, Find, get, grep, grepl, intersect,
##     is.unsorted, lapply, lengths, Map, mapply, match, mget, order,
##     paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind,
##     Reduce, rowMeans, rownames, rowSums, sapply, setdiff, sort,
##     table, tapply, union, unique, unsplit, which, which.max,
##     which.min

## Loading required package: S4Vectors

## Loading required package: stats4

##
## Attaching package: 'S4Vectors'

## The following object is masked from 'package:base':
##
##     expand.grid

## Loading required package: IRanges

##
## Attaching package: 'IRanges'

## The following object is masked from 'package:grDevices':
##
##     windows

## Loading required package: XVector

##
## Attaching package: 'Biostrings'

## The following object is masked from 'package:base':
##
##     strsplit

dnaSequence <- DNAStringSet( c("AAACTG", "CCCAACCA") )
dnaSequence

##   A DNAStringSet instance of length 2
##     width seq
## [1]     6 AAACTG
## [2]     8 CCCAACCA

接着,我们可以使用 complement() 函数来计算 DNA 序列的互补序列。

complement(dnaSequence)

##   A DNAStringSet instance of length 2
##     width seq
## [1]     6 TTTGAC
## [2]     8 GGGTTGGT

学习使用诸如 DNAStringSet 与 Biostrings 的包,有助于掌握其他包的用法。

在另一个示例中,我们将使用 GenomicRanges 包,它用于在 Bioconductor 项目中表示基因组位置。它通过提供三个类(GRanges、GPos 和 GRangesList)为基因组分析奠定基础,分别用于表示基因组范围、基因组位置和基因组范围的集合。

GRanges 类用于表示一组基因组范围,每个范围在基因组上都有唯一的起止位置。它可存储连续结合位点、转录本与外显子等基因组特征的位置。我们可以像下面这样通过 GRanges 构造函数访问这些对象:

library(GenomicRanges)

## Loading required package: GenomeInfoDb

grangeObj <-
  GRanges(seqnames =
            Rle(c("chr1", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)),
          ranges =
            IRanges(1:10, end = 7:16, names = head(letters, 10)),
          strand =
            Rle(strand(c("-", "+", "*", "+", "-")),
                c(1, 2, 2, 3, 2)),
          score = 1:10,
          GC = seq(1, 0, length=10))

grangeObj

## GRanges object with 10 ranges and 2 metadata columns:
##     seqnames    ranges strand |     score                GC
##        <Rle> <IRanges>  <Rle> | <integer>         <numeric>
##   a     chr1       1-7      - |         1                 1
##   b     chr2       2-8      + |         2 0.888888888888889
##   c     chr2       3-9      + |         3 0.777777777777778
##   d     chr2      4-10      * |         4 0.666666666666667
##   e     chr1      5-11      * |         5 0.555555555555556
##   f     chr1      6-12      + |         6 0.444444444444444
##   g     chr3      7-13      + |         7 0.333333333333333
##   h     chr3      8-14      + |         8 0.222222222222222
##   i     chr3      9-15      - |         9 0.111111111111111
##   j     chr3     10-16      - |        10                 0
##   -------
##   seqinfo: 3 sequences from an unspecified genome; no seqlengths

结果是一个包含 10 个基因组范围的 GRanges 对象。如上所示,输出被“|”分隔为左右两部分。左侧为基因组坐标(seqnames、ranges 与 strand),右侧为元数据列(注释)。在此示例中,元数据包括 score 与 GC 信息。当然,几乎任何信息都可以存储在 GRanges 对象的元数据部分。

seqnames()、ranges() 与 strand() 等辅助函数可用于提取 GRanges 对象中基因组坐标的各组成部分。例如,

seqnames(grangeObj)

## factor-Rle of length 10 with 4 runs
##   Lengths:    1    3    2    4
##   Values : chr1 chr2 chr1 chr3
## Levels(3): chr1 chr2 chr3

ranges(grangeObj)

## IRanges object with 10 ranges and 0 metadata columns:
##         start       end     width
##     <integer> <integer> <integer>
##   a         1         7         7
##   b         2         8         7
##   c         3         9         7
##   d         4        10         7
##   e         5        11         7
##   f         6        12         7
##   g         7        13         7
##   h         8        14         7
##   i         9        15         7
##   j        10        16         7

strand(grangeObj)

## factor-Rle of length 10 with 5 runs
##   Lengths: 1 2 2 3 2
##   Values : - + * + -
## Levels(3): + - *

最后一个示例展示了该项目统计分析部分的一个简单用法,即基因本体(GO)分类。GO 是一组用于描述基因功能及其相互关系的概念(分子功能、细胞组分与生物过程)。更具体地,GO 将基因及其产物按层级组织成一个称为本体的图。为获得可用示例,我们可以使用三个 Bioconductor 包:clusterProfiler、DOSE 和 org.Hs.eg.db。

library("clusterProfiler")

##

## clusterProfiler v3.8.1  For help: https://guangchuangyu.github.io/software/clusterProfiler
##
## If you use clusterProfiler in published research, please cite:
## Guangchuang Yu, Li-Gen Wang, Yanyan Han, Qing-Yu He. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology. 2012, 16(5):284-287.

library("DOSE")

## DOSE v3.6.1  For help: https://guangchuangyu.github.io/DOSE
##
## If you use DOSE in published research, please cite:
## Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an R/Bioconductor package for Disease Ontology Semantic and Enrichment analysis. Bioinformatics 2015, 31(4):608-609

library("org.Hs.eg.db")

## Loading required package: AnnotationDbi

## Loading required package: Biobase

## Welcome to Bioconductor
##
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

##

clusterProfiler 提供了一组用于分析与可视化基因及基因簇功能谱(如 GO)的特定方法。借此,您可以按相似性对不同基因进行聚类。作为示例,DOSE 包中的 geneList 数据集是不错的选择,它包含基因 ID 的信息量;另一方面,org.Hs.eg.db 包则为人类提供了全基因组注释,即富含与基因组位置相关信息的遗传物质序列表示。

data(geneList, package="DOSE")

gene <- names(geneList)[abs(geneList) > 2]

enrichGO() 函数对给定的基因向量执行 GO 富集分析。富集分析是一种方法,用于根据功能特征将一组基因分配到预定义类别中,并识别显著富集的类别。富集结果可能包含非常一般的术语。使用该函数时,您需要填写名为“OrgDb”的参数;为此,需要包含人类全基因组注释的 org.Hs.eg.db 包。该包由遗传物质的序列表示构成,并富含与基因组位置相关的信息。

ego <- enrichGO(gene          = gene,
                universe      = names(geneList),
                OrgDb         = org.Hs.eg.db,
                ont           = "CC",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.01,
                qvalueCutoff  = 0.05,
        readable      = TRUE)
head(ego)

##                    ID                              Description GeneRatio
## GO:0005819 GO:0005819                                  spindle    25/201
## GO:0000779 GO:0000779 condensed chromosome, centromeric region    15/201
## GO:0072686 GO:0072686                          mitotic spindle    14/201
## GO:0000775 GO:0000775           chromosome, centromeric region    18/201
## GO:0000776 GO:0000776                              kinetochore    15/201
## GO:0000793 GO:0000793                     condensed chromosome    18/201
##              BgRatio       pvalue     p.adjust       qvalue
## GO:0005819 256/11817 1.353146e-12 3.897060e-10 3.546666e-10
## GO:0000779  91/11817 3.013479e-11 3.393318e-09 3.088217e-09
## GO:0072686  77/11817 3.534706e-11 3.393318e-09 3.088217e-09
## GO:0000775 156/11817 1.429309e-10 1.029103e-08 9.365738e-09
## GO:0000776 106/11817 2.849245e-10 1.406873e-08 1.280378e-08
## GO:0000793 163/11817 2.970249e-10 1.406873e-08 1.280378e-08
##                                                                                                                                                         geneID
## GO:0005819 CDCA8/CDC20/KIF23/CENPE/ASPM/DLGAP5/SKA1/NUSAP1/TPX2/TACC3/NEK2/CDK1/MAD2L1/KIF18A/BIRC5/KIF11/TTK/AURKB/PRC1/KIFC1/KIF18B/KIF20A/AURKA/CCNB1/KIF4A
## GO:0000779                                                            CENPE/NDC80/HJURP/SKA1/NEK2/CENPM/CENPN/ERCC6L/MAD2L1/CDT1/BIRC5/NCAPG/AURKB/AURKA/CCNB1
## GO:0072686                                                                KIF23/CENPE/ASPM/NUSAP1/TPX2/TACC3/CDK1/MAD2L1/KIF18A/KIF11/AURKB/KIFC1/KIF18B/AURKA
## GO:0000775                                           CDCA8/CENPE/NDC80/HJURP/SKA1/NEK2/CENPM/CENPN/ERCC6L/MAD2L1/KIF18A/CDT1/BIRC5/TTK/NCAPG/AURKB/AURKA/CCNB1
## GO:0000776                                                             CENPE/NDC80/HJURP/SKA1/NEK2/CENPM/CENPN/ERCC6L/MAD2L1/KIF18A/CDT1/BIRC5/TTK/AURKB/CCNB1
## GO:0000793                                          CENPE/NDC80/TOP2A/NCAPH/HJURP/SKA1/NEK2/CENPM/CENPN/ERCC6L/MAD2L1/CDT1/BIRC5/NCAPG/AURKB/CHEK1/AURKA/CCNB1
##            Count
## GO:0005819    25
## GO:0000779    15
## GO:0072686    14
## GO:0000775    18
## GO:0000776    15
## GO:0000793    18

为可视化富集图谱,可以使用 emapplot() 函数。该函数将基因集展示为网络。简单来说,彼此重叠的基因集会聚成簇。

emapplot(ego)
enrichment map

获取帮助

在使用 Bioconductor 各类包的过程中,如需帮助,您可以使用以下命令:

class(dnaSequence) # Indicate the class of the input object

## [1] "DNAStringSet"
## attr(,"package")
## [1] "Biostrings"

methods(class = "DNAStringSet") # List all available methods for a S3 and S4 generic function or class

##   [1] !                               !=                             
##   [3] $                               $<-                            
##   [5] %in%                            [                              
##   [7] [[                              [[<-                           
##   [9] [<-                             <                              
##  [11] <=                              ==                             
##  [13] >                               >=                             
##  [15] aggregate                       alphabetFrequency              
##  [17] anyNA                           append                         
##  [19] as.character                    as.complex                     
##  [21] as.data.frame                   as.env                         
##  [23] as.factor                       as.integer                     
##  [25] as.list                         as.logical                     
##  [27] as.matrix                       as.numeric                     
##  [29] as.raw                          as.vector                      
##  [31] by                              c                              
##  [33] cbind                           chartr                         
##  [35] coerce                          compact                        
##  [37] compareStrings                  complement                     
##  [39] concatenateObjects              consensusMatrix                
##  [41] consensusString                 countOverlaps                  
##  [43] countPattern                    countPDict                     
##  [45] dinucleotideFrequencyTest       do.call                        
##  [47] droplevels                      duplicated                     
##  [49] elementMetadata                 elementMetadata<-              
##  [51] elementNROWS                    elementType                    
##  [53] eval                            expand                         
##  [55] expand.grid                     extractAt                      
##  [57] extractROWS                     Filter                         
##  [59] findOverlaps                    getListElement                 
##  [61] hasOnlyBaseLetters              head                           
##  [63] ifelse2                         intersect                      
##  [65] is.na                           is.unsorted                    
##  [67] isEmpty                         isMatchingEndingAt             
##  [69] isMatchingStartingAt            lapply                         
##  [71] length                          lengths                        
##  [73] letterFrequency                 match                          
##  [75] matchPattern                    matchPDict                     
##  [77] mcols                           mcols<-                        
##  [79] merge                           metadata                       
##  [81] metadata<-                      mstack                         
##  [83] names                           names<-                        
##  [85] nchar                           neditEndingAt                  
##  [87] neditStartingAt                 NROW                           
##  [89] nucleotideFrequencyAt           oligonucleotideFrequency       
##  [91] order                           overlapsAny                    
##  [93] PairwiseAlignments              PairwiseAlignmentsSingleSubject
##  [95] parallelSlotNames               parallelVectorNames            
##  [97] pcompare                        pcompareRecursively            
##  [99] PDict                           PWM                            
## [101] rank                            rbind                          
## [103] Reduce                          relist                         
## [105] relistToClass                   rename                         
## [107] rep                             rep.int                        
## [109] replaceAt                       replaceLetterAt                
## [111] replaceROWS                     rev                            
## [113] revElements                     reverse                        
## [115] reverseComplement               ROWNAMES                       
## [117] sapply                          seqinfo                        
## [119] seqinfo<-                       seqlevelsInUse                 
## [121] seqtype                         seqtype<-                      
## [123] setdiff                         setequal                       
## [125] setListElement                  shiftApply                     
## [127] show                            showAsCell                     
## [129] sort                            split                          
## [131] split<-                         stack                          
## [133] stringDist                      strsplit                       
## [135] subseq                          subseq<-                       
## [137] subset                          subsetByOverlaps               
## [139] table                           tail                           
## [141] tapply                          threebands                     
## [143] toString                        transform                      
## [145] translate                       trimLRPatterns                 
## [147] twoWayAlphabetFrequency         union                          
## [149] unique                          uniqueLetters                  
## [151] unlist                          unsplit                        
## [153] unstrsplit                      updateObject                   
## [155] values                          values<-                       
## [157] vcountPattern                   vcountPDict                    
## [159] vmatchPattern                   vwhichPDict                    
## [161] which.isMatchingEndingAt        which.isMatchingStartingAt     
## [163] whichPDict                      width                          
## [165] window                          window<-                       
## [167] windows                         with                           
## [169] within                          xtabs                          
## [171] xvcopy                          zipdown                        
## see '?methods' for accessing help and source code

?complement # Show the help page of the function

## starting httpd help server ... done

browseVignettes(package="Biostrings") # Refers to all available vignettes of the input package

此外,Bioconductor 站点上还提供了大量实用的帮助页面与工作流程:

总结

阅读本教程后,您已了解 Bioconductor 项目及其在生物研究中的基础性作用。同时也介绍了 Bioconductor 包的类型及部分重要包。随后简述了安装与使用该项目的简短工作流程。您学习了如何使用 Biostring 与 GenomicRanges 包,从而熟悉了 DNA 序列与基因组范围。在另一个示例中,您尝试对一组人类基因进行富集分析,因而了解了富集分析与 GO 的相关概念。最后也给出了在使用 Bioconductor 过程中获取帮助的多种途径。若您想进一步了解该项目,欢迎访问 Bioconductor_courses。

如果您想进一步学习 R,请参加 DataCamp 的 Experimental Design in R 课程。

主题