```{r setup, include=FALSE}
# note that in rmd setwd is not persistent through chunks in .Rmd.
# https://bookdown.org/yihui/rmarkdown-cookbook/working-directory.html
# instead use knitr::opts_knit$set(root.dir = '/tmp') to set a persistent working directory
# NOTE that this doesn't modify the console wd!
#knitr::opts_knit$set(root.dir = '/gpfs/data/yanailab/projects/vt2217/vinh_202308_jingye_ST/r_scripts/')

library(Seurat)
#library(SeuratData)
library(ggplot2)
library(patchwork)
# https://patchwork.data-imaginist.com/articles/guides/annotation.html
library(dplyr)
library(tidyr)
library(future)
# BiocManager::install("clusterProfiler")
# BiocManager::install("EnhancedVolcano")
# BiocManager::install("org.Mm.eg.db")
# BiocManager::install("MAST")
library(clusterProfiler)
#library(org.Mm.eg.db)
#library(msigdbr)
#library(cowplot)
library(harmony)
#library(MAST)
#library(glmGamPoi)
#library(raster)
library(RColorBrewer)
# modify the statement below to match the number of cores in my r session
plan("multisession", workers = 32)
options(future.globals.maxSize = 24000 * 1024^2)  # Set the limit to 16 GB, increases the maximum size for https://nbisweden.github.io/workshop-scRNAseq/labs/compiled/seurat/seurat_07_spatial.html#Integration
# globals exported for the future package.

#make sure you have at least 256 gb ram to run this whole code safely!

set.seed(100)

scripts_path = '/gpfs/data/yanailab/public/jingye_share/' # goes to the project directory
data_path = '/gpfs/data/yanailab/projects/jw7078/processed_data/20231129/' # goes to raw data
save_path = '/gpfs/data/yanailab/public/jingye_share/' # go to new folder or something that you've already made

knitr::opts_chunk$set(echo = TRUE)
```

```{r}
postMALDI.list <- list() # stores seurat objects for 359
# loading in 359- samples (post-MALDI)
for (x in c("A","B","C","D")){
  postMALDI.list[[paste0(x,'1')]]<-Load10X_Spatial(data.dir = paste0(data_path, '20231129-',x,"1",'/outs'))
  }
```

```{r}
for (x in c("A","B","C","D")){
  postMALDI.list[[paste0(x,'1')]]<-PercentageFeatureSet(postMALDI.list[[paste0(x,'1')]],
                                                         "^mt-",col.name="pct.mito",assay = "Spatial")
  postMALDI.list[[paste0(x,'1')]] <-
    postMALDI.list[[paste0(x,'1')]][,postMALDI.list[[paste0(x,'1')]]$nFeature_Spatial>300 
                                                   & postMALDI.list[[paste0(x,'1')]]$nCount_Spatial>200]
  postMALDI.list[[paste0(x,'1')]]$orig.ident<-paste0(x,'1')
  print(VlnPlot(postMALDI.list[[paste0(x,'1')]], 
                features = c("nFeature_Spatial", "nCount_Spatial", "pct.mito"), ncol = 3))
  print(SpatialFeaturePlot(postMALDI.list[[paste0(x,'1')]], 
                           features = c("nFeature_Spatial", "nCount_Spatial", "pct.mito"), ncol = 3))
}

```


```{r}
# wrapper function for normalization, reducing dimensions, finding neighbors and clustering.
# recommended that values for umap_dim, neighbor_dim, and k.param match.
NormSTNeighClust<-function(srt_obj, npcs=30, umap_dim=30, neighbor_dim=30, k.param = 30, clust.res = 0.8){
  srt_obj <- SCTransform(srt_obj, variable.features.n = 3000, assay="Spatial", verbose = F, vst.flavor="v2") 
  srt_obj <- RunPCA(srt_obj, assay ="SCT", npcs=npcs, verbose = F)
  srt_obj <- FindNeighbors(srt_obj, reduction = "pca", dims = 1:neighbor_dim, k.param=k.param, verbose = F) 
  srt_obj <- FindClusters(srt_obj, verbose = F, resolution = clust.res)
  srt_obj<-RunUMAP(srt_obj, dims=1:umap_dim, verbose = F)
  return(srt_obj)
}
for (x in c('A', 'B', 'C', 'D')) {
  postMALDI.list[[paste0(x,'1')]] <- NormSTNeighClust(postMALDI.list[[paste0(x,'1')]])
}

# rename cells to ensure that we can identify origin when merging noMALDI and postMALDI data.
for (x in c('A', 'B', 'C', 'D')) {
  postMALDI.list[[paste0(x,'1')]] <- RenameCells(postMALDI.list[[paste0(x,'1')]], 
                                                      add.cell.id = paste0(x,'1'))
}

```



