1.实验设计

保持简单

• 经典实验设计
• 时间序列
• 在可能的情况下，不要丢失值
• 预期的分析必须是可行的——可用的样本和感兴趣的假设是否可以结合起来，形成一个可检验的统计假设?

复制

• 复制的程度决定了生物学问题的细微差别。
• 无复制(每个处理1个样本):定性描述，统计选项有限。
• 每个处理3-5个重复:用细胞系或其他定义良好的实体设计实验操作;两组之间的平均表达变化为2倍(?)
• 每次治疗10-50个重复:群体研究，例如癌细胞系。
• 1000个重复:前瞻性研究，例如，SNP发现
• 一个资源:RNASeqPower

避免将实验因素与其他因素混淆

• 常见问题:同一处理的样品都在同一流池上;处理1的样品先处理，处理2的样品再处理，以此类推。

记录co-variates

注意批处理的影响

• 已知的

  • 表型协变量，如年龄、性别
  • 实验协变量，如实验室或处理日期
  • 至少近似地纳入线性模型

• 未知的

  • 或者只是意想不到/未被发现
  • 描述使用，例如，股东价值分析．

• 替代变量分析

  • Leek等人，2010，自然评论遗传学11733 - 739，韭菜和故事PLoS Genet 3(9):e161．
  • 科学发现:普遍批量效应
  • 统计洞察:替代变量分析:识别和构建替代变量;删除已知的批处理效果
  • 好处:减少依赖，稳定错误率估计，并提高可重复性
  • 战斗软件/股东价值分析Bioconductor包

HapMap样品来自一个工厂，按加工日期订购。

2.湿实验室

的混杂因素

• 记录或避免

你的文物特定的协议

• 序列的污染物
• 富集偏向，例如，不均匀的转录本表现。
• PCR伪影-适配器污染物，序列特异性扩增偏置，…

3.测序

变化轴

• 单端与配对端
• 长度:50 - 200元
• 每个样本的读取数

特定于应用程序,例如,

• ChIP-seq:短的单端读取通常就足够了
• RNA-seq，已知基因:单端或成对读取
• RNA-seq，转录本或新变体:配对末端读取
• 拷贝数:单端或成对端读取
• 结构变体:成对的末端读取
• 变体:深度通过更长的、成对的末端读取
• 微生物组:长配对末端读取(末端重叠)

4.对齐

定位策略

• 新创
• 无参考基因组;可观的排序和计算资源
• 基因组
• 已建立的参考基因组
• Splice-aware对准器
• 新型转录本发现
• 转录组
• 建立参考基因组;可靠的基因模型
• 简单的调整器
• 已知基因/转录本表达

感知拼接的对齐器(和Bioconductor包装器)

• Bowtie2（Rbowtie）
• 明星（doi）
• subread（Rsubread）
• 系统评价(Engstrom et al.， 2013，doi）

(5。领结/大礼帽/袖扣/袖扣等)

• 大礼帽使用Bowtie2来执行基本的单端和配对端对齐，然后使用算法将难以对齐的读取数据放置在对齐良好的读取数据附近。
• 袖扣（doi)大礼帽输出和估计现有的和新的转录丰度。袖扣是如何工作的
• Cuffdiff评估实验组之间估计丰度的统计显著性
• RSEM包括从头组装和量化

5.简化为“计数表”

• 使用已知的基因模型来计数感兴趣的重叠区域/基因模型的对齐reads
• 基因模型可以是公开的(如UCSC, NCBI, ENSEMBL)或特别的(人造石铺地面文件)
• GenomicAlignments: summarizeOverlaps ()
• Rsubread: featureCount ()
• HTSeq，htseq-count

6.分析

独特的统计方面

• 大数据，少样本
• 跨样本的每个基因的比较;单变量措施
• 对每个基因进行分析相同实验设计，下相同零假设

摘要

