Signature Markers(Malignat_vs_Normal_CD4Tcells) filtered on MeanExpression

1. load libraries

2. Load the filtered upregulated and Downregulated genes

3. Initial test

df <- df[!(df$Malignant.exp < 0.20 & df$Control.exp < 0.20), ]

df <- as.data.frame(df)

write.csv(df, "../18-03-25_top-markers/Psedibulk_Deseq2_filtered_on_mean.csv", row.names = FALSE)

# Ensure the dataframe has necessary columns: "gene", "avg_logFC", "p_val_adj"
# If gene names are stored as row names, move them to a column
if (!"gene" %in% colnames(df)) {
  df <- rownames_to_column(df, "gene")
}

4. Initial test and visualization

library(dplyr)
library(ggplot2)

# Filter significant genes (p_val_adj < 0.05)
significant_genes <- df %>%
  filter(p_val_adj < 1e-4)

# Select top 30 upregulated genes (highest avg_logFC)
top_30_up <- significant_genes %>%
  arrange(desc(avg_logFC), p_val_adj) %>%
  head(30)

# Filter out genes starting with "TRVB" and "TRAV" and select top 30 downregulated genes (lowest avg_logFC)
top_30_down <- significant_genes %>%
  filter(!grepl("^TRBV", gene) & !grepl("^TRAV", gene)) %>%  # Exclude genes starting with "TRVB" or "TRAV"
  arrange(avg_logFC, p_val_adj) %>%
  head(30)

# Combine both into one dataframe with a new column indicating direction
top_genes <- bind_rows(
  top_30_up %>% mutate(direction = "Upregulated"),
  top_30_down %>% mutate(direction = "Downregulated")
)

# Save both separately and combined
write.csv(top_30_up, "top_30_upregulated.csv", row.names = FALSE)
write.csv(top_30_down, "top_30_downregulated.csv", row.names = FALSE)
write.csv(top_genes, "top_30_biomarkers_combined.csv", row.names = FALSE)

# Print selected columns
print(top_genes %>% select(gene, avg_logFC, p_val_adj, Malignant.pct, Control.pct, Malignant.exp, Control.exp, direction))

