高通量测序数据差异分析（DESeq2） / 开普饭

今天我们来学习R语言DESeq2做差异分析，第一次我推送这个差异分析到现在已经过去一年多了，我重新排版，更加了一些感悟，重新推送给大家。基于高通量测序数据的差异分析，为了矫正测序平台，批次，深度等差异，所以这里做DEsep2差异分析。

这个包适用于：

高通量数据分析过程中，基于count数据，对其进行标准化处理，并对两个样本的差异做定量比较。注意这里是基于read数进行分析的，卖的是标准化方法。

说道安装DEsep2这个包，这段时间在ubuntu中安装更加困难，还好最后都解决了。这个包绝对算得上的较难安装的包之一。当然现在我目前安装最为困难的包是microbiomeSeq，大家有时间可以试试，做扩增子的同学可以看看这个包，是有作者的心意在里面的。

#安装这个包，安装过程不太顺利，请多试几次，还是没问题的，我也没有出现意料之外的问题

source("https://bioconductor.org/biocLite.R") biocLite("DESeq2") library(DESeq2)

读取我们的数据，这里使用两个文件，一个是分组文件，另一个是OTU表格。

有的人曾问过，这个OTU表格需要抽平吗？这里我直接回答：不需要。

# 读入实验设计，Qiime的mapping文件删去第一行的“#”即可使用


design =read.table("map_lxdjhg.txt", header=T, row.names= 1,sep="\t")
head(design)
# 读取OTU表,全部otu表没有抽平，基于count的数据，不可用相对丰度数据
otu_table =read.delim("otu_table.txt", row.names= 1,sep="\t",header=T,check.names=F)
#做判别，做排序为了使变量可以一一对应上
idx =rownames(design) %in% colnames(otu_table)
design = design[idx,]
#我们后面利用sub_design文件做分组文件
sub_design =design[idx,]
#我们后面利用count文件做差异分析文件
count =otu_table[, rownames(sub_design)]
#这里我们将count转化为矩阵，当然可以不用转化
count=as.matrix(count)
head(count)
第一步构建两个矩阵，第一个OTU矩阵，第二个分组矩阵
dds <-DESeqDataSetFromMatrix(countData = count,
colData =sub_design,

design = ~ SampleType)

countData：即为差异分析otu表格，格式为行名otu名称，列名变量名为样品名

colData：即为分组文件；

design：即为分组文件中指定分组的列的列名；

第二步

dds2 <-DESeq(dds) #直接用DESeq函数即可 resultsNames(dds2) #这里我们指定需要比较的两组样本 # 将结果用results()函数来获取，赋值给res变量,这里我设定显著性阈值为alpha=0.05 res <- results(dds2,contrast=c("SampleType","G0", "GC1"),alpha=0.05)


#看一下结果的概要信息

summary(res)

提取结果进行保存

#这里我们可以更据结果文件提取我们需要的行，使用subset函数，也是也个很有用的函数


WT <-subset(res,padj< 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
#当然可以只得到想要的OTU名称

WT2 <-row.names(WT)

# res就是最后的显著性比较的文件，这里我们按照显著性排序
res <-res[order(res$padj),]
#将差异和OTU相对丰度做一个合并
wt3 <- merge(as.data.frame(res),as.data.frame(counts(dds2,normalize=TRUE)),by="row.names",sort=FALSE)

head(wt3)

# 得到csv格式的差异表达分析结果

write.table(wt3,"otu差异统计表格GC1-G0.txt",quote= FALSE,row.names = T,

col.names = T,sep = "\t")这里提到了相对丰度提取的命令 ## 获得normalize的count dds <-estimateSizeFactors(dds) wt6 <-counts(dds, normalized=T) head(wt6)

这是上面我们合并采用的命令 counts(dds2,normalize=TRUE)

学习永无止境，分享永不停歇！

高通量测序数据差异分析（DESeq2）

相关推荐