代码展示:
asv <- read.delim(paste0(input,'_0.5wen.10050.asv_table.txt'), row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
group <- read.delim(paste0('group2_',input,'.txt'),row.names = 1,sep = '\t')
asv <- asv+1
#将变量转换为因子型
group$group1 <- as.factor(group$group1)
dds <- DESeqDataSetFromMatrix(countData = asv,
colData = group,
design = ~group1)#构建 DESeqDataSet 对象
dim(dds)
##过滤
#dds <- dds[rowSums(counts(dds)) > 1,]
## 差异比较
#dds <- DESeq(dds) #差异分析dds <- DESeq(dds) #asv差异分析
dep <- DESeq(dds)
res <- results(dep)
suppressMessages(dds)
diff = res
diff <- na.omit(diff) ## 去除缺失值NA
dim(diff)
#此时,diff就是差异分析的总分析结果,我们可以先将其导出保存
write.csv(diff,paste0(input,'_all_diff.csv'))
##进一步筛选差异ASV,使用Padj值和log2FC进行筛选
#Padj是P值矫正之后的数值,一般选取小于等于0.05(显著差异)的基因;同时log2FC是基因表达量的差异倍数。
#例如log2FC为1,证明这个基因在两种不同处理中的表达量相差了一倍,通常以大于1或小于-1为标准,大于1的为上调表达,少于-1的为下调表达。
foldChange = 1
padj = 0.05
diffsig <- diff[(diff$pvalue < padj & abs(diff$log2FoldChange) > foldChange),]
dim(diffsig)
write.csv(diffsig, paste0(input,'_all_diff.sig.csv'))
参考:
[1]DESeq2筛选差异OTU及绘制火山图
[2]RStudio | 利用DESeq2包进行差异表达分析
[3]NC:基于38组数据评估14种微生物组差异分析方法—结果差异很大