ggplot(de_df, aes(x = log2FC, y = -log10(padj), color = sig)) +
  geom_point(size = 0.6, alpha = 0.5) +
  geom_text_repel(data = top, aes(label = symbol), size = 3) +
  scale_color_manual(values = c(Up="#dc2626", Down="#2563eb", NS="#9ca3af")) +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
  theme_classic()

ggplot(ma_df, aes(x = baseMean, y = log2FC, color = sig)) +
  geom_point(size = 0.4, alpha = 0.5) +
  scale_x_log10() +
  scale_color_manual(values = c(DE="#dc2626", NS="#9ca3af")) +
  geom_hline(yintercept = 0, color = "grey40") +
  labs(x = "Mean expression (log10)", y = "log2 Fold Change") +
  theme_classic()

mat_z <- t(scale(t(expr_mat)))
pheatmap(mat_z,
  color = colorRampPalette(c("#2563eb","white","#dc2626"))(100),
  breaks = seq(-3, 3, length.out = 101),
  annotation_col = anno,
  show_colnames = FALSE, fontsize_row = 8)

d <- dist(t(assay(vsd)))
pheatmap(as.matrix(d),
  clustering_distance_rows = d,
  clustering_distance_cols = d,
  annotation_row = anno,
  col = colorRampPalette(rev(brewer.pal(9,"Blues")))(255))

ego <- enrichGO(gene = entrez_ids, OrgDb = org.Hs.eg.db,
                ont = "BP", pvalueCutoff = 0.05, readable = TRUE)
dotplot(ego, showCategory = 15) +
  scale_color_gradient(low = "#dc2626", high = "#2563eb") +
  theme_classic()

gsea_res <- gseGO(geneList = ranked_genes, OrgDb = org.Hs.eg.db,
                   ont = "BP", pvalueCutoff = 0.25)
gseaplot2(gsea_res, geneSetID = top_term_id,
          pvalue_table = TRUE, color = "#dc2626")

ggplot(top_terms, aes(x = reorder(Description, NES), y = NES, fill = direction)) +
  geom_col(width = 0.7) + coord_flip() +
  scale_fill_manual(values = c(Activated="#dc2626", Suppressed="#2563eb")) +
  labs(x = NULL, y = "NES") + theme_classic()

vsd <- vst(dds, blind = FALSE)
pca <- plotPCA(vsd, intgroup = "condition", returnData = TRUE)
pct <- round(100 * attr(pca, "percentVar"))
ggplot(pca, aes(PC1, PC2, color = condition)) +
  geom_point(size = 4) +
  geom_text_repel(aes(label = name), size = 3) +
  labs(x = paste0("PC1 (",pct[1],"%)"), y = paste0("PC2 (",pct[2],"%)")) +
  theme_classic()

DimPlot(seurat_obj, reduction = "umap",
        label = TRUE, label.size = 4, repel = TRUE) +
  labs(title = "UMAP Clustering") +
  theme_classic()

ggplot(df, aes(x = condition, y = expression, fill = condition)) +
  geom_violin(width = 0.8, alpha = 0.7, color = NA) +
  geom_boxplot(width = 0.15, fill = "white", outlier.shape = NA) +
  geom_jitter(width = 0.1, size = 1.5, alpha = 0.6) +
  facet_wrap(~ gene, scales = "free_y") +
  scale_fill_manual(values = c("#475569","#dc2626"), guide = "none") +
  theme_classic()

fit <- survfit(Surv(time, status) ~ group, data = df)
ggsurvplot(fit, data = df,
  pval = TRUE, risk.table = TRUE,
  palette = c("#0f766e", "#dc2626"),
  xlab = "Time (months)", ylab = "Survival Probability")

library(maftools)
laml <- read.maf(maf = "somatic.maf")
oncoplot(laml, top = 10, fontSize = 0.7)

lollipopPlot(maf_obj, gene = "DNMT3A",
             AACol = "Protein_Change",
             showMutationRate = TRUE)

library(ChIPseeker)
peak <- readPeakFile("peaks.narrowPeak")
anno <- annotatePeak(peak, TxDb = txdb, tssRegion = c(-3000, 3000))
plotAnnoPie(anno)

cnetplot(ego, showCategory = 5, foldChange = fc_vec) +
  scale_color_gradient2(low = "#2563eb", mid = "grey90",
                        high = "#dc2626", midpoint = 0)

library(DiffBind)
data(tamoxifen_counts)
tamoxifen <- dba.contrast(tamoxifen, categories = DBA_CONDITION)
tamoxifen <- dba.analyze(tamoxifen)
dba.plotMA(tamoxifen)