```{r}
# Code for changing cluster resolutions (this was copied up here since more cluster resolutions are calculated later) 

for (x in c('A', 'B', 'C', 'D')) {
  Idents(postMALDI.list[[paste0(x,'1')]]) <- "SCT_snn_res.0.8"
  postMALDI.list[[paste0(x,'1')]]$seurat_clusters <- postMALDI.list[[paste0(x,'1')]]$SCT_snn_res.0.8
  #since some plots are grouping by seurat_clusters set this too
}

pdf("myplot.pdf")
for (x in c("A","B","C","D")){
  p1<-DimPlot(postMALDI.list[[paste0(x,'1')]], reduction = "umap", label = T)
  p2<-SpatialDimPlot(postMALDI.list[[paste0(x,'1')]], label.size = 3,image.alpha = 1)
  print(p1)
  print(p2)
#  print(p1+p2+plot_annotation(title = paste0(x,'1')))
}
dev.off()
```


```{r}
# Define the specific idents for each object
# Example: You might want to change these based on your analysis
cancer_idents <- list(
  A1 = Idents(postMALDI.list[[paste0("A",'1')]]),  # replace with actual idents for A1
#A1 =c("1", "6"),
    B1 = c("1", "7"),  # replace with actual idents for B1
  C1 = c("2", "6","7"),  # replace with actual idents for C1
  D1 = c("2", "6","5")   # replace with actual idents for D1
)

color_idents <- list(
  A1 =c("0"='#F68282','3'='#D4D915','1'='#ff9a36','2'='#2FF18B','5'='#CCB1F1',"4"="pink",'6'='#B95FBB'),  # replace with actual idents for A1
  B1 = c('1'='brown1','7'='cyan1'),  # replace with actual idents for B1
  C1 =  c('2'='brown1','6'='cyan1','7'="green"),  # replace with actual idents for C1
  D1 =c('2'='brown1','6'='cyan1','5'="green")   # replace with actual idents for D1
)

# Loop through each object and process
pdf("myplot.pdf")
for (x in c("A", "B", "C", "D")) {
  # Retrieve the specific idents for this object
  idents_to_use <- cancer_idents[[paste0(x, '1')]]
  colo=color_idents[[paste0(x, '1')]]
  # Create a subset of the Seurat object based on these idents
  cancer_subset <- subset(postMALDI.list[[paste0(x, '1')]], idents = idents_to_use)

  
  # Create a spatial dimension plot
  p <- SpatialDimPlot(cancer_subset, label = F,cols=colo)

  # Print the plot
  print(p)
}
dev.off()
# Now you have A1.cancer, B1.cancer, C1.cancer, D1.cancer and their plots stored.

```


```{r}
postMALDI.list$D1 %>% subset(ident=c("9","10")) %>% SpatialDimPlot()

postMALDI.list$B1 %>% SpatialFeaturePlot(c("Cyp3a4","Cyp2f2","Cyp2s1"))
```



```{r}

FindMarkers(postMALDI.list$B1,ident.1 = "1" ,logfc.threshold = 0.3)
```
```{r}
SpatialFeaturePlot(postMALDI.list$A1,features = c("Ltbp2","Nrep","Ppic"))
SpatialFeaturePlot(postMALDI.list$B1,features = c("Ltbp2","Nrep","Ppic"))
SpatialFeaturePlot(postMALDI.list$C1,features = c("Ltbp2","Nrep","Ppic"))
SpatialFeaturePlot(postMALDI.list$D1,features = c("Ltbp2","Nrep","Ppic"))

```

