通过RSEM我们获取了样本中每个基因的counts和表达量，接下来使用tximport校正不同样本间基因长度的差异。
代码语言：javascript复制## 安装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分隔。
准备好输入文件后，开始执行代码：
代码语言：javascript复制## 加载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进行差异表达分析。
代码语言：javascript复制## 安装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来设置比较的样本。
代码语言：javascript复制## 输出 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差异表达分析实操了
