搜索
查看: 2822|回复: 0

[software] Count Based RNAseq 就是差异表达的两个包的简介

[复制链接]

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-7-11 08:54:13 | 显示全部楼层 |阅读模式
[AppleScript] 纯文本查看 复制代码
```{r dependencies, warning=FALSE, message=FALSE}
library(DESeq2)
library(edgeR)
```


Use the following commands to install these packages in R.

[AppleScript] 纯文本查看 复制代码
```{r biocLite, eval=FALSE}
source("http://www.bioconductor.org/biocLite.R")
biocLite(c("DESeq2", "edgeR))
```


## Overview

通常通过创建每个样品的基因计数的计数矩阵来分析RNA seq数据。 该矩阵使用基于计数的模型进行分析,通常建立在负二项分布上。 包含这个模型的包有“r Biocpkg(”edgeR“)`和`r Biocpkg(”DESeq“)`/`r Biocpkg(”DESeq2“)`。

这种类型的分析丢弃了RNA排序读取中的部分信息 。
## RNA-seq count data

分析RNA测序数据的一种简单方法是使其看起来像微阵列数据。这是通过计算每个样本中的几个读段与基因重叠来完成的。有很多方法可以做到这一点。它显然取决于使用的注释,还取决于如何确定读段与区域重叠。特别关注的是哪些基因组区域是具有多个转录本的基因的一部分。

对这个过程没有一致意见,并且知道影响结果的不同选择。

进行基因计数的工具包括

- 通过使用`r Biocpkg(“Rsubread”)`包中的`featureCounts()`。
- [HTSeq](http://www-huber.embl.de/users/anders/HTSeq/)包(这是一个python包,而不是Bioconductor包)。
- 通过使用`r Biocpkg(“GenomicAlignments”)`包中的`summarizeOverlaps()`。


将RNA测序数据减少到每个基因的单个整数是一个简化。实际上,它忽略了进行RNA测序的一些主要原因,包括评估选择性剪接。另一方面,我们很好地了解这个程序的统计特性,并且能够深入了解大多数研究者想要了解的内容。最后,这种方法需要预先知道不同的基因组区域

## Statistical issues

在RNA-seq数据分析中,我们经常看到许多基因(高达50%)很少或没有表达。 在分析前预先过滤(除去)这些基因是很常见的。 一般来说,基因组过滤可能会对您的分析有所帮助,但此讨论超出了本文档的范围。

## The Data


我们将使用“Biocexptpkg(”airway“)”数据集,其中包含“SummarizedExperiment”形式的RNA-seq数据。
[AppleScript] 纯文本查看 复制代码
```{r data}
library(airway)
data(airway)
airway
assay(airway, "counts")[1:3, 1:3]
airway$dex
```

感兴趣的主要变量是“dex”,它具有“trt”(被处理)和“untrt”(未处理)的级别。 第一级将是此因子的参考级别,因此我们使用`relevel()`来设置`untrt`级别作为引用; 这很容易解释。
[AppleScript] 纯文本查看 复制代码
```{r relevel}
airway$dex <- relevel(airway$dex, "untrt")
airway$dex
```


There is rich information about which gene model was used for each gene:
有关每种基因使用哪种基因模型的丰富信息:
```{r granges}
granges(airway)
```

## edgeR

`r Biocpkg(“edgeR”)`在“r Biocpkg(”limma“)`的数据结构和功能方面非常相似。 而`r Biocpkg(“limma”)`允许我们在“ExpressionSet”上直接操作,edgeR不能直接与`SummarizedExperiment`一起工作。 我们首先需要将数据放入edgeR特定的容器中。
[AppleScript] 纯文本查看 复制代码
```{r edgeRsetup}
library(edgeR)
dge <- DGEList(counts = assay(airway, "counts"),
               group = airway$dex)
dge$samples <- merge(dge$samples,
                     as.data.frame(colData(airway)),
                     by = 0)
dge$genes <- data.frame(name = names(rowRanges(airway)),
                        stringsAsFactors = FALSE)
```


这个对象有一个叫做“group”的东西,它是每个样本的基本实验组。 当您创建“DGEList”对象时,它还具有“$ samples”(pheno数据 )

设置输入对象后,我们现在进行如下操作。

首先我们估计标准化因素或有效的library 规模
[AppleScript] 纯文本查看 复制代码
```{r calcNormFactors}
dge <- calcNormFactors(dge)
dge
```


接下来,我们设置设计矩阵并估计 (方差)。 有多种方法可以做到这一点, 有两步是必须的
[AppleScript] 纯文本查看 复制代码
```{r disp}
design <- model.matrix(~dge$samples$group)
dge <- estimateGLMCommonDisp(dge, design)
dge <- estimateGLMTagwiseDisp(dge, design)
```


`glmFit()`, 类似 `r Biocpkg("limma")`

[AppleScript] 纯文本查看 复制代码
```{r edgeRdesign}
fit <- glmFit(dge, design)
```


做检验并且选出 top hits

[AppleScript] 纯文本查看 复制代码
```{r glmLRT}
lrt <- glmLRT(fit, coef = 2)
topTags(lrt)
```


## DESeq2


像`r Bioc Package(“edgeR”)`,DESeq2要求我们将数据放入特定于程序包的容器(“DESEQ DataSet”)中。 但是,与edgeR不同,这很简单。
[AppleScript] 纯文本查看 复制代码
```{r DESeq2setup}
library(DESeq2)
dds <- DESeqDataSet(airway, design = ~ dex)
```


请注意,实验的设计存储在对象内。 最后一个变量(在多个变量为列表的情况下)将是不同结果输出中报告的感兴趣的变量。

模型很简单
[AppleScript] 纯文本查看 复制代码
```{r deseqfit}
dds <- DESeq(dds)
```


然后我们需要做的就是得到结果。 请注意,结果没有排序的,所以下面用order。
[AppleScript] 纯文本查看 复制代码
```{r deseqResults}
res <- results(dds)
res <- res[order(res$padj),]
res
```




回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-16 09:12 , Processed in 0.027822 second(s), 24 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.