```{r}
plan("multisession", workers = 40)

```

```{r}
B1.cancer %>% SpatialFeaturePlot("Trp53")+scale_y_reverse()
B1.cancer %>% SpatialFeaturePlot("Cyp2f2")+scale_y_reverse()
B1.cancer %>% SpatialFeaturePlot("Cyp2a1")+scale_y_reverse()
postMALDI.list$B1 %>% SpatialFeaturePlot("Trp53")
postMALDI.list$B1 %>% subset(ident=c("5","7")) %>% SpatialFeaturePlot(features = "drug_conc")
postMALDI.list$B1 %>%DimPlot(reduction = "umap")
```

```{r}
library(Seurat)
library(enrichR)
library(ggplot2)

# Define your cancer_idents list
cancer_idents <- list(
    A1 = Idents(postMALDI.list[["A1"]]) %>% levels(),  # Replace with actual idents for A1
    B1 = c("1", "7"),  # Replace with actual idents for B1
    C1 = c("2", "7", "6"),  # Replace with actual idents for C1
    D1 = c("2", "6", "5")   # Replace with actual idents for D1
)
pdf("Combined_Plots.pdf")

# 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, ".")
  }
}

# Loop through each object and process
for (x in c("A", "B", "C", "D")) {
    # Retrieve the specific idents for this object
    idents_to_use <- cancer_idents[[paste0(x, '1')]]

    # Create a subset of the Seurat object based on these idents
    cancer_subset <- subset(postMALDI.list[[paste0(x, '1')]], idents = idents_to_use)

    for (ident in idents_to_use) {
        # Find markers for the current identity class compared to all others
         plan("sequential")

        find_markers_result <- FindMarkers(cancer_subset, ident.1 = ident, ident.2 = NULL)
        plan("multisession", workers = 32)

        # Proceed with enrichment analysis if markers are found
        if (nrow(find_markers_result) > 0) {
            gene_list <- find_markers_result$avg_log2FC
            names(gene_list) <- rownames(find_markers_result)
            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)
        }
    }

}
dev.off()
```


```{r}
# Additional library for pathway gene retrieval
library(clusterProfiler)

# Function to find and print overlapped genes for multiple pathways
find_and_print_overlap <- function(gene_list, pathway_names, kegg_organism) {
    # Iterate through each pathway
    for (pathway_name in pathway_names) {
        # Get the KEGG pathway gene set
        kegg_pathways <- enrichKEGG(gene = names(gene_list), organism = kegg_organism)
        pathway_genes <- kegg_pathways@result[pathway_name, "geneID"]

        # Find overlapping genes
        overlap_genes <- intersect(names(gene_list), unlist(strsplit(pathway_genes, split = "/")))

        # Print overlapping genes
        if (length(overlap_genes) > 0) {
            cat("Overlapping genes in", pathway_name, "pathway:", paste(overlap_genes, collapse = ", "), "\n")
        } else {
            cat("No overlapping genes found in the", pathway_name, "pathway.\n")
        }
    }
}

# Example usage:
# Define a list of pathways of interest
pathway_names <- c("Pathway1", "Pathway2", "Pathway3")  # Replace with actual pathway names

# Specify the cancer ident and organism
cancer_ident <- "A1"  # Replace with the actual cancer ident
kegg_organism <- "mmu"  # Use 'mmu' for mouse

# Retrieve gene list for the specific cancer ident and sample
# Replace with your function or code to get the gene list
         plan("sequential")

        find_markers_result <- FindMarkers(cancer_subset, ident.1 = ident, ident.2 = NULL)
        plan("multisession", workers = 32)

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


# Find and print overlapped genes for multiple pathways
find_and_print_overlap(gene_list, pathway_names, kegg_organism)

```