• 计数本身，而不是一个总结(RPKM, FRPKM，…)，是相关的分析
• 对于一个给定的基因，数量越大意味着信息越多;RPKM等，将所有估计视为同等信息。
• 比较是在每一个感兴趣的区域;所有样本都有相同的兴趣区域，因此模库大小不需要校正，例如，基因长度或可映射性。

归一化

• 库的大小(每个样本的总读取计数)因一些无趣的原因而不同;我们需要在统计分析中考虑到图书馆大小的差异。
• 每个样本的总计数读数为不一个很好的库大小估计。它不一定会受到数量大的区域的影响，而且可能会在基因之间引入偏差和相关性。相反，使用考虑计数分布的倾斜的库大小的稳健度量(最简单:修剪的几何平均值;更高级/更适合在实验室遇到)。
• 库大小(已计数的读取总数)在不同的样本之间是不同的，并且应该包括在内作为一个统计偏差在差分表达式的分析中，而不是在分析的早期“除以”库的大小。

合适的误差模型

• 计数数据为不正态分布或泊松分布，而是负二项分布。
• 结合泊松(“射击”噪声，即样本内技术和读取计数的采样变化)与生物样品之间的变化的结果。
• 负二项式模型需要估计一个额外的参数(“离散度”)，在小样本中估计得很差。
• 基本策略是通过从具有相似表达值的基因中获得更可靠的局部估计来调节每个基因的估计(下面将提供更多关于借用信息的信息)。

预滤器

• 简单地说，统计检验(例如t检验)可以应用于计数表的每一行。然而，我们有相对较少的样本(10个)和非常多的比较(10,000个)，所以一个幼稚的方法可能是非常不足的，导致一个非常高的错误发现率
• 一个简单的方法是通过删除不可能产生统计显著性的区域来进行更少的检验，而不管考虑的假设是什么。
• 例如:在所有样本中计数为0的区域不可能是显著的假设再梯度，因此排除在进一步分析之外。
• 基本方法:' K / A '风格的过滤器-要求在至少K个样本中有最小的A(标准化)读取计数。方差滤波器，如IQR(四分位间范围)提供了可变性的稳健估计;可用于对变化最小的区域进行排序和丢弃。
• 更细致入微的方法:刨边机装饰图案;今天的工作流程。

贷款信息

• 为什么低统计能力会提高错误发现率?
• 培养直觉的一种方法是将t检验(例如)视为方差之比。分子是特定于治疗的，但分母是总体可变性的衡量标准。
• 用不确定度测量方差;高估或低估分母方差会对t统计量或类似比率产生不对称影响膨胀这一数据比高估数据更能显著地降低数据的价值。从而提高了错误发现率。
• 在微阵列或RNA-seq实验中使用的典型零假设下，每个基因对处理的反应可能不同(分子方差)，但基因的总体变异性是相同的，至少对于平均表达相似的基因是如此
• 策略是将分母方差估计为基因的组间方差，主持通过所有基因的组间平均方差。
• 这一策略利用了这样一个事实，即相同的实验设计已应用于所有分析的基因，并有效地调节错误发现率。

7.理解

将不同表达的区域置于上下文中

• 与基因组范围相关的基因名称
• 基因集富集及相似分析
• 接近监管标志
• 结合其他分析，例如，甲基化，拷贝数，变异，…

基因组拷贝数与mRNA表达之间的相关性在TCGA卵巢癌Affymetrix微阵列数据集中鉴定了38个错误标记的样本。

归一化

DESeq2estimateSizeFactors ()，安德斯和胡贝尔，2010

• 对于每个基因:所有样本的几何平均值。
• 每个样本:样本基因与所有样本几何平均值之比的中位数
• 可以使用除中位数以外的其他函数;控制基因可以代替

刨边机calcNormFactors ()Robinson和Oshlack的TMM方法，2010

• 确定参考样本:上四分位数最接近所有库中平均上四分位数的库
• 计算每个基因的m值(相对于参考值的对数倍变化)
• 将库大小总结为m值的加权修剪平均值。

分散

DESeq2estimateDispersions ()

