16s分析之差异分析（DESeq2）

今天我们来学习R语言DESeq2包，原理什么的后不说，在操作过程中点缀一下，等四个差异包推送完成后，咱们在对这四个包做差异分析的原理做一个比较分析；

这个包适用于：

高通量数据分析过程中，基于count数据，对其进行标准化处理，并对两个样本的差异做定量比较

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

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

biocLite("DESeq2")

library(DESeq2)

读取我们的数据

# 读入实验设计，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)

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