################################
## 2020/10/05  DES analysis

setwd("E:/RNAseq/")
library(DESeq2)
library(apeglm)
library(ashr)

# 1. 构建dds
mycounts <- read.table(file = paste0( "gene_all.counts.matrix"), 
                         header = T, 
                         sep = "\t", 
                         quote = "", 
                         check.names = F)

head(mycounts)
rownames(mycounts) <- mycounts[,1]
mycounts <- mycounts[,-1] ## 这里要注意一定要删掉第一列第一行的那个空格或文字

condition <- factor(c(rep("h",3),rep("L",3)), levels = c("h","L"))

colData <- data.frame(row.names = colnames(mycounts), condition)

# 2. DESeq流程
dds <- DESeqDataSetFromMatrix(mycounts, colData, design= ~ condition)
dds <- DESeq(dds)

res <- results(dds)
head(res)
head(as.data.frame(res))

# sort by p.adj
res <- res[order(res$padj),]
# get diff_gene ，导出差异基因
diff_gene <- subset(res, padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))

resdata <- merge(as.data.frame(res), 
                 as.data.frame(counts(dds, normalized=F)), 
                 by="row.names", sort=FALSE)

resdata$significant <- "unchanged"
resdata$significant[resdata$padj <= 0.05 & resdata$log2FoldChange >= 1 ] <- "upregulated"
resdata$significant[resdata$padj <= 0.05 & resdata$log2FoldChange <= -1 ] <- "downregulated"

head(resdata)
write.csv(resdata,file= "resdata",row.names = F)
write.csv(diff_gene, file = "diff_gene.csv",row.names = T)