```{r}
original_gene_list <- find_markers_result$avg_log2FC

# name the vector
names(original_gene_list) <-rownames(find_markers_result)
gene_list<-na.omit(original_gene_list)

gene_list = sort(gene_list, decreasing = TRUE)
gene_list=gene_list[gene_list>0]

# Connect to Enrichr and choose a database
db <- "KEGG_2019_Mouse" # Select an appropriate database

# 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 using ggplot2
    library(ggplot2)
    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 = "Top 10 Pathways from Enrichment Analysis")
} else {
    message("No enrichment results to plot or empty results.")
}
original_gene_list <- find_markers_result$avg_log2FC

# name the vector
names(original_gene_list) <-rownames(find_markers_result)
gene_list<-na.omit(original_gene_list)

gene_list = sort(gene_list, decreasing = TRUE)
gene_list=gene_list[gene_list<0]
# DEG list

# Install and load enrichR
if (!requireNamespace("enrichR", quietly = TRUE))
    install.packages("enrichR")
library(enrichR)

# Connect to Enrichr and choose a database
db <- "KEGG_2019_Mouse" # Select an appropriate database

# 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 using ggplot2
    library(ggplot2)
    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 = "Top 10 Pathways from Enrichment Analysis")
} else {
    message("No enrichment results to plot or empty results.")
}


```

```{r}

```


```{r}
gene_list <- setNames(find_markers_result$avg_log2FC, find_markers_result %>% rownames())
gene_list=gene_list %>% sort(decreasing = TRUE)

entrez_ids <- bitr(names(gene_list), fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Mm.eg.db)
named_entrez_list <- setNames(gene_list[match(entrez_ids$SYMBOL, names(gene_list))], entrez_ids$ENTREZID)
# Perform GSEA
gsea_results <- gseKEGG(geneList = named_entrez_list, organism = "mmu", nPerm = 1000, minGSSize = 10, maxGSSize = 500)
```

```{r}
library(Seurat)
library(infercnv)

# Define reference identifiers for each sample
ref_idents <- list(
    A1 = NULL,  # Replace with actual idents for A1
    B1 = c("0", "2", "4", "6", "8"),  # Replace with actual idents for B1
    C1 = c("0", "1", "3", "4", "9"),  # Replace with actual idents for C1
    D1 = c("0", "7", "8", "4", "1", "9", "10")   # Replace with actual idents for D1
)

# Prepare a combined reference matrix
combined_ref <- NULL
combined_ref_metadata <- data.frame(cell = character(), group = character())

for (sample in names(ref_idents)) {
    if (!is.null(ref_idents[[sample]])) {
        ref_cells <- subset(postMALDI.list[[sample]], idents = ref_idents[[sample]])
        ref_cells_metadata <- data.frame(cell = colnames(ref_cells), 
                                         group = paste0(sample, "_", Idents(ref_cells)))
        combined_ref_metadata <- rbind(combined_ref_metadata, ref_cells_metadata)
        
        if (is.null(combined_ref)) {
            combined_ref <- ref_cells
        } else {
            combined_ref <- merge(combined_ref, ref_cells)
        }
    }
}

# Extract the counts matrix for the combined reference
combined_ref_counts <- GetAssayData(combined_ref, assay = "Spatial", slot = "counts")

for (sample in names(postMALDI.list)) {
    # Create output directory
    out_path <- paste0("/gpfs/data/yanailab/public/jingye_share/20231129/refagg", sample)
    dir.create(out_path, recursive = TRUE)

    # Extract counts matrix for the current sample
    sample_counts_matrix <- GetAssayData(postMALDI.list[[sample]], assay = "Spatial", slot = "counts")

    # Create a vector of group names for each cell in the sample
    sample_group_names <- paste0(sample, "_", Idents(postMALDI.list[[sample]]))

    # Determine the group names for the reference cells not in the current sample
    ref_cells_not_in_sample <- setdiff(combined_ref_metadata$cell, colnames(sample_counts_matrix))
    ref_group_names <- combined_ref_metadata$group[combined_ref_metadata$cell %in% ref_cells_not_in_sample]

    # Combine the current sample's counts with the reference counts
    all_counts_matrix <- cbind(sample_counts_matrix, combined_ref_counts[, colnames(combined_ref_counts) %in% ref_cells_not_in_sample])

    # Create the annotations file
    annotations_file <- data.frame(
        group = c(sample_group_names, ref_group_names)
    )
rownames(annotations_file)=colnames(all_counts_matrix)
    # Create inferCNV object with the annotations
    infercnv_obj <- CreateInfercnvObject(
        raw_counts_matrix = all_counts_matrix,
        annotations_file = annotations_file %>% as.matrix(),
        gene_order_file = 'mouse_gencode.GRCm38.p6.vM25.basic.annotation.by_gene_name.infercnv_positions',
        ref_group_names = ref_group_names %>% unique()
    )
    # Run inferCNV analysis
    options(scipen = 100)
 infercnv_obj <- infercnv::run(
    infercnv_obj,
    resume_mode = FALSE,
    cutoff = 0.1,
    out_dir = out_path,
    cluster_by_groups = FALSE,
    denoise = TRUE,
    HMM = FALSE,
 #   HMM_type = "i6",
   # analysis_mode = "subclusters",
#    HMM_report_by = "subcluster",
    scale_data = TRUE,
    tumor_subcluster_partition_method = "qnorm",
   # k_nn = 10,
    plot_steps = FALSE,
    debug = FALSE
  )
}




```



