作者:马丁·摩根(mtmorgan@fredhutch.org）
时间:2015年9月7日
回车间大纲

本文件中的材料要求R3.2版及Bioconductor版本3.1

stopifnot(getRversion() >= '3.2' && getRversion() < '3.3'， BiocInstaller::biocVersion() >= "3.1")

总体工作流程

1. 实验设计
   • 保持简单!
   • 复制!
   • 避免或跟踪批效应
2. 湿实验室制备
3. 高通量测序
   • 输出:读取的FASTQ文件及其质量分数
4. 对齐
   • 许多不同的对齐器，一些专门用于不同的目的
   • 输出:对齐读取的BAM文件
5. 总结
   • 例如,数感兴趣的重叠区域(如基因)
6. 统计分析
7. 理解

笔记

• 新创组装在…范围之外Bioconductor

两个简单的闪亮的应用程序

• BAMSpector -跨BAM(对齐读取)文件显示基因模型和底层支持

  App <- system。file(package="BiocAsiaPacific2015"， "BAMSpector") shiny::runApp(app)

• MAPlotExplorer -总结两组差异表达，包括个别基因的钻取。基于CSAMA 2015实验室的Andrzej Oles。

  App <- system。file(package="BiocAsiaPacific2015"， "MAPlotExplorer") shiny::runApp(app)

如何Bioconductor帮助

注释

• 基因(例如，外显子内转录)模型
• 标识符映射(例如，'HNRPC'到基因模型中使用的Entrez标识符)

标准(大)文件输入和操作，例如，BAM文件的对齐读取

统计分析

• 微分表达式
• ChIPSeq
• 变体
• 流cytometery
• 蛋白质组学
• .．.

关键资源

• 网站:http://biocondcutor.org

  • biocViews
  • 包的“着陆页”，例如:刨边机
  • 帮助
  • 开发人员资源:新包提交，匈牙利瑞士比分开发人员邮件列表等。

• 支持:https://support.bioconductor.org

R复习

统计计算和图形的语言和环境

• 全功能的编程语言
• 互动和解释-方便和宽容
• 连贯的、广泛的文档
• 统计。因素()，NA
• 可扩展- CRAN, Bioconductor, github，…

向量，类，对象

• 非常高效。矢量化原子向量的计算逻辑，整数，数字，复杂的，字符，字节
• 原子向量是更复杂的构造块对象
  • 矩阵-带有'dim'属性的原子向量
  • data.frame-长度相等的原子向量列表
• 正式的类表示向量的复杂组合，例如的返回值lm (),低于

函数，泛型，方法

• 函数将输入转换为输出，可能会有副作用，例如，rnorm (1000)
  • 参数先按名称匹配，然后按位置匹配
  • 函数可以定义(一些)参数为默认值
• 通用的函数分派到特定的方法基于参数的类别，例如，print ()．
• 方法是实现特定泛型的函数，例如，print.factor；方法被调用间接，通过通用。
• 许多但不是所有能够操作特定类的函数都是方法，例如:abline ()下面使用的是一个普通的旧函数。

编程

迭代:

• 拉普兰人()

  args(拉普兰人)

  ##函数(X, FUN，…

  • 意义:对于一个向量X(通常是一个列表())，应用一个函数有趣的到每个向量元素，返回结果为l坚持。.．.是附加参数有趣的．
  • 有趣的可以是内置的，也可以是用户定义的函数

  LST <- list(a=1:2, b=2:4) lapply(LST, log) # 'base'参数default;自然对数

  ## $ [1] 0.0000000 0.6931472 ## ## $ b# # [1] 0.6931472 1.0986123 1.3862944

  Lapply (lst, log, 10) # '10'是'log()'的第二个参数，即log以10为底

  ## $ [1] 0.00000 0.30103 ## ## $b ## [1] 0.3010300 0.4771213 0.6020600

  • 酸式焦磷酸钠()——就像拉普兰人()，但如果可能的话，将结果简化为向量、矩阵或数组。
  • vapp ()——就像酸式焦磷酸钠()，但要求返回类型为有趣的指定;这可能更安全——当结果为意外类型时出现错误。

• 宾州()(也Map ()）

  args (mapp)

  ##函数(FUN，…，MoreArgs = NULL, SIMPLIFY = TRUE, USE.NAMES = TRUE) ## NULL

  • .．.是一个或多个向量，循环为相同长度。有趣的一个函数的参数和它的组成部分一样多吗.．.．宾州返回应用的结果有趣的向量的元素.．.．

  mapply(seq, 1:3, 4:6, SIMPLIFY=FALSE) # seq(1,4);seq(2、5);seq(3、6)

  # # # # [1] ([1]) 1 2 3 4 # # # # # # [1] ([2]) 2 3 4 5 # # # # # # [1] [[3]] 3 4 5 6

  • 应用()

  args(应用)

  ##函数(X, MARGIN, FUN，…

  • 对于矩阵或数组X,应用有趣的对每一个保证金(尺寸,例如,利润= 1方法应用有趣的对于每一行，利润= 2方法应用有趣的到每一列)

• 传统的迭代编程结构重复{}，对于(){}

  • 几乎总是更容易出错，效率更低，更难理解拉普兰人()！

有条件的

if (test) {## code if test == TRUE} else {## code if test == FALSE}

函数(参见下表中的一些收藏)

• 很容易定义自己的函数

fun <- function(x) {length(unique(x))} ##长度为5的列表，每个包含一个字母的样本(替换)let <- replication (5, sample(letters, 50, TRUE)， simplify=FALSE) sapply(lets, fun)

## [1] 25 23 21 22 20

反省和帮助

自省

• 一般性质，例如:类()，str ()
• 类特定的属性，例如，昏暗的()

帮助

• ?“打印”:帮助一般打印
• " ? print.data.frame:帮助打印方法类的对象data.frame。
• 帮助(包=“GenomeInfoDb”)
• browseVignettes(“GenomicRanges”)
• 方法(“阴谋”)
• 方法(class = " lm”)

例子

R向量，向量化运算，data.frame ()、公式、函数、对象、类和方法发现(自省)。

X <- rnorm(1000) #原子向量y <- X + rnorm(1000, sd=.5) df <- data.frame(X = X, y=y) #类'data.frame'的对象plot(y ~ X, df) #泛型plot，方法plot。公式fit <- lm(y ~x, df) #类'lm'方法的对象(class=class(fit)) #自省

##[1]添加一个别名anova case.names强制限制##[7]厨师。距离偏差dfbeta dfbetas drop1 dummy。coef ## [13] effects extractAIC family formula hatvalues influence ## [19] initialize kappa labels logLik model.frame模型。矩阵## [25] nobs plot predict print proj qr ## [31] residuals rstandard rstudent show simulate slotsFromS3 ## [37] summary variable.names vcov ## see '?methods' for accessing help and source code

方差分析(适合)

##方差分析表## ##响应:y ## Df Sum Sq Mean Sq F value Pr(>F) ## x 1 963.82 963.82 3649.1 < 2.2e-16 *** #残差998 263.60 0.26 ##—##显著值。代码:0 '***' 0.001 '**' 0.01 '*' 0.05 '。' 0.1 ' ' 1

Plot (y ~ x, df) #方法(Plot);情节吗?。公式abline(fit, col="red"， lwd=3, lty=2) #一个函数，而不是generic.method

编程示例-将1000个符号分组为GO标识符

##示例数据fl <- file.choose() ## symbol .csv .csv

symgo <- read.csv(fl, row.names=1, stringsAsFactors=FALSE)

##符号走证据本体## 1 ppiap28 < na > < na > < na > ## 2 pptlah < na > < na > < na > ## 3 hist1h2bc go:0000786 nas cc ## 4 hist1h2bc go:0000788 iba cc ## 5 hist1h2bc go:0002227 IDA bp ## 6 hist1h2bc go:0003677 iba mf

暗(symgo)

## [1] 5041 4

长度(独特(symgo $符号))

## [1] 1000

## split-sapply go2sym <- split(symgo$SYMBOL, symgo$GO) len1 <- sapply(go2sym, length) #比较lapply, vapply ##内置函数的通用操作len2 <-长度(go2sym)相同(len1, len2)

##[1]真

##更智能的内置函数，例如，省略NAs len3 <- aggregate(SYMBOL ~ GO, symgo, length) head(len3)

## go符号## 1 go:0000049 3 ## 2 go:0000050 2 ## 3 go:0000060 1 ## 4 go:0000077 1 ## 5 go:0000086 3 ## 6 go:0000118 1

##使用aggregate() head(aggregate(GO ~ SYMBOL, symgo, length))

##符号go ## 1 abcd4 15 ## 2 abcg2 22 ## 3 ace 57 ## 4 adamtsl2 6 ## 5 aldh1l2 11 ## 6 alox5 19

head(aggregate(SYMBOL ~ GO, symgo, c))

## go符号## 1 go:0000049 yars2, yars2, eef1a1 ## 2 go:0000050 asl, asl ## 3 go:0000060 oprd1 ## 4 go:0000077 pea15 ## 5 go:0000086 tubb4a, cenpf, clasp1 ## 6 go:0000118 cir1

uidfun <- function(x) {unique(tolower(x))} head(aggregate(SYMBOL ~ GO, symgo, uidfun))

## GO符号## 1 GO:0000049 yars2, eef1a1 ## 2 GO:0000050 asl ## 3 GO:0000060 oprd1 ## 4 GO:0000077 pea15 ## 5 GO:0000086 tubb4a, cenpf, clasp1 ## 6 GO:0000118 cir1

##作为“匿名”函数头(aggregate(SYMBOL ~ GO, symgo, function(x) {unique(tolower(x))})))

## GO符号## 1 GO:0000049 yars2, eef1a1 ## 2 GO:0000050 asl ## 3 GO:0000060 oprd1 ## 4 GO:0000077 pea15 ## 5 GO:0000086 tubb4a, cenpf, clasp1 ## 6 GO:0000118 cir1

案例研究:所有表型数据

这些案例研究可以作为复习R数据的输入和操作。

输入一个包含ALL(急性淋巴细胞白血病)患者信息的文件

fname <- file.choose() ## " all表型数据。Tsv " stopifnot(file.exists(fname)) pdata <- read.delim(fname)

查看帮助页面read.delim ?对于输入选项，并探索您所创建的对象的基本属性，例如…

类(pdata)

## [1] "data.frame"

colnames (pdata)

# #[1]“id”“诊断”“性”“年龄”“转基因”# #[6]“缓解”“CR”“日期。cr”“t.4.11。”“t.9.22。”## [11] "cyto。正常" "citog" "mol.biol" "融合。蛋白质"mdr" ## [16] "kinet" "ccr" "复发" "移植" "f.u" ##[21] "日期。最后。见"

暗(pdata)

## [1] 127 21

头(pdata)

id诊断性别年龄BT缓解CR日期。cr t.4.11。t.9.22。阶段。正常citog # # 1 1005 5/21/1997 53 B2 CR CR 8/6/1997假真的假t (9; 22) # # 2 1010 3/29/2000 19 B2 CR CR 6/27/2000假假假简单的alt。# # 3 3002 6/24/1998 F 52 B4 CR CR 8/17/1998 NA NA NA < NA > # # 4 4006 7/17/1997 38 B1 CR CR 9/8/1997真的假假的t (4, 11) # # 5 4007 7/22/1997 57 B2 CR CR 9/17/1997假假假del q (6) # # 6 4008 7/30/1997 17 B1 CR CR 9/27/1997假假假复杂的alt。# # mol.biol融合。蛋白质mdr kinet ccr复发移植f.u日期。最后见## 1 BCR/ABL p210阴性异倍体假假真BMT /死亡CR  ## 2 NEG  POS异倍体假真假REL 2000年8月28日## 3 BCR/ABL p190阴性异倍体假真假REL 1999年10月15日## 4 ALL1/AF4 阴性异倍体假真假REL 1998年1月23日## 5 NEG 阴性异倍体假真假REL 1997年4月11日## 6 NEG 阴性超倍体。假真假rel 1997年12月15日

总结(pdata性美元)

## F M NA的## 42 83 2

总结(pdata cyto.normal美元)

##模式FALSE TRUE NA的##逻辑69 24 34

提醒自己有各种方法来子集和访问data.frame的列

pdata [1:5, 3:4)

性别年龄## 1岁53 ## 2岁19 ## 3岁52 ## 4岁38 ## 5岁57

pdata [1:5]

id诊断性别年龄BT缓解CR日期。cr t.4.11。t.9.22。阶段。正常citog mol.biol # # 1 1005 5/21/1997 53 B2 CR CR 8/6/1997假真的假t (9; 22) BCR / ABL # # 2 1010 3/29/2000 19 B2 CR CR 6/27/2000假假假简单的alt。NEG # # 3 3002 6/24/1998 F 52 B4 CR CR 8/17/1998 NA NA NA < NA > BCR / ABL # # 4 4006 7/17/1997 38 B1 CR CR 9/8/1997真的假假的t (4, 11) ALL1 / AF4 # 57 # 5 4007 7/22/1997 B2 CR CR 9/17/1997假假假德尔(6问)底片# #融合。蛋白质mdr kinet ccr复发移植f.u date. lastseen ## 1 p210阴性dyploid假假真BMT / DEATH IN CR  ## 2  POS dyploid假真假REL 2000年8月28日## 3 p190阴性dyploid假真假REL 1999年10月15日## 4 阴性dyploid假真假REL 1998年1月23日## 5 阴性dyploid假真假REL 1997年11月4日

头(pdata [3:5])

##性别年龄BT ## 1 M 53 B2 ## 2 M 19 B2 ## 3 F 52 B4 ## 4 M 38 B1 ## 5 M 57 B2 ## 6 M 17 B1

Tail (pdata[， 3:5]， 3)

性别年龄BT ## 125 M 19 T2 ## 126 M 30 T3 ## 127 M 29 T2

头(pdata时代美元)

## [1] 53 19 52 38 57 17

头(pdata性美元)

## [1] M M F M M M M ##级别:F M

Head (pdata[pdata$年龄> 21，])

id诊断性别年龄BT缓解CR日期。cr t.4.11。t.9.22。阶段。正常citog # # 1 1005 5/21/1997 53 B2 CR CR 8/6/1997假真的假t (9; 22) # # 3 3002 6/24/1998 F 52 B4 CR CR 8/17/1998 NA NA NA < NA > # # 4 4006 7/17/1997 38 B1 CR CR 9/8/1997真的假假的t (4, 11) # # 5 4007 7/22/1997 57 B2 CR CR 9/17/1997假假假德尔(6问)# # 10 8001 1/15/1997 40 B2 CR CR 3/26/1997假假假德尔(p15) # # 11 8011 8/21/1998 33 B3 CR CR 10/8/1998假假假德尔(p15 / p16) # # mol.biol融合。mdr蛋白其ccr复发移植f.u date.last.seen # # 1 BCR / ABL p210 NEG dyploid假假真BMT死亡/ CR < NA > # # 3 BCR / ABL p190 NEG dyploid假真的假REL 10/15/1999 # # 4 ALL1 / AF4 < NA > NEG dyploid假真的假REL 1/23/1998 # # 5 NEG < NA > NEG dyploid假真的假REL 11/4/1997 # # 10 BCR / ABL p190 NEG < NA >假真的假REL 7/11/1997 # # 11 BCR / ABL p190 / p210 NEG dyploid假假真BMT死亡/ CR < NA >

从下面看，数据集中有17位40岁以上的女性，但当细分pdata为了包含这些个体，我们选择了19行。为什么?我们能做些什么来纠正这一点呢?

idx <- pdata$sex == "F" & pdata$age > 40表(idx)

## idx ##假的真实的## 108 17

暗(pdata [idx])

## [1] 19 21

使用mol.biol列，以将数据子集化，只包含“BCR/ABL”或“NEG”的个体，例如:

Bcrabl <- pdata[pdata$mol.]生物%in% c(“BCR/ABL”，“NEG”)，]

的mol.biol列是一个因子，即使在子集设置之后，它仍然保留所有级别。如何去掉未使用的因子级别?

bcrabl美元摩尔。Biol <- factor(bcrabl$mol.biol)

的英国电信列是描述B细胞和t细胞亚型的因子

水平(bcrabl BT美元)

# #[1]“B”“B1”“B2”“B3”“B4”“T”“T1”“T2”“T3”“T4”

如何将B1 B2…分解为单一类型B, T1 T2…也是如此，所以只有两种子类型B和T

表(bcrabl BT美元)

## ## b b1 b2 b3 b4 t t1 t2 t3 t4 ## 4 9 35 22 9 4 1 15 9 2

level (bcrabl$BT) <- substring(level (bcrabl$BT)， 1,1) table(bcrabl$BT)

## ## b t ## 79 31

使用xtabs ()(交叉制表)计算BCR/ABL组和NEG组中B-和t细胞类型的样本数量

xtabs(~ BT + mol.biol, bcrabl)

## mol.biol ## BT BCR/ABL NEG ## B 37 42 ## T 0 31

使用总()计算BCR/ABL和NEG治疗组男性和女性的平均年龄。

聚合(年龄~ mol.biol +性别，字母，平均值)

##摩尔生物学性别年龄## 1 BCR/ABL F 39.93750 ## 2阴性F 30.42105 ## 3 BCR/ABL M 40.50000 ## 4阴性M 27.21154

使用t.test ()比较BCR/ABL组与NEG组个体的年龄;使用以下方法将结果可视化箱线图()．在这两种情况下，都使用公式接口。咨询帮助页t.test ?并重新做测试，假设两组的年龄方差相同。测试输出的哪些部分会发生变化?

T.test(年龄~ mol.biol, bcrabl)

## ## Welch Two样本t-检验## ##数据:年龄由摩尔。biol ## t = 4.8172, df = 68.529, p-value = 8.401e-06 ##替代假设:均值的真实差异不等于0 ## 95%置信区间:## 7.13507 17.22408 ##样本估计:## BCR组均值/ ## ABL组均值NEG ## 40.25000 28.07042

箱线图(年龄~ mol.biol, bcrabl)