# 整理表型文件
pdata_col <- pdata %>% apply(2,function(x){unique(x) %>% length()}) %>% data.frame() %>% filter(.>1) %>% rownames() 
pdata <- pdata %>% select(pdata_col)

library(DESeq2)

## 3.1表达矩阵


## 3.2分组矩阵
condition = pdata$`sample type:ch1`

# 配对分析要加上这段代码，知道谁和谁是一对，比如1,1是一对，5,5是一对,需要手动写出
# subject <- factor(c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5))  

coldata <- data.frame(row.names = colnames(exprSet), condition)
coldata$condition <- as.factor(coldata$condition)

# 注意在design中加上配对信息,这里没加
dds <- DESeqDataSetFromMatrix(countData = exprSet,
                              colData = coldata,
                              design = ~condition)

nrow(dds)
rownames(dds) %>% head()
dds <- dds[rowSums(counts(dds)) > 1, ] 
nrow(dds)

dds$condition<- relevel(dds$condition, ref = "normal") # 指定哪一组作为对照组

## 3.3差异表达矩阵

dds <- DESeq(dds)  
res <- results(dds) %>% data.frame()
# 正确写法是：实验组在前，对照组在后
res1 <- results(dds,contrast = c("condition","diseased","normal")) %>% data.frame()                     
res2 <- results(dds,contrast = c("condition","normal","diseased")) %>% data.frame()



# 定义加上下调标记函数,适用于DEseq2包结果
# 用p
add_sign_p <- function(df){
  df <- df %>% mutate(sign = dplyr::case_when(df$pvalue < p_cutoff & df$log2FoldChange> logFC_cutoff ~ "up",
                                              df$pvalue < p_cutoff & df$log2FoldChange< -logFC_cutoff ~ "down",
                                              TRUE ~ "stable")
  )
  return(df)
}
# 用padj
add_sign <- function(df){
  df <- df %>% mutate(sign = dplyr::case_when(df$padj < p_cutoff & df$log2FoldChange> logFC_cutoff ~ "up",
                                              df$padj < p_cutoff & df$log2FoldChange< -logFC_cutoff ~ "down",
                                              TRUE ~ "stable")
  )
  return(df)
}


logFC_cutoff <- 2
p_cutoff <- 0.05

res1 <- add_sign(res1)
# res2 <- add_sign(res2)
# res3 <- add_sign(res3)

diffSig1 <- res1 %>% filter(sign == "up" | sign == "down")