```{r}
drug_B1.1=drug_B1[match(colnames(B1.cancer.1), names(drug_B1))]
hist(drug_B1.1,breaks=100)
abline(v = 4000, col = "red", lwd = 2)

```



```{r}
library(ggrepel)

# Add drug_B1.1 to the Seurat object's metadata (if not already added)
B1.cancer.1[["drug_B1.1"]] <- drug_B1.1

# Create a binary grouping variable based on drug_B1.1
B1.cancer.1$group <- ifelse(B1.cancer.1[["drug_B1.1"]] > 4000, "above_4000", "below_4000")

# Convert 'group' to a factor and set it as the active identity class
B1.cancer.1 <- SetIdent(B1.cancer.1, value = as.factor(B1.cancer.1$group))
#B1.cancer.1=B1.cancer.1  %>% UpdateSeuratObject()
#B1.cancer.1=B1.cancer.1  %>% SCTransform(assay = "Spatial", verbose = FALSE)
# Perform differential expression analysis
markers <- FindMarkers(B1.cancer.1, ident.1 = "above_4000", ident.2 = "below_4000",logfc.threshold =0.25)
#markers=markers %>% filter(p_val_adj<0.05)

markers=markers[,-1]
markers %>% print()
top_up <- head(markers[order(markers$p_val_adj, decreasing = FALSE), ], 10)

# Top 5 downregulated
top_down <- head(markers[order(markers$avg_log2FC), ], 5)

# Combine top genes for labeling
top_genes <- rbind(top_up, top_down)
ggplot(markers, aes(x = avg_log2FC, y = -log10(p_val_adj))) +
  geom_point(aes(fill = avg_log2FC), color = "black", shape = 21, alpha = 0.5, stroke = 0.5) + 
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0) +  # Adjust for fill
  geom_text_repel(data = top_up, 
                  aes(label = rownames(top_up), x = avg_log2FC, y = -log10(p_val_adj)),
                  box.padding = 0.35, 
                  point.padding = 0.5,
                  size = 3) +
  labs(x = "Log2 Fold Change", y = "-Log10 P-value") +
  theme_minimal() +
  ggtitle("Volcano Plot of Differentially Expressed Genes")




markers
```
```{r}
SpatialFeaturePlot(B1.cancer.1,c("Inmt","drug_conc"))
SpatialFeaturePlot(postMALDI.list$B1,"Cd38")

```


```{r}

# Add drug_B1.1 to the Seurat object's metadata (if not already added)
B1.cancer.1[["drug_B1.1"]] <- drug_B1.1

# Create a binary grouping variable based on drug_B1.1
B1.cancer.1$group <- ifelse(B1.cancer.1[["drug_B1.1"]] > 4000, "above_4000", "below_4000")

# Convert 'group' to a factor and set it as the active identity class
B1.cancer.1 <- SetIdent(B1.cancer.1, value = as.factor(B1.cancer.1$group))
#B1.cancer.1=B1.cancer.1  %>% UpdateSeuratObject()
#B1.cancer.1=B1.cancer.1  %>% SCTransform(assay = "Spatial", verbose = FALSE)
# Perform differential expression analysis
markers <- FindMarkers(B1.cancer.1, ident.1 = "above_4000", ident.2 = "below_4000",logfc.threshold = 0.5,test.use = "t")
markers=markers %>% filter(p_val_adj<0.05)
options(scipen = 5)

```

```{r}
library(fgsea)
library(data.table)

# Assuming you have a dataframe 'genes' with columns 'geneSymbol', 'foldChange', and 'pValue'
markers$rank <- -log10(markers$p_val_adj) * sign(markers$avg_log2FC)
rankedGenes <- markers[order(markers$rank, decreasing = TRUE),]
# If you have gene sets in GMT format
geneSets <- gmtPathways("/gpfs/data/yanailab/public/jingye_share/20231129/m5.go.v2023.2.Mm.entrez.gmt.txt")
a= rankedGenes$rank
names=mapIds(org.Mm.eg.db,
                    keys = rownames(markers),
                    column = "ENTREZID",
                    keytype = "SYMBOL",
                    multiVals = "first")
a=setNames(a,names)
str(a)
# Ensure the names of the ranked list match the gene sets
fgseaRes <- fgsea(pathways = geneSets, 
                  stats = a, 
                  minSize = 15, # Minimum size of gene set to consider
                  maxSize = 500 )# Maximum size of gene set to consider) # Number of permutations for estimating p-values
head(fgseaRes)
```