ggplot(qc_df, aes(x = reorder(sample, lib_size), y = lib_size / 1e6, fill = group)) +
  geom_col(width = 0.7) +
  geom_hline(yintercept = 10, linetype = "dashed", color = "red") +
  labs(x = NULL, y = "Library size (millions)") +
  coord_flip() + theme_classic()

dds <- DESeq(dds)
plotDispEsts(dds, main = "Dispersion Estimates")

top_genes <- head(rownames(counts(dds, normalized = TRUE)
  [order(-rowMeans(counts(dds, normalized = TRUE))), ]), 20)
plotCounts_df <- reshape2::melt(counts(dds, normalized = TRUE)[top_genes, ])
ggplot(plotCounts_df, aes(x = Var1, y = value, fill = Var1)) +
  geom_boxplot() + coord_flip() +
  labs(x = NULL, y = "Normalized counts") +
  theme_classic() + theme(legend.position = "none")

ego_all <- enrichGO(gene = entrez_ids, OrgDb = org.Hs.eg.db,
                    ont = "ALL", pvalueCutoff = 0.05, readable = TRUE)
barplot(ego_all, showCategory = 5, split = "ONTOLOGY") +
  facet_grid(ONTOLOGY ~ ., scales = "free_y", space = "free_y") +
  theme_classic()

ego <- pairwise_termsim(ego)
emapplot(ego, showCategory = 30, cex_label_category = 0.6) +
  theme_void()

ekegg <- enrichKEGG(gene = entrez_ids, organism = "hsa",
                    pvalueCutoff = 0.05)
dotplot(ekegg, showCategory = 15) +
  scale_color_gradient(low = "#dc2626", high = "#2563eb") +
  theme_classic()

ck <- compareCluster(geneCluster = gene_list, fun = "enrichKEGG",
                      organism = "hsa")
dotplot(ck, showCategory = 10, label_format = 40) +
  scale_color_gradient(low = "#dc2626", high = "#2563eb") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplot(gsea_top, aes(x = NES, y = reorder(Description, NES))) +
  geom_point(aes(size = setSize, color = p.adjust)) +
  geom_errorbarh(aes(xmin = NES - 0.3, xmax = NES + 0.3), height = 0.2) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_color_gradient(low = "#dc2626", high = "#2563eb") +
  labs(x = "NES", y = NULL) + theme_classic()

ggplot(traj_df, aes(x = time, y = expression, group = gene, color = cluster)) +
  geom_line(alpha = 0.4) +
  geom_point(size = 1.5) +
  stat_summary(aes(group = cluster), fun = mean, geom = "line", linewidth = 1.5) +
  facet_wrap(~ cluster) +
  labs(x = "Time point", y = "Scaled expression") +
  theme_classic()

library(Mfuzz)
eset_s <- standardise(eset)
m <- mestimate(eset_s)
cl <- mfuzz(eset_s, c = 6, m = m)
mfuzz.plot2(eset_s, cl = cl, mfrow = c(2, 3),
            time.labels = paste0("T", 1:5), colo = "fancy")

# Before correction
p1 <- plotPCA(vsd_before, intgroup = "batch") + ggtitle("Before ComBat")
# After correction
combat_mat <- ComBat_seq(counts, batch = batch, group = condition)
vsd_after <- vst(DESeqDataSetFromMatrix(combat_mat, colData, ~condition))
p2 <- plotPCA(vsd_after, intgroup = "condition") + ggtitle("After ComBat")
p1 + p2

ggplot(compare_df, aes(x = lfc_deseq2, y = lfc_edger)) +
  geom_point(size = 0.8, alpha = 0.4, color = "#475569") +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +
  geom_smooth(method = "lm", se = FALSE, color = "#2563eb") +
  annotate("text", x = -3, y = 4, label = paste0("r = ", round(cor_val, 3))) +
  labs(x = "log2FC (DESeq2)", y = "log2FC (edgeR)") +
  theme_classic()

cds <- learn_graph(cds)
cds <- order_cells(cds, root_cells = root_id)
plot_cells(cds, color_cells_by = "pseudotime",
           cell_size = 0.8, trajectory_graph_color = "grey30",
           label_cell_groups = FALSE) +
  scale_color_viridis_c() + theme_void()

library(CellChat)
cellchat <- computeCommunProb(cellchat)
cellchat <- computeCommunProbPathway(cellchat)
netVisual_chord_gene(cellchat, signaling = "WNT",
                     color.use = cell_colors, legend.pos.x = 8)

SpatialDimPlot(seurat_spatial, label = TRUE, label.size = 3,
               pt.size.factor = 1.6, image.alpha = 0.6) +
  theme(legend.position = "right")

SpatialFeaturePlot(seurat_spatial, features = c("EPCAM", "VIM"),
                    pt.size.factor = 1.8, alpha = c(0.1, 1),
                    max.cutoff = "q95", image.alpha = 0.5) +
  scale_fill_viridis_c()