• 估计每个基因的分散度
• 拟合离散度和丰度之间的平滑关系

刨边机estimateDisp ()

• 共性:所有基因单一分散;适用于小型实验(<10?样品)
• Tagwise:所有基因的分散度不同;适用于较大/表现良好的实验
• 趋势:基于丰度的仓，估计仓内常见的分散度，拟合仓内分散度和丰度之间的黄土平滑关系

例如:学习任务

t.test ()

• x:单变量测量向量
• y：因素描述实验设计
• var.equal = TRUE:适用于没有额外信息的相对较小的实验?
• 公式:替代表示法，Y ~ x．

(睡眠)

##额外组ID ## 1 0.7 1 1 ## 2 -1.6 1 2 ## 3 -0.2 1 3 ## 4 -1.2 1 4 ## 5 -0.1 1 5 ## 6 3.4 1 6

图(额外~组，数据=睡眠)

##与(sleep, t.test(extra[group == 1]， extra[group == 2]))的传统接口

## ## Welch Two样本t-检验## ##数据:额外的[组== 1]和额外的[组== 2]## t = -1.8608, df = 17.776, p-value = 0.07939 ##替代假设:真正的均值差异不等于0 ## 95%置信区间:## -3.3654832 0.2054832 ##样本估计:## x的均值y的均值0.75 2.33

t.test(额外的组，睡眠)

## ## Welch Two样本t检验## ##数据:额外的组## t = -1.8608, df = 17.776, p值= 0.07939 ##替代假设:均值的真实差异不等于0 ## 95%置信区间:## -3.3654832 0.2054832 ##样本估计:##组1的均值在组2的均值0.75 2.33

t.test(额外~组，睡眠，var.equal=TRUE)

## ##两个样本t-检验## ##数据:额外由组## t = -1.8608, df = 18, p-value = 0.07919 ##替代假设:均值的真实差异不等于0 ## 95%置信区间:## -3.363874 0.203874 ##样本估计值:##组1的均值，组2的均值，## 0.75 2.33

lm ()而且方差分析()

• lm ():适合线性模型．
• 方差分析():统计评价。

线性模型;比较t.test(var.equal=TRUE)拟合<- lm(额外~组，睡眠)方差分析(拟合)

##方差分析表## ##响应:额外## Df Sum Sq Mean Sq F value Pr(>F) ##组1 12.482 12.4820 3.4626 0.07919。##剩余18 64.886 3.6048 ##—##符号。代码:0 '***' 0.001 '**' 0.01 '*' 0.05 '。' 0.1 ' ' 1

• 引擎盖下:公式:翻译成模型矩阵，用于lm.fit ()．
• 在(隐式)截距为1时，模型矩阵的最后一个系数反映了群体效应
• 截距为0时，对比系数1和系数2的效应之间反映了群体效应

##底层模型，用于' lm.fit() '模型。矩阵(额外~组，睡眠)#最后一列表示组效应

# #(拦截)group2 # # 1 1 0 # # 2 1 0 # # 3 1 0 # # 4 1 0 0 5 # # 1 # # 6 0 1 0 7 # # 1 # # 8 1 0 0 9 # # 1 # # 10 1 0 11 # # 1 # # 12 1 1 # # 13 1 1 14 # # 1 # # 15 1 1 # # 16 1 1 # # 17 # # 18 1 1 # # 19 1 1 20 # # 1 # # attr(“转让”),# # [1]0 1 # # attr(“对比”)# # attr(“对比”)$ # #集团[1]“contr.treatment”

模型。矩阵(额外~ 0 +组，睡眠)#列之间的对比

# # group1 group2 # # 1 1 0 # # 2 1 0 # # 3 1 0 # # 4 1 1 0 # # 5 0 1 0 # # 7 # # 6 0 1 0 8 # # 1 # # 9 1 0 # # 10 1 0 0 1 # # 12 0 1 # # 11 # # 13 0 1 # # 14 15 0 0 1 # # 1 # # 16 0 1 0 1 # # 18 0 1 # # 17 # # 19 20 0 0 1 # # 1 # # attr(“转让”),# # [1]1 1 # # attr(“对比”)# # attr(“对比”)$ # #集团[1]“contr.treatment”