```{r}
library(clusterProfiler)
#gene_list <-markers$avg_log2FC %>% sign()*markers$p_val_adj
gene_list <-markers$avg_log2FC
names(gene_list) <- rownames(markers)

# Ensure the gene_list is in decreasing order of association
gene_list <- sort(gene_list, decreasing = TRUE)

library(org.Mm.eg.db)

# Run GSEA
gsea_results <- gseGO(geneList = gene_list, 
                      ont = "BP", 
                      keyType = "SYMBOL", 
                      OrgDb = org.Mm.eg.db, 
                      minGSSize = 10,
                      maxGSSize = 500,
                      pAdjustMethod = "BH",
                      pvalueCutoff = 0.1)

# View results
gsea_results@result[["Description"]] %>% print()
gsea_results@result[["core_enrichment"]] %>% print()
#SpatialFeaturePlot(B1.cancer.1, features = c("Rgs9","drug_conc"))

```


```{r}

library(enrichR)



# perform_enrichment function as previously defined
perform_enrichment <- function(gene_list, db, title_suffix) {


  # Perform Enrichment Analysis
  enrich_results <- enrichR::enrichr(names(gene_list), db)

  # Check if enrich_results is not empty
  if (length(enrich_results) > 0 && nrow(enrich_results[[1]]) > 0) {
    # Select the first result (or a specific one if you prefer)
    result <- enrich_results[[1]]


    # Plotting the top 10 results
    ggplot(head(result, 10), aes(x = reorder(Term, -Adjusted.P.value), y = -log10(Adjusted.P.value))) +
      geom_bar(stat = "identity", fill = "steelblue") +
      theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
      labs(x = "Pathway", y = "-log10(Adjusted P-value)", title = paste("Top 10 Pathways from", title_suffix))
  } else {
    message("No enrichment results to plot or empty results for", title_suffix, ".")
  }
}


        # Proceed with enrichment analysis if markers are found
        if (nrow(markers) > 0) {
            gene_list <- markers$avg_log2FC
            names(gene_list) <- rownames(markers)
            gene_list <- na.omit(gene_list)
            gene_list <- sort(gene_list[gene_list > 0], decreasing = TRUE)

            # Perform enrichment analysis
            db <- "KEGG_2019_Mouse"  # Adjust the database as needed
            title_suffix <- paste0("Subset ", x, "1 - ", ident)
            perform_enrichment(gene_list, db, title_suffix)
        }
        if (nrow(markers) > 0) {
            gene_list <- markers$avg_log2FC
            names(gene_list) <- rownames(markers)
            gene_list <- na.omit(gene_list)
            gene_list <- sort(gene_list[gene_list < 0], decreasing = TRUE)

            # Perform enrichment analysis
            db <- "KEGG_2019_Mouse"  # Adjust the database as needed
            title_suffix <- paste0("Subset ", x, "1 - ", ident)
            perform_enrichment(gene_list, db, title_suffix)
}

gene_list %>% print()
```