```{r}
library(Seurat)
library(infercnv)

# Define reference identifiers for each sample
ref_idents <- list(
    A1 = NULL,  # Replace with actual idents for A1
    B1 = c("0", "2", "4", "6", "8"),  # Replace with actual idents for B1
    C1 = c("0", "1", "3", "4", "9"),  # Replace with actual idents for C1
    D1 = c("0", "7", "8", "4", "1", "9", "10")   # Replace with actual idents for D1
)

# Prepare a combined reference matrix
combined_ref <- NULL
combined_ref_metadata <- data.frame(cell = character(), group = character())

for (sample in names(ref_idents)) {
    if (!is.null(ref_idents[[sample]])) {
        ref_cells <- subset(postMALDI.list[[sample]], idents = ref_idents[[sample]])
        ref_cells$sample_origin <- paste0(sample, "_", Idents(ref_cells))
        ref_group_names_list <- c(ref_group_names_list, unique(ref_cells$sample_origin))
        ref_cells_metadata <- data.frame( group = paste0(sample, "_", Idents(ref_cells)))
        rownames( ref_cells_metadata)=colnames(ref_cells)
        combined_ref_metadata <- rbind(combined_ref_metadata, ref_cells_metadata)
        
        if (is.null(combined_ref)) {
            combined_ref <- ref_cells
        } else {
            combined_ref <- merge(combined_ref, ref_cells)
        }
    }
}

# Extract the counts matrix for the combined reference
combined_ref_counts <- GetAssayData(combined_ref, assay = "Spatial", slot = "counts")

# Loop through each sample for inferCNV analysis
for (sample in names(postMALDI.list)) {
    # Create output directory
    out_path <- paste0("/path/to/output/", sample)
    dir.create(out_path, recursive = TRUE)

    # Extract counts matrix and cell identities for the current sample
    sample_counts_matrix <- GetAssayData(postMALDI.list[[sample]], assay = "Spatial", slot = "counts")
    sample_metadata <- data.frame( group = Idents(postMALDI.list[[sample]]))
rownames(sample_metadata)=colnames(sample_counts_matrix)
    # Combine the current sample's counts with the reference counts
    all_counts_matrix <- cbind(sample_counts_matrix, combined_ref_counts[, !(colnames(combined_ref_counts) %in% colnames(sample_counts_matrix))])

    # Combine metadata for the inferCNV analysis
    all_metadata <- rbind(sample_metadata, combined_ref_metadata)

    # Create inferCNV object
    infercnv_obj <- CreateInfercnvObject(
        raw_counts_matrix = all_counts_matrix,
        annotations_file = all_metadata,
        gene_order_file = 'mouse_gencode.GRCm38.p6.vM25.basic.annotation.by_gene_name.infercnv_positions'
        ref_group_names=
    )

    # Run inferCNV analysis
    infercnv_obj <- infercnv::run(
        infercnv_obj,
        out_dir = out_path,
        # Other inferCNV parameters as needed
    )
}

```