• 协变量-仅包含协变量的拟合基础模型，当模型包含感兴趣因素时，测试拟合的改进

fit0 <- lm(额外~ ID, sleep) fit1 <- lm(额外~ ID + group, sleep)方差分析(fit0, fit1)

##方差分析表## ##模型1:额外~ ID ##模型2:额外~ ID +组## Res.Df RSS Df Sq F Pr的总和(>F) ## 1 10 19.290 ## 29 6.808 1 12.482 16.501 0.002833 ** ##—##有意义。代码:0 '***' 0.001 '**' 0.01 '*' 0.05 '。' 0.1 ' ' 1

t.test(extra ~ group, sleep, var.equal=TRUE, paired=TRUE)

## ##配对t-检验## ##数据:额外由组## t = -4.0621, df = 9, p-value = 0.002833 ##替代假设:均值的真实差异不等于0 ## 95%置信区间:## -2.4598858 -0.7001142 ##样本估计:##差异的均值## -1.58

genefilter: rowttests ()

• 基因表达数据的t检验
• 用于探索性分析，但在统计上并不理想
• x:表达式值矩阵

• 特征x样本(与“统计学家”表示数据的方式相反——样本x特征)

• 前沿空中管制官:描述实验设计的一两个层次的因子

限制

• 假设特征是独立的
• 忽略了常见的实验设计
• 忽略多重测试

后果

• 对每个特征的组间方差估计不足
• 错误发现率上升

常用实验设计

• 学习任务:计数因数．选择:计数~ 0 +因子和对比
• 共:计数~协变量+因子
• 单因素，多水平(单向方差分析)-统计对比:指定模型为计数因数或计数~ 0 +因子
• 析因设计-主要效果，计数~因子1 +因子2；主要效应和相互作用，计数~因子1 *因子2．对比提出具体问题
• 配对设计:包括ID作为协变量(近似，因为ID是随机效应);limma方法:duplicateCorrelation ()

1.实验设计

本工作流程中使用的数据是用地塞米松(一种具有抗炎作用的合成糖皮质激素类固醇)处理的气道平滑肌细胞的RNA-Seq实验。例如，糖皮质激素被用于哮喘患者，以预防或减轻呼吸道炎症。在实验中，四株原代人气道平滑肌细胞系用1微摩尔地塞米松处理18小时。对于每一种细胞系，我们都有一个处理过的和一个未处理过的样本。实验参考如下:

Himes BE, Jiang X, Wagner P, Hu R, Wang Q, Klanderman B, Whitaker RM, Duan Q, Lasky-Su J, Nikolos C, Jester W, Johnson M, Panettieri R Jr, Tantisira KG, Weiss ST, Lu Q.“RNA-Seq转录组分析发现CRISPLD2是调节气道平滑肌细胞细胞因子功能的糖皮质激素响应基因。”公共科学图书馆，2014年6月13日;9(6):e99625。PMID:24926665．地理:GSE52778．

2、3和4:湿实验室、测序和比对

• 配对端测序导致FASTQ阅读文件及其质量分数。

• 与参考基因组或转录组排列的Reads导致BAM文件。本实验的读数与Ensembl发布的75人类参考基因组进行了比对明星调整器

5.减少

我们使用气道包装说明减少。该包提供了示例信息、八个BAM文件的子集以及计算读取所需的已知基因模型。

库(气道)路径<-系统。文件(package="气道"，"extdata")

“GSE52778_series_matrix.txt”##[2]“homo_sapiens . grch37.75 _子集。gtf“##[3]”srr1039508_子集。“##[4]”srr1039509_子集。“##[5]”srr1039512_子集。“##[6]”srr1039513_子集。“##[7]”srr1039516_子集。“##[8]”srr1039517_子集。“##[9]”srr1039520_子集。“##[10]”srr1039521_子集。SraRunInfo_SRP033351.csv" ## [12] "sample_table.csv"

