DESeq2 差异表达分析

library(airway)
data(airway)
exprSet=assay(airway)
group_list=colData(airway)[,3]
group_list =c("untrt", "trt", "untrt", "trt","untrt","trt","untrt", "trt")

my-R/10-RNA-seq-3-groups/hisat2_mm10_htseq.R https://github.com/jmzeng1314/my-R/blob/master/10-RNA-seq-3-groups/hisat2_mm10_htseq.R

suppressMessages(library(DESeq2)) 
(colData <- data.frame(row.names=colnames(exprSet), group_list=group_list) )
dds <- DESeqDataSetFromMatrix(countData = exprSet,
                              colData = colData,
                              design = ~ group_list)
dds <- DESeq(dds)

res <- results(dds, contrast=c("group_list","treat_2","control"))
resOrdered <- res[order(res$padj),]
head(resOrdered)
DEG_treat_2=as.data.frame(resOrdered)
write.csv(DEG_treat_2,"DEG_treat_2_deseq2.results.csv")

volcano plot

GEO/airway_RNAseq/run_DEG_RNA-seq.R
https://github.com/jmzeng1314/GEO/blob/master/airway_RNAseq/run_DEG_RNA-seq.R

logFC_cutoff <- with(need_DEG,mean(abs( log2FoldChange)) + 2*sd(abs( log2FoldChange)) )
# logFC_cutoff=1

need_DEG$change = as.factor(ifelse(need_DEG$pvalue < 0.05 & abs(need_DEG$log2FoldChange) > logFC_cutoff,
                                   ifelse(need_DEG$log2FoldChange > logFC_cutoff ,'UP','DOWN'),'NOT')
)
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
                    '\nThe number of up gene is ',nrow(need_DEG[need_DEG$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(need_DEG[need_DEG$change =='DOWN',])
)
library(ggplot2)
g = ggplot(data=need_DEG, 
           aes(x=log2FoldChange, y=-log10(pvalue), 
               color=change)) +
  geom_point(alpha=0.4, size=1.75) +
  theme_set(theme_set(theme_bw(base_size=20)))+
  xlab("log2 fold change") + ylab("-log10 p-value") +
  ggtitle( this_tile ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
  scale_colour_manual(values = c('blue','black','red')) ## corresponding to the levels(res$change)
print(g)
ggsave(g,filename = paste0(n,'_volcano.png'))

DESeq2结果p-value和padj设为NA的理由：https://blog.csdn.net/linkequa/article/details/83116789

DESeq2的baseMean和log2FoldChange是如何得到的？ https://www.jianshu.com/p/2ece602d8519

转录组入门（7）：差异表达分析 https://www.jianshu.com/p/5f94ae79f298
简单使用DESeq2/EdgeR做差异分析 https://www.bioinfo-scrounger.com/archives/113/
RNA-seq(7): DEseq2筛选差异表达基因并注释(bioMart) https://www.jianshu.com/p/3a0e1e3e41d0
DESeq2差异分析 https://www.jianshu.com/p/a27dce71f6ea

DESeq2标准化
[转载] DEseq2归一化 https://www.jianshu.com/p/562c4be23b1d

library(DESeq2)
lib.size<-estimateSizeFactorsForMatrix(mydata)
ed<-t(t(mydata)/lib.size)

或者，从差异表达分析结果中提取

normlized_counts <- counts(dds, normalized=TRUE)