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

library(DESeq2)

## 3.1表達(dá)矩陣


## 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差異表達(dá)矩陣

dds <- DESeq(dds)  
res <- results(dds) %>% data.frame()
# 正確寫法是：實(shí)驗(yàn)組在前，對照組在后
res1 <- results(dds,contrast = c("condition","diseased","normal")) %>% data.frame()                     
res2 <- results(dds,contrast = c("condition","normal","diseased")) %>% data.frame()



# 定義加上下調(diào)標(biāo)記函數(shù),適用于DEseq2包結(jié)果
# 用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")