转录组分析 | 使用DESeq2进行基因差异表达分析

生信小王子

发布于 2020-08-10 17:17:29

2.8K0

发布于 2020-08-10 17:17:29

通过RSEM我们获取了样本中每个基因的counts和表达量，接下来使用tximport校正不同样本间基因长度的差异。

## 安装R包
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("tximport")
## 加载R包
library("tximport")

安装好R包后，准备两个输入文件：

sample.txt：第一列为样本名，第二列为处理方式，以制表符Tab分隔。

gene_trans.txt：第一列为转录本ID，第二列为基因ID，以制表符Tab分隔。

准备好输入文件后，开始执行代码：

## 加载R包
library("tximport")
## 导入文件
samples <- read.table("./sample.txt", header = T)
tx2gene <- read.table("./gene_trans.txt",sep="\t",header = T)
## 导入表达量文件
files <- paste(samples$Sample,".isoforms.results",sep="")
## 校正样本间基因长度的差异
txi.rsem <- tximport(files, type = "rsem", tx2gene = tx2gene,countsFromAbundance = c("lengthScaledTPM"))

接下来使用DESeq2进行差异表达分析。

## 安装R包
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("DESeq2")
## 加载R包
library("DESeq2")
## 导入数据
dds <- DESeqDataSetFromTximport(txi.rsem, colData = samples, design = ~ Treatment)
## 过滤低表达基因
dds <- dds[rowSums(counts(dds)) > 1,]
## 进行差异表达分析
dds <- DESeq(dds)

完成差异表达分析后，通过results函数提取结果。

通过contrast来设置比较的样本。

## 输出 CK 与 H30 处理样本的差异表达结果
res_ck_30 <- results(dds, contrast = c("Treatment","H30","CK"))
ouf <- paste("CK_30.txt",sep ="\t")
write.table(res_ck_30,ouf)

获得差异表达分析结果后，就可以根据我们的需求制定标准筛选差异表达基因啦！

参考资料：

http://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

是时候来一波RNA-Seq差异表达分析实操了

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2020-02-13，如有侵权请联系 cloudcommunity@tencent.com 删除