```{r}
library(NMF)
DefaultAssay(postMALDI.list$B1) <- "Spatial"

# Normalize, find variable features, and scale
postMALDI.list$B1 <- postMALDI.list$B1 %>% 
                     NormalizeData() %>% 
                     FindVariableFeatures(nfeatures=3000) %>% 
                     ScaleData()
topgenes=head(VariableFeatures(postMALDI.list$B1), 3000)
nmf_expr<-postMALDI.list$B1@assays$Spatial@scale.data # raw = @counts, normalized = @data
#nmf_expr<-pmax(nmf_expr,0)
nmf_expr <- as.matrix(nmf_expr)
nmf_expr<-nneg(nmf_expr,method = "pmax")
nmf_expr_est2<-nmf_expr
nmf_survey<-nmf(nmf_expr_est2, rank=c(6,10,12,15), method='nsNMF', seed=123456,.options='p+v')
nmf_10<-nmf(nmf_expr_est2, rank=10, method='nsNMF', seed=123456,.options='p+v')
print(plot(nmf_10))

```
```{r}
library(NMF)
DefaultAssay(postMALDI.list$A1) <- "Spatial"

# Normalize, find variable features, and scale
postMALDI.list$A1 <- postMALDI.list$A1 %>% 
                     NormalizeData() %>% 
                     FindVariableFeatures(nfeatures=3000) %>% 
                     ScaleData()
topgenes=head(VariableFeatures(postMALDI.list$A1), 3000)
nmf_expr<-postMALDI.list$A1@assays$Spatial@scale.data # raw = @counts, normalized = @data
#nmf_expr<-pmax(nmf_expr,0)
nmf_expr <- as.matrix(nmf_expr)
nmf_expr<-nneg(nmf_expr,method = "pmax")
nmf_expr_est2<-nmf_expr
#nmf_survey<-nmf(nmf_expr_est2, rank=c(6,10,12,15), method='nsNMF', seed=123456,.options='p+v')
nmf_12_A1<-nmf(nmf_expr_est2, rank=12, method='nsNMF', seed=123456,.options='p+v')

DefaultAssay(postMALDI.list$B1) <- "Spatial"

# Normalize, find variable features, and scale
postMALDI.list$B1 <- postMALDI.list$B1 %>% 
                     NormalizeData() %>% 
                     FindVariableFeatures(nfeatures=3000) %>% 
                     ScaleData()
topgenes=head(VariableFeatures(postMALDI.list$B1), 3000)
nmf_expr<-postMALDI.list$B1@assays$Spatial@scale.data # raw = @counts, normalized = @data
#nmf_expr<-pmax(nmf_expr,0)
nmf_expr <- as.matrix(nmf_expr)
nmf_expr<-nneg(nmf_expr,method = "pmax")
nmf_expr_est2<-nmf_expr
#nmf_survey<-nmf(nmf_expr_est2, rank=c(6,10,12,15), method='nsNMF', seed=123456,.options='p+v')
nmf_12_B1<-nmf(nmf_expr_est2, rank=12, method='nsNMF', seed=123456,.options='p+v')

DefaultAssay(postMALDI.list$C1) <- "Spatial"

# Normalize, find variable features, and scale
postMALDI.list$C1 <- postMALDI.list$C1 %>% 
                     NormalizeData() %>% 
                     FindVariableFeatures(nfeatures=3000) %>% 
                     ScaleData()
topgenes=head(VariableFeatures(postMALDI.list$C1), 3000)
nmf_expr<-postMALDI.list$C1@assays$Spatial@scale.data # raw = @counts, normalized = @data
#nmf_expr<-pmax(nmf_expr,0)
nmf_expr <- as.matrix(nmf_expr)
nmf_expr<-nneg(nmf_expr,method = "pmax")
nmf_expr_est2<-nmf_expr

nmf_12_C1<-nmf(nmf_expr_est2, rank=12, method='nsNMF', seed=123456,.options='p+v')
print(plot(nmf_10))

DefaultAssay(postMALDI.list$C1) <- "Spatial"

# Normalize, find variable features, and scale
postMALDI.list$C1 <- postMALDI.list$C1 %>% 
                     NormalizeData() %>% 
                     FindVariableFeatures(nfeatures=3000) %>% 
                     ScaleData()
topgenes=head(VariableFeatures(postMALDI.list$C1), 3000)
nmf_expr<-postMALDI.list$C1@assays$Spatial@scale.data # raw = @counts, normalized = @data
#nmf_expr<-pmax(nmf_expr,0)
nmf_expr <- as.matrix(nmf_expr)
nmf_expr<-nneg(nmf_expr,method = "pmax")
nmf_expr_est2<-nmf_expr

nmf_12_C1<-nmf(nmf_expr_est2, rank=12, method='nsNMF', seed=123456,.options='p+v')


DefaultAssay(postMALDI.list$D1) <- "Spatial"

# Normalize, find variable features, and scale
postMALDI.list$D1 <- postMALDI.list$D1 %>% 
                     NormalizeData() %>% 
                     FindVariableFeatures(nfeatures=3000) %>% 
                     ScaleData()
topgenes=head(VariableFeatures(postMALDI.list$D1), 3000)
nmf_expr<-postMALDI.list$D1@assays$Spatial@scale.data # raw = @counts, normalized = @data
#nmf_expr<-pmax(nmf_expr,0)
nmf_expr <- as.matrix(nmf_expr)
nmf_expr<-nneg(nmf_expr,method = "pmax")
nmf_expr_est2<-nmf_expr

nmf_12_D1<-nmf(nmf_expr_est2, rank=12, method='nsNMF', seed=123456,.options='p+v')

```



