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)
学习永无止境，分享永不停歇！