```{r}
library(Seurat)
library(infercnv)
library(dplyr)

ref_idents <- list(
    A1 = NULL,  # Replace with actual idents for A1
    B1 = c("0", "2", "4", "6", "8"),  # Replace with actual idents for B1
    C1 = c("0", "1", "3", "4", "9"),  # Replace with actual idents for C1
    D1 = c("0", "7", "8", "4", "1", "9", "10")   # Replace with actual idents for D1
)

# Prepare a combined reference matrix
combined_ref <- NULL
combined_ref_metadata <- data.frame()

# Process each sample
for (sample in names(ref_idents)) {
    if (!is.null(ref_idents[[sample]])) {
        ref_cells <- subset(postMALDI.list[[sample]], idents = ref_idents[[sample]])
        ref_cells_metadata <- data.frame(cell = colnames(ref_cells), 
#                                         group = paste0(sample, "_", Idents(ref_cells)))
      group = rep("ref",length(colnames(ref_cells))))
        combined_ref_metadata <- rbind(combined_ref_metadata, ref_cells_metadata)
        
        if (is.null(combined_ref)) {
            combined_ref <- ref_cells
        } else {
            combined_ref <- merge(combined_ref, ref_cells)
        }
    }
}

# Extract the counts matrix for the combined reference
combined_ref_counts <- GetAssayData(combined_ref, assay = "Spatial", slot = "counts")

# Subsetting process for reference cells
subset_ref_cells <- list()
for (cluster in unique(combined_ref_metadata$group)) {
    cells_in_cluster <- which(combined_ref_metadata$group == cluster)
    if (length(cells_in_cluster) > 30) {
        subset_ref_cells[[cluster]] <- sample(cells_in_cluster, 30)
    } else {
        subset_ref_cells[[cluster]] <- cells_in_cluster
    }
}
subset_ref_cells <- unlist(subset_ref_cells)
subset_combined_ref <- combined_ref[, subset_ref_cells, drop = FALSE]

# Updated reference counts matrix
combined_ref_counts <- GetAssayData(subset_combined_ref, assay = "Spatial", slot = "counts")

# Loop through each sample for inferCNV analysis
for (sample in names(postMALDI.list)) {
    # Create output directory
    out_path <- paste0("/gpfs/data/yanailab/public/jingye_share/20231129/refagg", sample)
    dir.create(out_path, recursive = TRUE)

    # Subset process for each sample
    subset_sample_cells <- list()
    for (cluster in unique(Idents(postMALDI.list[[sample]]))) {
        cells_in_cluster <- WhichCells(postMALDI.list[[sample]], idents = cluster)
        if (length(cells_in_cluster) > 30) {
            subset_sample_cells[[cluster]] <- sample(cells_in_cluster, 30)
        } else {
            subset_sample_cells[[cluster]] <- cells_in_cluster
        }
    }
    subset_sample_cells <- unlist(subset_sample_cells)
    subset_seurat <- postMALDI.list[[sample]][, subset_sample_cells]


    # Extract counts matrix for the subset sample
    subset_counts_matrix <- GetAssayData(subset_seurat, assay = "Spatial", slot = "counts")

    # Combine the subset sample's counts with the reference counts
    all_counts_matrix <- cbind(subset_counts_matrix, combined_ref_counts[, colnames(combined_ref_counts) %in% setdiff(combined_ref_metadata$cell, colnames(subset_counts_matrix))])

    # Create the annotations file
    annotations_file <- data.frame( group = c(Idents(subset_seurat), combined_ref_metadata$group[combined_ref_metadata$cell %in% colnames(all_counts_matrix)])
    )
rownames(annotations_file)=colnames(all_counts_matrix)
    # Create inferCNV object with the annotations
    infercnv_obj <- CreateInfercnvObject(
        raw_counts_matrix = all_counts_matrix ,
        annotations_file = annotations_file %>% as.matrix(),
        gene_order_file = 'mouse_gencode.GRCm38.p6.vM25.basic.annotation.by_gene_name.infercnv_positions',
#        ref_group_names = combined_ref_metadata$group %>% unique()
ref_group_names ="ref"
    )

    # Run inferCNV analysis
    options(scipen = 100)
 infercnv_obj <- infercnv::run(
    infercnv_obj,
    resume_mode = FALSE,
    cutoff = 0.1,
    out_dir = out_path,
    cluster_by_groups = FALSE,
    denoise = TRUE,
    HMM = FALSE,
 #   HMM_type = "i6",
    analysis_mode = "subclusters",
  #  HMM_report_by = "subcluster",
    scale_data = TRUE,
    tumor_subcluster_partition_method = "qnorm",
    k_nn = 10,
    plot_steps = FALSE,
    debug = FALSE
  )

}

```
```{r}
library(Seurat)
library(infercnv)
library(dplyr)

# Define reference identifiers for each sample
ref_idents <- list(
  A1 = NULL,  # Replace with actual idents for A1
  B1 = c("0", "2", "4", "6", "8"),  # Replace with actual idents for B1
  C1 = c("0", "1", "3", "4", "9"),  # Replace with actual idents for C1
  D1 = c("0", "7", "8", "4", "1", "9", "10")   # Replace with actual idents for D1
)


# Prepare a combined reference matrix
combined_ref <- NULL
### from linda
for(sam in names(ref_idents)){
  if(!is.null(ref_idents[[sam]])){
    ref_cells <- subset(postMALDI.list[[sam]], idents = ref_idents[[sam]])
    ref_cells$inferCNV_group<-paste0(sam, "_", Idents(ref_cells))
    if (is.null(combined_ref)) {
      combined_ref <- ref_cells
    } else {
      combined_ref <- merge(combined_ref, ref_cells)
    }
  }
}

for (sam in names(postMALDI.list)) {
#for (sam in c("B1")) {
  # Create output directory
  out_path <- paste0("/gpfs/data/yanailab/public/jingye_share/20231129/refaggg", sam)
  dir.create(out_path, recursive = TRUE)
  if(!is.null(ref_idents[[sam]])){
    obj<-subset(postMALDI.list[[sam]], invert=T, cell=WhichCells(postMALDI.list[[sam]], 
 cells=colnames(combined_ref)[combined_ref$orig.ident==sam]))
  } else{
    obj<-postMALDI.list[[sam]]
  }
  obj$inferCNV_group<-paste0(sam, "_", Idents(obj))
  obj<-merge(obj, combined_ref)
  st_meta<-as.data.frame(matrix(nrow=ncol(obj),ncol=1))
  rownames(st_meta)<-colnames(obj)
  st_meta$V1<-obj$inferCNV_group
  infercnv_obj <- CreateInfercnvObject(
    raw_counts_matrix = GetAssayData(obj, assay = "Spatial", slot = "counts"),
    annotations_file = st_meta,
    gene_order_file = '/gpfs/data/yanailab/public/jingye_share/20231129/mouse_gencode.GRCm38.p6.vM25.basic.annotation.by_gene_name.infercnv_positions', #human genecode
    ref_group_names = unique(combined_ref$inferCNV_group)
  )
  options(scipen = 100)
  infercnv_obj <- infercnv::run(
    infercnv_obj,
    resume_mode = FALSE,
    cutoff = 0.1,
    out_dir = out_path,
    cluster_by_groups = FALSE,
    denoise = F,
    HMM = FALSE,
    #   HMM_type = "i6",
    # analysis_mode = "subclusters",
    #    HMM_report_by = "subcluster",
    scale_data = TRUE,
    tumor_subcluster_partition_method = "qnorm",
    # k_nn = 10,
  #  k_obs_groups = 4,
    plot_steps = FALSE,
    debug = FALSE
  )
  
}
```