```{r}
#save.image("/gpfs/data/yanailab/public/jingye_share/20231129/correlation.RData")
# Extracting NMF factor scores (assuming nmf_10@fit@H is factor by cells)
factor_scores <- nmf_12_B1@fit@H

# Transpose matrix if necessary to ensure it's cells by factors
if (nrow(factor_scores) != ncol(postMALDI.list$B1)) {
  factor_scores <- t(factor_scores)
}



# Assuming factor_scores is a matrix with columns corresponding to cells
for (i in 1:ncol(factor_scores)) {
    # Create a name for the factor
    factor_name <- paste("Factor", i, sep = "_")
    
    # Add factor scores as metadata to the Seurat object
    postMALDI.list$B1 <- AddMetaData(postMALDI.list$B1, metadata = factor_scores[, i], col.name = factor_name)
}

# Now, you can create spatial feature plots for each factor
for (i in 1:ncol(factor_scores)) {
    factor_name <- paste("Factor", i, sep = "_")
    print(SpatialFeaturePlot(postMALDI.list$B1, features = factor_name) + 
        theme_minimal() +
        ggtitle(paste("Spatial Plot for", factor_name)))
}
SpatialFeaturePlot(postMALDI.list$B1, features = "drug_conc") %>% print()
SpatialFeaturePlot(B1.cancer.1, features = "nUMI_Spatial") %>% print()

SpatialFeaturePlot(B1.cancer.1, features = "drug_conc") %>% print()
SpatialFeaturePlot(B1.cancer.7, features = "drug_conc") %>% print()

```

```{r}
library(Seurat)

# Extract the drug_conc metadata
drug_conc <- postMALDI.list$B1[["drug_conc"]]

# Initialize a vector to store correlation results
correlations <- numeric(ncol(factor_scores))

# Compute correlations for each factor
for (i in 1:ncol(factor_scores)) {
    factor_name <- paste("Factor", i, sep = "_")
    factor_scores <- postMALDI.list$B1[[factor_name]]

    # Compute correlation and store the result
    correlations[i] <- cor(drug_conc, factor_scores, use = "complete.obs")
}


# Print the correlation results
print(correlations)
VlnPlot(postMALDI.list$B1, features = "Trp53",pt.size = 0.1) +
  theme_minimal()
```

