数据挖掘—DESeq2差异基因分析将部分基因padj设置为NA

rm(list = ls())
library(dplyr)
library(openxlsx)
library(DESeq2)
#大鼠
exp  <- read.xlsx("../data/Rattus norvegicus/Expression_with_annotation.xlsx")
colnames(exp)

select_name <- c("Name","cg_1:read.count","cg_2:read.count","cg_3:read.count",
                 "cg_Icariin_1:read.count","cg_Icariin_2:read.count","cg_Icariin_3:read.count")


exp <- exp[,select_name]
dim(exp) #22016     7

exp[1:2,]

table(duplicated(exp$Name))


#数据清洗
# ---------- 1. 去除重复基因（保留总表达量最高的记录） ----------
exp <- exp %>%
  mutate(TotalCount = rowSums(across(-Name))) %>%
  arrange(desc(TotalCount)) %>%
  distinct(Name, .keep_all = TRUE) %>%
  select(-TotalCount)
exp <- as.data.frame(exp)
rownames(exp) <- exp$Name
exp <- exp[, -1]
colnames(exp) <- c(
  "Control1", "Control2", "Control3",
  "Icariin1", "Icariin2", "Icariin3"
)
exp <- as.matrix(exp)
storage.mode(exp) <- "integer"
dim(exp)#21812     6
head(exp)
sum(duplicated(rownames(exp)))


#过滤
keep_feature <- rowSums (exp > 1) > 1 
table(keep_feature)
exp_clean <- exp[keep_feature, ]  
dim(exp_clean)#7767    6
colnames(exp_clean)
class(exp_clean)
exp_clean <- as.data.frame(exp_clean)


#分组信息
library(stringr)
k = str_detect(colnames(exp),"Control")
table(k) 
Group = ifelse(k,"Control","Icariin")
Group = factor(Group,levels = c("Control","Icariin"))
Group
# [1] Control Control Control Icariin Icariin Icariin
# Levels: Control Icariin


#差异基因分析
group_list <- Group
exprSet = exp_clean

(colData <- data.frame(row.names=colnames(exprSet), 
                       group_list=group_list) )

exprSet[, 1:ncol(exprSet)] <- lapply(exprSet[, 1:ncol(exprSet)], as.integer)

class(exprSet)
dds <- DESeqDataSetFromMatrix(countData = exprSet,
                              colData = colData,
                              design = ~ group_list)
dds <- DESeq(dds) 
res <- results(dds, 
               contrast=c("group_list",
                          levels(group_list)[2],
                          levels(group_list)[1]))
resOrdered <- res[order(res$padj),]
head(resOrdered)
DEG =as.data.frame(resOrdered)
dim(DEG)#7767    6
deg = na.omit(DEG)


#差异基因分析
logFC_t = 0.585
p_t = 0.05

k1 = (deg$pvalue < p_t)&(deg$log2FoldChange < -logFC_t)
k2 = (deg$pvalue < p_t)&(deg$log2FoldChange > logFC_t)
deg = mutate(deg,change = ifelse(k1,"Down",ifelse(k2,"Up","Not")))
table(deg$change)
# Down  Not   Up 
# 13 2469   15 

> dim(DEG)#7767    6
[1] 7767    6
> deg = na.omit(DEG)
> dim(deg)
[1] 2497    6
> # 1. 原始表达矩阵基因数
> nrow(exp_clean)
[1] 7767
> # 2. DESeq2完整结果基因数
> nrow(DEG)
[1] 7767
> # 3. na.omit后剩余基因数
> nrow(na.omit(DEG))
[1] 2497
> # 4. 分别检查各列NA数量
> colSums(is.na(DEG))
      baseMean log2FoldChange          lfcSE           stat         pvalue           padj 
             0              0              0              0              0           5270 
> # 5. 检查P值和校正P值
> table(is.na(DEG$pvalue))

FALSE 
 7767 
> table(is.na(DEG$padj))

FALSE  TRUE 
 2497  5270 

#使用原始P值筛选
logFC_t = 0.585
p_t = 0.05

k1 = (deg$pvalue < p_t) &
     (deg$log2FoldChange < -logFC_t)

k2 = (deg$pvalue < p_t) &
     (deg$log2FoldChange > logFC_t)

> deg <- DEG
> dim(deg)
[1] 7767    6

...

logFC_t = 0.585
p_t = 0.05

k1 = (deg$pvalue < p_t)&(deg$log2FoldChange < -logFC_t)
k2 = (deg$pvalue < p_t)&(deg$log2FoldChange > logFC_t)
deg = mutate(deg,change = ifelse(k1,"Down",ifelse(k2,"Up","Not")))
table(deg$change)
# Down  Not   Up 
# 77 7615   75