Seurat DEG로 MA plot 그리기

오늘은 Seurat에서 FindMarker로 DEG구하고 결과를 MA plot으로 그려보려함.

Volcano plot이랑 다르게 mean expression을 고려할 수 있어서 fold change가 같더라도 더 발현이 많이 되는 유전자는 골라낼 수 있어 유용함.

Volcano plot은 p-val 와 foldchange여서 필요한 값을 FindMarker() 돌리면서 한번에 얻을 수 있었는데, MA plot그리기 위해서는 averge expression이 필요해서.. seurat의 AverageExpression() 펑션으로 따로 구해줘야함

# Seurat 객체에서 평균 발현량 계산
Idents(seurat.obj) <- "orig.ident"
avg_expression <- AverageExpression(seurat.obj)$RNA

# DEG 찾기
Idents(seurat.obj) <- "Condition"
de <- FindMarkers(seurat.obj, ident.1 = "treated", ident.2 = "ctrl", test.use = "MAST", min.pct = 0, logfc.threshold = -Inf)

# 필요한 데이터 결합
# 이 부분은 데이터 구조에 따라 조정이 필요할 수 있습니다.   avgExp <- AverageExpression(subset, slot="data")$RNA
merged_data <- merge(de, avg_expression@assays$RNA@data, by = "row.names", all.x = TRUE)
colnames(merged_data)[which(colnames(merged_data) == "Row.names")] <- "gene"
colnames(merged_data)[7] <- "avgExp"
merged_data$logAvgExp <- log2(merged_data$avgExp+1)

# MAplot 준비
merged_data$significance <- 'not significant'
merged_data$significance[merged_data$p_val_adj < 0.05 & merged_data$avg_log2FC > 0] <- 'upDEG'
merged_data$significance[merged_data$p_val_adj < 0.05 & merged_data$avg_log2FC < 0] <- 'downDEG'
color <- c("upDEG" = "red", "downDEG" = "#00B2FF", "not significant" = "darkgrey")
merged_data$delabel <- merged_data$gene
merged_data[merged_data$significance=="not significant",'delabel'] <- NA


# MAplot 그리기
g <- ggplot(data=merged_data, aes(x=logAvgExp, y=avg_log2FC, color=significance, label=delabel))+
        ggtitle(label = paste0("treated-ctrl ",i," MA Plot: MAST")) +
        geom_point(alpha=0.7, size=2)+
        geom_text_repel()+
        geom_hline(aes(yintercept=0), colour = "black", size = 0.5) +
      ylim(c(min(merged_data$avg_log2FC), max(merged_data$avg_log2FC))) + 
      xlab("Log2 Mean expression") + 
      ylab("Log2 Fold Change") + 
      theme(axis.title.x = element_text(face = "bold", size = 15),
            axis.text.x = element_text(face = "bold", size = 12)) +
      theme(axis.title.y = element_text(face = "bold", size = 15),
            axis.text.y = element_text(face = "bold", size = 12)) +
        scale_color_manual(values=color) +
      theme(legend.title = element_text(face = "bold", size = 15)) +
      theme(legend.text = element_text(size = 14))+
        theme_minimal()
ggsave("MAplot.png")

Results

예시가 못생겼지만.. 대충 이렇게 그려짐 어휴