```{r}
# Extract the 5th column (factor 5 loadings) from the W matrix
factor_5_loadings <- nmf_10@fit@W[, 5]

# Sort the loadings to find the top ones
# Adjust 'n' for the number of top loadings you want to retrieve
top_genes_indices <- order(factor_5_loadings, decreasing = TRUE)[1:1000]

# Assuming the row names of W are gene names
top_genes <- rownames(nmf_10@fit@W)[top_genes_indices]

# Combine gene names with their loadings
top_genes_loadings <- data.frame(Gene = top_genes, Loading = factor_5_loadings[top_genes_indices])

# Print the top genes and their loadings
print(top_genes_loadings)
top_genes_loadings$Gene

```
```{r}
library(clusterProfiler)
library(org.Mm.eg.db)

# Assuming nmf_10 is your NMF result
# Extract the 5th column (factor 5 loadings) from the W matrix
factor_5_loadings <- nmf_10@fit@W[, 5]

# Get top 1000 genes
top_genes_indices <- order(factor_5_loadings, decreasing = TRUE)[1:1000]
top_genes <- rownames(nmf_10@fit@W)[top_genes_indices]

# Create a named vector for GSEA, using loadings as the values
gene_list <- factor_5_loadings[top_genes_indices]
names(gene_list) <- top_genes

# Sort the gene list in decreasing order of loadings
gene_list <- sort(gene_list, decreasing = TRUE)

# Run GSEA
gsea_results <- gseGO(geneList = gene_list, 
                      ont = "BP", 
                      keyType = "SYMBOL", 
                      OrgDb = org.Mm.eg.db, 
                      minGSSize = 10,
                      maxGSSize = 500,
                      pAdjustMethod = "BH",
                      pvalueCutoff = 0.05)

# View results
print(gsea_results@result[["Description"]])
print(gsea_results@result[["core_enrichment"]])

```


```{r}
SpatialFeaturePlot(postMALDI.list$B1, features = "Epcam")
current_idents <- Idents(postMALDI.list$B1)

# Create a violin plot for the Epcam gene using the current identities
VlnPlot(postMALDI.list$B1, features = "Epcam",pt.size = 0.1) +
  theme_minimal()
VlnPlot(postMALDI.list$B1, features = "nFeature_Spatial",pt.size = 0.1) +
  theme_minimal()

```
```{r}
SpatialDimPlot(postMALDI.list$B1)
```

```{r}
SpatialFeaturePlot(postMALDI.list$B1, features = "nFeature_Spatial") %>% print()
SpatialFeaturePlot(B1.cancer.1, features = "nFeature_Spatial") %>% print()
SpatialFeaturePlot(postMALDI.list$B1, features = c("Cyp2f2","drug_conc")) %>% print()

```

C1
C1
C1
C1
C1
C1
```{r}
drug_C1.2=drug_C1[match(colnames(C1.cancer.2), names(drug_C1))]
hist(drug_C1.2,breaks=100)
abline(v = 4500, col = "red", lwd = 2)
```

```{r}

# Add drug_C1.2 to the Seurat object's metadata (if not already added)
C1.cancer.2[["drug_C1.2"]] <- drug_C1.2

# Create a binary grouping variable based on drug_C1.2
C1.cancer.2$group <- ifelse(C1.cancer.2[["drug_C1.2"]] > 4500, "above_4500", "below_4500")

# Convert 'group' to a factor and set it as the active identity class
C1.cancer.2 <- SetIdent(C1.cancer.2, value = as.factor(C1.cancer.2$group))
#C1.cancer.2=C1.cancer.2  %>% UpdateSeuratObject()
#C1.cancer.2=C1.cancer.2  %>% SCTransform(assay = "Spatial", verbose = FALSE)
# Perform differential expression analysis
markers <- FindMarkers(C1.cancer.2, ident.1 = "above_4500", ident.2 = "below_4500",logfc.threshold = 0.25)
#markers=markers %>% filter(p_val_adj<0.05)

markers
markers=markers[,-1]
markers %>% print()
top_up <- head(markers[order(markers$p_val_adj, decreasing = FALSE), ], 10)

# Top 5 downregulated
top_down <- head(markers[order(markers$avg_log2FC), ], 5)

# Combine top genes for labeling
top_genes <- rbind(top_up, top_down)

ggplot(markers, aes(x = avg_log2FC, y = -log10(p_val_adj))) +
  geom_point(aes(fill = avg_log2FC), color = "black", shape = 21, alpha = 0.5, stroke = 0.5) + 
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0) +  # Adjust for fill
  geom_text_repel(data = top_up, 
                  aes(label = rownames(top_up), x = avg_log2FC, y = -log10(p_val_adj)),
                  box.padding = 0.35, 
                  point.padding = 0.5,
                  size = 3) +
  labs(x = "Log2 Fold Change", y = "-Log10 P-value") +
  theme_minimal() +
  ggtitle("Volcano Plot of Differentially Expressed Genes")

markers

```


