library(Seurat)
obj.raw = Read10X('raw_data/Path/filtered_feature_bc_matrix/')
obj.srt = CreateSeuratObject(counts = obj.raw, project = 'project_name')

import numpy as np
import pandas as pd
import scanpy as sc
import scrublet as scr

def run_scrublet(input_file, output_file, expected_doublet_rate=0.1):
    df = pd.read_csv(input_file, header=0, index_col=0)
    adata = sc.AnnData(df)
    
    # Set the expected_doublet_rate parameter
    sc.external.pp.scrublet(adata, expected_doublet_rate=expected_doublet_rate)
    
    # Save the observation (results) to CSV
    adata.obs.to_csv(output_file)

input_file = 'input.count.csv'
output_file = 'output.scr.csv'

# Adjust the expected_doublet_rate parameter
expected_doublet_rate = 0.1
run_scrublet(input_file, output_file, expected_doublet_rate)

## import scrublet result
df= read.csv('path/output.scr.csv', row.names = 1)

# If the previous process changes the name of cell ID, 
# Use the following code to correct it. 
rownames(df) = gsub(pattern = '_', replacement = '-' ,rownames(df))

## add doublet info to the srt obj

obj.srt[['doublet_score']] = df[rownames(obj.srt@meta.data),]$doublet_score
obj.srt[['predicted_doublet']] =df[rownames(obj.srt@meta.data),]$predicted_doublet

obj.srt %>% saveRDS("path/saved.obj.rds")

# You might want to check the number of doublets in the data
obj.srt@meta.data %>% select(predicted_doublet) %>% table()

## perform default analysis
perform_default_analysis <- function(obj.srt, n_features = 2000, n_pcs = 20, 
                                     dims_for_neighbors = 1:20, 
                                     resolutions = c(0.1,0.5,1), 
                                     umap_dims = 1:20) {
  # Step 1: Find variable features
  obj.srt <- FindVariableFeatures(obj.srt, 
                                  selection.method = 'vst', 
                                  nfeatures = n_features)
  
  # Step 2: Scale and normalize data
  all_genes <- rownames(obj.srt)
  obj.srt <- NormalizeData(obj.srt)
  obj.srt <- ScaleData(obj.srt, features = all_genes)
  
  # Step 3: Run PCA
  obj.srt <- RunPCA(obj.srt, 
                    features = VariableFeatures(object = obj.srt), npcs = n_pcs)
  
  # Step 4: Find neighbors
  obj.srt <- FindNeighbors(obj.srt, dims = dims_for_neighbors)
  
  # Step 5: Find clusters
  obj.srt <- FindClusters(obj.srt, resolution = resolutions)
  
  # Step 6: Run UMAP
  obj.srt <- RunUMAP(obj.srt, dims = umap_dims)
  
  # Return the Seurat object with analysis results
  return(obj.srt)
}

# apply
obj.srt <- perform_default_analysis(obj.srt)

# Define what groups to be integrated 
groups_to_be_integrated = "Timepoint"

# Split object and make them a list 
obj.list <- SplitObject(obj.srt, split.by = groups_to_be_integrated)

# Normalize and identify variable features for each dataset independently
# selection method and number of features can be modified 