csvfile <- dir(path， "sample_table.csv"， full=TRUE)

SRR1039508 GSM1275862 N61311 untrt untrt SRR1039508 126 SRX384345 ## SRR1039509 GSM1275863 N61311 trt untrt SRR1039509 126 SRX384346 ## SRR1039512 GSM1275866 N052611 trt untrt SRR1039513 126 SRX384349 ## SRR1039513 GSM1275867 N052611 trt untrt SRR1039513 87 SRX384350 ## SRR1039516 GSM1275870 N080611 trt untrt SRR1039516 # SRR1039517 GSM1275871 N080611 trt untrt SRR1039517 SRR1039508取样生物样本## #srr508568 samn02422669 ## srr1039509 srs508567 samn02422675 ## srr1039512 srs508571 samn02422678 ## srr1039513 srs508572 samn02422670 ## srr1039516 srr508575 samn02422682 ## srr1039517 srs508576 samn02422673

库(Rsamtools)文件名<- dir(路径，”。bam$"， full=TRUE) bamfiles <- BamFileList(文件名，yieldSize=1000000) names(bamfiles) <- sub(" _子集。Bam "， ""， basename(文件名))

库(基因组特征)gtffile <- file。path(path， "Homo_sapiens.GRCh37.75_subset.gtf") txdb <- makeTxDbFromGFF(gtffile, format="gtf"， circ_seqs=character())

##准备元数据帧元数据:好的

- exonsBy(txdb, by="gene")

library(GenomicAlignments) se <- summarizeOverlaps(features=基因，reads=bamfiles, mode="Union"， singleEnd=FALSE，忽略。strand=TRUE, fragments=TRUE) colData(se) <- as(sampleTable， "DataFrame") se

##类:summarizeexperiment ## dim: 20 8 ## exptData(0): ## assays(1): counts ## rownames(20): ENSG00000009724 ENSG00000116649…ENSG00000271794 ## ENSG00000271895 ## rowRanges元数据列名(0):## colnames(8): SRR1039508 SRR1039509…SRR1039520 SRR1039521 ## colData names(9): SampleName cell…样本BioSample

colData (se)

运行avgLength ##       ## SRR1039508 GSM1275862 N61311 untrt untrt SRR1039509 126 ## SRR1039509 GSM1275866 N052611 untrt untrt SRR1039512 126 ## SRR1039513 GSM1275867 N052611 trt untrt SRR1039516 120 ## SRR1039517 GSM1275871 N080611 trt untrt SRR1039517 126 #### SRR1039508 SRX384345 SRS508568 SAMN02422669 ## SRR1039509 SRX384346 SRS508567 SAMN02422675 ## SRR1039513 SRX384349 SRS508571 SAMN02422678 ## SRR1039516 SRX384353 SRS508575 SAMN02422682 ## SRR1039517 SRX384354 SRS508576 SAMN02422673 ## SRR1039520 SRX384357 SRS508579Samn02422683 ## srr1039521 srx384358 srs508580 samn02422677

rowData (se)

##警告:'rowData'已弃用。##使用“rowRanges”代替。##查看帮助(“已弃用”)

##长度为20的GRangesList对象## $ENSG00000009724## seqnames ranges strand | exon_id exon_name ##    |   ## [1] 1 [11086580,11087705] - | 98 ENSE00000818830 ## [2] 1 [11090805,11090939] - | 100 ENSE00000743084 ## [4] 1 [11094885,11094963] - | 101 ENSE00000743085 ## [5] 1 [11097750,11097868] - | 103 ENSE00003520086 ## ... ... ... ... ... ... ...## | 111 ense00003467404 ## [15] 1 [11106948, 11107176] - | 112 ense00003489217 ## [16] 1 [11107260, 11107280] - | 113 ense00001833377 ## [17] 1 [11107260, 11107284] - | 114 ense00001472289 ## [18] 1 [11107260, 11107290] - | 115 ense00001881401 ## ##…## <19个更多的元素> ## ------- ## seqinfo: 1个来自未指定基因组的序列;没有seqlengths

(试验(se))

# # SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 24 42 # # # # 28 ENSG00000009724 66 ENSG00000116649 1004 1255 1122 1313 1100 218 256 233 252 269 # # ENSG00000120942 # # ENSG00000120948 2751 2080 3353 1614 3519 # # ENSG00000171819 4 50 19 543 1 # # ENSG00000171824 869 1075 1115 1051 944 # # SRR1039517 SRR1039520 SRR1039521 # # ENSG00000009724 41 47 36 # # ENSG00000116649 1879 745 1536 # # ENSG00000120942 465 207 400 # # ENSG00000120948 3716 2220 1990 1067 # # # # ENSG00000171819 10 14Ensg00000171824 1405 748 1590

6.分析使用DESeq2

前一节说明了数据子集上的约简步骤;这是完整的数据集

资料(气道)se <-气道

该对象包含一个信息类colData槽位-按照手册中的描述准备气道装饰图案。特别是colData ()包括描述单元系的列细胞和治疗敏捷对于每个样本

colData (se)

运行avgLength ##       ## SRR1039508 GSM1275862 N61311 untrt untrt SRR1039509 126 ## SRR1039509 GSM1275866 N052611 untrt untrt SRR1039512 126 ## SRR1039513 GSM1275867 N052611 trt untrt SRR1039516 120 ## SRR1039517 GSM1275871 N080611 trt untrt SRR1039517 126 #### SRR1039508 SRX384345 SRS508568 SAMN02422669 ## SRR1039509 SRX384346 SRS508567 SAMN02422675 ## SRR1039513 SRX384349 SRS508571 SAMN02422678 ## SRR1039516 SRX384353 SRS508575 SAMN02422682 ## SRR1039517 SRX384354 SRS508576 SAMN02422673 ## SRR1039520 SRX384357 SRS508579Samn02422683 ## srr1039521 srx384358 srs508580 samn02422677

DESeq2使分析特别容易，只需添加实验设计，运行管道，并提取结果

library(DESeq2) dds <- DESeqDataSet(se, design = ~ cell + dex) dds <- DESeq(dds)

##估计大小因素##估计分散##基因分散估计##均值-分散关系##最终分散估计##拟合模型和测试

Res <-结果(dds)

简单的可视化/健全检查包括

• 看看强分化基因的计数，了解计数是如何转化为结果表中报告的汇总统计数据的

  topGene <- rownames(res)[where .min(res$padj)] res[topGene，]

  ## log2 fold change (MAP): dex untrt vs trt# # Wald test p-value: dex untrt vs trt# #数据帧1行6列## baseMean log2FoldChange lfcSE stat pvalue ## <数字> <数字> <数字> <数字> <数字> ## ENSG00000152583 997.4398 -4.3161 0.1724125 -25.03357 2.636198e-138 ## padj ## <数字> ## ENSG00000152583 4.624155e-134

  plotCounts(dds, gene=topGene, intgroup=c("dex"))

  

• 一个“MA”图显示每个基因的组间对数折叠变化与平均对数计数的关系;它应该是漏斗形，周围近似对称y = 0，对于低计数的基因，在处理之间有很多变化。

  plotMA (res ylim = c (5,5))

  

• 绘制(未调整的)P值的分布，它应该是均匀的(在null下)，但在小P值处有一个峰值(希望是真阳性!)

  嘘(res pvalue美元,减免= 50)

  

• 观察调整后的p值与对数折叠变化的“火山图”，以了解上调与下调基因的比例

  plot(-log10(padj) ~ log2FoldChange, as.data.frame(res)， pch=20)

  

许多额外的诊断方法在DESeq2(和edgeR)小插图中描述，在RNA-seq基因差异表达工作流程中描述。

7.理解

参见E部分，基因集富集