```{r}
library(fgsea)
library(data.table)
markers=markers %>% filter(p_val_adj<0.05)
# Assuming you have a dataframe 'genes' with columns 'geneSymbol', 'foldChange', and 'pValue'
markers$rank <- -log10(markers$p_val_adj) * sign(markers$avg_log2FC)
rankedGenes <- markers[order(markers$rank, decreasing = TRUE),]
# If you have gene sets in GMT format
geneSets <- gmtPathways("/gpfs/data/yanailab/public/jingye_share/20231129/m5.go.v2023.2.Mm.entrez.gmt.txt")
a= rankedGenes$rank
names=mapIds(org.Mm.eg.db,
                    keys = rownames(markers),
                    column = "ENTREZID",
                    keytype = "SYMBOL",
                    multiVals = "first")
a=setNames(a,names)
str(a)
# Ensure the names of the ranked list match the gene sets
fgseaRes <- fgsea(pathways = geneSets, 
                  stats = a, 
                  minSize = 15, # Minimum size of gene set to consider
                  maxSize = 500 )# Maximum size of gene set to consider) # Number of permutations for estimating p-values
head(fgseaRes)
```

```{r}


library(enrichR)



# perform_enrichment function as previously defined
perform_enrichment <- function(gene_list, db, title_suffix) {


  # Perform Enrichment Analysis
  enrich_results <- enrichR::enrichr(names(gene_list), db)

  # Check if enrich_results is not empty
  if (length(enrich_results) > 0 && nrow(enrich_results[[1]]) > 0) {
    # Select the first result (or a specific one if you prefer)
    result <- enrich_results[[1]]


    # Plotting the top 10 results
    ggplot(head(result, 10), aes(x = reorder(Term, -Adjusted.P.value), y = -log10(Adjusted.P.value))) +
      geom_bar(stat = "identity", fill = "steelblue") +
      theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
      labs(x = "Pathway", y = "-log10(Adjusted P-value)", title = paste("Top 10 Pathways from", title_suffix))
  } else {
    message("No enrichment results to plot or empty results for", title_suffix, ".")
  }
}


        # Proceed with enrichment analysis if markers are found
        if (nrow(markers) > 0) {
            gene_list <- markers$avg_log2FC
            names(gene_list) <- rownames(markers)
            gene_list <- na.omit(gene_list)
            gene_list <- sort(gene_list[gene_list > 0], decreasing = TRUE)

            # Perform enrichment analysis
            db <- "KEGG_2019_Mouse"  # Adjust the database as needed
            title_suffix <- paste0("Subset ", x, "1 - ", ident)
            perform_enrichment(gene_list, db, title_suffix)
        }
        if (nrow(markers) > 0) {
            gene_list <- markers$avg_log2FC
            names(gene_list) <- rownames(markers)
            gene_list <- na.omit(gene_list)
            gene_list <- sort(gene_list[gene_list < 0], decreasing = TRUE)

            # Perform enrichment analysis
            db <- "KEGG_2019_Mouse"  # Adjust the database as needed
            title_suffix <- paste0("Subset ", x, "1 - ", ident)
            perform_enrichment(gene_list, db, title_suffix)
}

gene_list %>% print()

```

```{r}
drug_D1.2=drug_D1[match(colnames(D1.cancer.2), names(drug_D1))]
hist(drug_D1.2,breaks=100)
```