intgr.list <- lapply(X = obj.list, FUN = function(x) {
  x <- NormalizeData(x)
  x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

# select features that are repeatedly variable across datasets for integration
features <- SelectIntegrationFeatures(object.list = intgr.list)

# Perform integration
obj.anchors <- FindIntegrationAnchors(object.list = intgr.list, 
                                      anchor.features = features)

# this command creates an 'integrated' data assay
obj.intgr <- IntegrateData(anchorset = obj.anchors)

# Perform an integrated analysis

# specify that we will perform downstream analysis on the corrected data note that the
# original unmodified data still resides in the 'RNA' assay
# To use integrated data, change the default assay to "integrated"
DefaultAssay(obj.intgr) <- "integrated"

# Run the standard workflow for visualization and clustering
obj.intgr <- ScaleData(obj.intgr, verbose = FALSE)
obj.intgr <- RunPCA(obj.intgr, npcs = 30, verbose = FALSE)
obj.intgr <- RunUMAP(obj.intgr, reduction = "pca", dims = 1:30)
obj.intgr <- FindNeighbors(obj.intgr, reduction = "pca", dims = 1:30)
obj.intgr <- FindClusters(obj.intgr, resolution = c(0.1,0.2,0.5,0.8))

obj.intgr %>% saveRDS("path/saved.obj.rds")

# Visualization
DimPlot(obj.intgr, reduction = "umap", group.by = groups_to_be_integrated) 

s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes
obj.srt <- CellCycleScoring(obj.srt, 
                            s.features = s.genes, g2m.features = g2m.genes)

# load gene set preparation function
source("https://raw.githubusercontent.com/IanevskiAleksandr/sc-type/master/R/gene_sets_prepare.R")
# load cell type annotation function
source("https://raw.githubusercontent.com/IanevskiAleksandr/sc-type/master/R/sctype_score_.R")

# DB file
db_ <- "https://raw.githubusercontent.com/IanevskiAleksandr/sc-type/master/ScTypeDB_full.xlsx";
tissue <- "Immune system" # e.g. Immune system,Pancreas,Liver,Kidney,Brain,Lung

# prepare gene sets
gs_list <- gene_sets_prepare(db_, tissue)

# extract scaled scRNA-seq matrix
scRNAseqData_scaled <- as.matrix(obj.srt[["RNA"]]@scale.data)

# run ScType
es.max <- sctype_score(scRNAseqData = scRNAseqData_scaled, 
                       scaled = TRUE, 
                       gs = gs_list$gs_positive, 
                       gs2 = gs_list$gs_negative)

# Merge by cluster
cL_results <- do.call("rbind", 
  lapply(unique(obj.srt@meta.data$seurat_clusters), function(cl) {
    cl_data <- obj.srt@meta.data[obj.srt@meta.data$seurat_clusters == cl, ]
    es_sum <- rowSums(es.max[, rownames(cl_data)])
    es_sorted <- sort(es_sum, decreasing = TRUE)
    top_es <- head(data.frame(cluster = cl, 
                              type = names(es_sorted), 
                              scores = es_sorted, 
                              ncells = nrow(cl_data)), 10)
    return(top_es)
  })
)

# Calculate top scores by cluster
sctype_scores <- cL_results %>% 
  group_by(cluster) %>% 
  top_n(n = 1, wt = scores)

# Set low-confidence clusters to "Unknown"
threshold <- sctype_scores$ncells / 4
sctype_scores$type[sctype_scores$scores < threshold] <- "Unknown"

# overlay the identified cell types on UMAP plot
obj.srt@meta.data$sctype_classification = ""

# Assign cell type classification to meta data
unique_clusters <- unique(sctype_scores$cluster)

for (j in unique_clusters) {
  # Extract the type for the current cluster
  cl_type <- sctype_scores[sctype_scores$cluster == j, ]
  type_label <- as.character(cl_type$type[1])
  
  # Assign type label to all cells in the cluster
  obj.srt@meta.data$sctype_classification[obj.srt@meta.data$seurat_clusters == j] <- type_label
}

# Draw plot 
DimPlot(obj.srt, reduction = "umap", label = TRUE, repel = TRUE, 
        group.by = 'sctype_classification')        

library(Seurat)
library(AUCell)
library(GSEABase)

# import data from single-cell experiment
exprMatrix <- obj.srt@assays$RNA@data
# Convert to sparse:
exprMatrix <- as(exprMatrix, "dgCMatrix")

# Create genesets to test
# getting them from msigDB or make custom ones 
gmt 

library(clusterProfiler)
hallmark <- msigdbr::msigdbr(species = "Mus musculus", category = "H") %>% 
  dplyr::select(gs_name, gene_symbol)

perform_GSEA <- function(res, ref, pvalueCutoff = 1) {
  ranking <- function(res) {
    df <- res$avg_log2FC
    names(df) <- rownames(res)
    df <- sort(df, decreasing = TRUE)
    return(df)
  }
  
  ranked.res <- ranking(res)
  
  x <- clusterProfiler::GSEA(geneList = ranked.res,
                             TERM2GENE = ref,
                             pvalueCutoff = pvalueCutoff,
                             pAdjustMethod = "BH",
                             verbose = TRUE,
                             seed = TRUE)
  
  result <- x@result %>% arrange(desc(NES))
  result <- result[, c('NES', 'pvalue', 'p.adjust', 'core_enrichment', 'ID')]
  return(result)
}

# Application 
gsea.res = perform_GSEA(res = deg.cluster3, ref = hallmark)

# Import the necessary libraries.
import pandas as pd
import magic
import csv

# Load the CSV file into a DataFrame.
data_path = '~/path/data.csv'
x = pd.read_csv(data_path)

# Remove the 'Unnamed: 0'
del x['Unnamed: 0']

# Initialize the MAGIC algorithm. the number of nearest neighbors (knn) is set to 5, and the diffusion time (t) is set to 3.
magic_operator = magic.MAGIC(knn=5, t=3)

# Apply the MAGIC algorithm to the data.
x_magic = magic_operator.fit_transform(x)

# Save the data to a CSV file.
output_path = '~/path/data.magic.csv'
x_magic.to_csv(output_path)