# For the combined plot of upregulated and downregulated genes with custom colors
ggplot(top_genes, aes(x = reorder(gene, avg_logFC), y = avg_logFC, fill = avg_logFC > 0)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  scale_fill_manual(values = c("blue", "red")) +  # Red for upregulated, Blue for downregulated
  theme_minimal() +
  labs(title = "Top 30 Upregulated and Downregulated Genes",
       x = "Gene",
       y = "Log2 Fold Change",
       fill = "Expression")

# For separate plots for upregulated and downregulated genes
ggplot(top_30_up, aes(x = reorder(gene, avg_logFC), y = avg_logFC, fill = avg_logFC > 0)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  scale_fill_manual(values = c("red", "blue")) +
  theme_minimal() +
  labs(title = "Top 30 Upregulated Genes",
       x = "Gene",
       y = "Log2 Fold Change",
       fill = "Expression")

ggplot(top_30_down, aes(x = reorder(gene, avg_logFC), y = avg_logFC, fill = avg_logFC < 0)) +  # < 0 for downregulated genes
  geom_bar(stat = "identity") +
  coord_flip() +
  scale_fill_manual(values = c("blue", "red")) +  # Red for upregulated, Blue for downregulated
  theme_minimal() +
  labs(title = "Top 30 Downregulated Genes",
       x = "Gene",
       y = "Log2 Fold Change",
       fill = "Expression")

NA
NA

5. Volcano Plot

library(ggplot2)
library(ggrepel)

# Volcano plot with labeled top genes
volcano_plot <- ggplot(top_genes, aes(x = avg_logFC, y = -log10(p_val_adj))) +
  geom_point(aes(color = avg_logFC), size = 3, alpha = 0.7) +
  scale_color_gradient2(
    low = "blue", 
    mid = "grey", 
    high = "red", 
    midpoint = median(top_genes$avg_logFC)
  ) +
  geom_text_repel(
    aes(label = gene),  # Label genes
    size = 3, 
    max.overlaps = 20  # Adjust to avoid label crowding
  ) +
  theme_minimal() +
  labs(
    title = "Volcano Plot of Top 30 Biomarkers",
    x = "log2 Fold Change (avg_logFC)",
    y = "-log10(Adjusted p-value)",
    color = "Fold Change"
  )

print(volcano_plot)

6. Bar Plot of Top Genes by Fold Change

# Bar plot (sorted by fold change)
bar_plot <- ggplot(top_genes, 
       aes(x = reorder(gene, avg_logFC), y = avg_logFC, 
           fill = avg_logFC)) +
  geom_col() +
  scale_fill_gradient2(
    low = "blue", 
    mid = "white", 
    high = "red", 
    midpoint = median(top_genes$avg_logFC)
  ) +
  coord_flip() +  # Horizontal bars
  theme_minimal() +
  labs(
    title = "Top 30 Biomarkers by log2 Fold Change",
    x = "Gene",
    y = "avg_logFC"
  )

print(bar_plot)

7. dotplot of Expression Prevalence

# Dot plot comparing pct.1 vs pct.2
dot_plot <- ggplot(top_genes, aes(x = Malignant.pct, y = Control.pct, 
                  size = avg_logFC, color = avg_logFC)) +
  geom_point(alpha = 0.8) +
  geom_text_repel(aes(label = gene), size = 3, max.overlaps = 20) +
  scale_color_gradient2(
    low = "blue", 
    mid = "grey", 
    high = "red", 
    midpoint = median(top_genes$avg_logFC)
  ) +
  theme_minimal() +
  labs(
    title = "Expression Prevalence: Target Group (pct.1) vs Control (pct.2)",
    x = "Expression in Target Group (pct.1)",
    y = "Expression in Control (pct.2)",
    color = "Fold Change",
    size = "Fold Change"
  )

print(dot_plot)

8. dotplot of Mean Expression Prevalence

# Dot plot comparing pct.1 vs pct.2
dot_plot <- ggplot(top_genes, aes(x = Malignant.exp, y = Control.exp, 
                  size = avg_logFC, color = avg_logFC)) +
  geom_point(alpha = 0.8) +
  geom_text_repel(aes(label = gene), size = 3, max.overlaps = 20) +
  scale_color_gradient2(
    low = "blue", 
    mid = "grey", 
    high = "red", 
    midpoint = median(top_genes$avg_logFC)
  ) +
  theme_minimal() +
  labs(
    title = "Expression Prevalence: Target Group (exp.1) vs Control (exp.2)",
    x = "Expression in Target Group (exp.1)",
    y = "Expression in Control (exp.2)",
    color = "Fold Change",
    size = "Fold Change"
  )

print(dot_plot)

9. FeaturePlot for Top30 UP

FeaturePlot(All_samples_Merged, 
             features = top_30_up$gene[1:10], 
             reduction = "umap", 
             cols = c("lightblue", "red"),  # Custom color gradient from light blue to red
             label = TRUE)

FeaturePlot(All_samples_Merged, 
             features = top_30_up$gene[11:20], 
             reduction = "umap", 
             cols = c("lightblue", "red"),  # Custom color gradient from light blue to red
             label = TRUE)

FeaturePlot(All_samples_Merged, 
             features = top_30_up$gene[21:30], 
             reduction = "umap", 
             cols = c("lightblue", "red"),  # Custom color gradient from light blue to red
             label = TRUE)

NA
NA

10. FeaturePlot for Top30 DOWN

FeaturePlot(All_samples_Merged, 
             features = top_30_down$gene[1:10], 
             reduction = "umap", 
             cols = c("lightblue", "red"),  # Custom color gradient from light blue to red
             label = TRUE)

FeaturePlot(All_samples_Merged, 
             features = top_30_down$gene[11:20], 
             reduction = "umap", 
             cols = c("lightblue", "red"),  # Custom color gradient from light blue to red
             label = TRUE)

FeaturePlot(All_samples_Merged, 
             features = top_30_down$gene[21:30], 
             reduction = "umap", 
             cols = c("lightblue", "red"),  # Custom color gradient from light blue to red
             label = TRUE)