```{r}

# Add drug_D1.2 to the Seurat object's metadata (if not already added)
D1.cancer.2[["drug_D1.2"]] <- drug_D1.2

# Create a binary grouping variable based on drug_D1.2
D1.cancer.2$group <- ifelse(D1.cancer.2[["drug_D1.2"]] > 5000, "above_5000", "below_5000")

# Convert 'group' to a factor and set it as the active identity class
D1.cancer.2 <- SetIdent(D1.cancer.2, value = as.factor(D1.cancer.2$group))
D1.cancer.2=D1.cancer.2  %>% UpdateSeuratObject()
#D1.cancer.2=D1.cancer.2  %>% SCTransform(assay = "Spatial", verbose = FALSE)
# Perform differential expression analysis
markers <- FindMarkers(D1.cancer.2, ident.1 = "above_5000", ident.2 = "below_5000",logfc.threshold = 0.5)
markers=markers %>% filter(p_val_adj<0.05)

markers


```
```{r}
SpatialFeaturePlot(postMALDI.list$D1, features = c("Cyp2f2","drug_conc")) %>% print()
SpatialFeaturePlot(postMALDI.list$D1, features = "nFeature_Spatial") %>% print()
```

```{r}
library(NMF)
library(proxy)
library(ComplexHeatmap)

# Step 1: Extract and threshold gene sets
get_gene_sets <- function(nmf_object, threshold = 0.001) {
  apply(nmf_object@fit@W, 2, function(x) names(x[x > threshold]))
}

gene_sets_A1 <- get_gene_sets(nmf_12_A1)
gene_sets_B1 <- get_gene_sets(nmf_12_B1)
gene_sets_C1 <- get_gene_sets(nmf_12_C1)

# Combine all gene sets into a list
all_gene_sets <- c(gene_sets_A1, gene_sets_B1, gene_sets_C1)

# Step 2: Calculate Jaccard similarity
jaccard_similarity <- function(set1, set2) {
  if(length(set1) == 0 & length(set2) == 0) {
    return(1) # If both sets are empty, return similarity of 1
  }
  intersection_size <- length(intersect(set1, set2))
  union_size <- length(union(set1, set2))
  similarity <- intersection_size / union_size
  return(similarity)
}

# Initialize an empty matrix for Jaccard similarity
num_factors <- length(all_gene_sets)
similarity_matrix <- matrix(0, nrow = num_factors, ncol = num_factors,
                            dimnames = list(names(all_gene_sets), names(all_gene_sets)))

# Calculate Jaccard similarity for every pair of gene sets
for(i in 1:num_factors) {
  for(j in i:num_factors) {
    similarity <- jaccard_similarity(all_gene_sets[[i]], all_gene_sets[[j]])
    similarity_matrix[i, j] <- similarity
    similarity_matrix[j, i] <- similarity # Symmetric
  }
}

# Step 3: Prepare data for heatmap and sample identity
factor_names <- c(paste0("A1_Factor_", seq_along(gene_sets_A1)),
                  paste0("B1_Factor_", seq_along(gene_sets_B1)),
                  paste0("C1_Factor_", seq_along(gene_sets_C1)))

sample_identity <- c(rep("A1", length(gene_sets_A1)),
                     rep("B1", length(gene_sets_B1)),
                     rep("C1", length(gene_sets_C1)))

# Create a color mapping for sample identity
sample_colors <- setNames(c("red", "blue", "green"), c("A1", "B1", "C1"))

# Step 4: Plot the heatmap
library(ComplexHeatmap)
library(circlize) # colorRamp2 is available in the circlize package, which ComplexHeatmap depends on

# Define breakpoints and colors
breaks <- c(0, 0.25, 0.5, 1)
colors <- c("blue", "white", "red", "red") # Adjust colors as needed

# Create color mapping function
color_mapping <- colorRamp2(breaks, colors)

# Plot the heatmap with the custom color scale
Heatmap(similarity_matrix, name = "Jaccard Similarity",
        row_labels = factor_names, column_labels = factor_names,
        row_title = "Factors", column_title = "Factors",
        top_annotation = HeatmapAnnotation(Sample = sample_identity, col = list(Sample = sample_colors)),
        col = color_mapping) # Use the custom color mapping


```

```{r}
factor_names
```

