miRNA seq差异表达分析练习（二）——DESeq差异表达分析

数据准备

接上一篇《miRNA seq差异表达分析练习（一）——GEO样本数据下载》，下载的数据保存在文件夹miRNA_seq中

有60个txt.gz文件

从后缀名可以看出，这些都是压缩文件，直接打开是乱码的，解压后打开就可以看到具体内容了

事实上，用R语言是可以直接读取txt.gz文件而不需要解压的。

接下来，我将用R语言批量读取该文件夹下所有的txt.gz文件，并合并成一个包含60个样本数据的数据框。

# 读取D:\用户\桌面\\miRNA_seq文件夹下所有文件
setwd('D:\\用户\\桌面\\miRNA_seq')
fileName <- dir()  # 获取所有文本名

# 先读取第一个文件
mydata <- read.table(fileName[1], header = T, row.names = 1)
names(mydata) <- c('sample_1')

# 读取其余文件，并逐一并到第一个文件上
for(i in 2:length(fileName)){
  data <- read.table(fileName[i],header = T, row.names = 1)
  names(data) <- c(paste('sample',as.character(i),sep = '_'))
  mydata <- cbind(mydata,data)
}
names(mydata) <- c(paste('N-',1:20,sep = ''),paste('T-',1:20,sep = ''),paste('P-',1:20,sep = ''))

# 对所有数据四舍五入取整
mydata <- round(mydata,0)

# 查看mydata数据框前6行
head(mydata)

合并后的数据框mydata，行名即miRNA名称，列名为各样本名称。

因为DESeq包要求的表达矩阵数据为整数，所以在正式分析前，应将所有数据四舍五入取整数，一行代码搞定！

 

# 对所有数据四舍五入取整
mydata <- round(mydata,0)  # 0表示保留0给小数，即整数

DESeq差异表达分析

差异表达分析需要两个数据：表达矩阵和分组信息

表达矩阵简单获得，只需将前面的数据框转化为矩阵即可。

mymatrix <- as.matrix(mydata)

分组信息也很简单，60个样本，三个组织即分为3组，每组20个,其中,N组是正常组织,T和P组都属于癌变组织。

condition <- factor(c(rep('normal',20),rep('cancerous',40)))
org <- factor(c(rep('N',20),rep('T',20),rep('P',20)))
coldata <- data.frame(row.names = colnames(mydata),org)

OK，万事具备，差异分析走起！

# 差异表达分析
library(DESeq2)
dds <- DESeqDataSetFromMatrix(mymatrix,DataFrame(coldata), design = ~org)
# 标准化
dds = DESeq(dds)
# 提取结果
res <- results(dds,alpha = 0.001)
res = res[order(res$padj),]
summary(res) #查看结果
# 输出所有结果
# write.csv(res,file="D:\\用户\\桌面\\All_results.csv")

搞定！差异表达分析似乎也没有想象中那么难哈，下一篇进行结果分析与保存。

啰嗦一句，DESeq2包与其他包的安装有些不同，不过方法已经给你备好了！

source("http://bioconductor.org/biocLite.R")

biocLite("DESeq2")

安装时间还是挺久的，不过，磨刀不误砍柴工，一次安装，永久使用！