library(Seurat)
library(SeuratDisk)
library(SeuratData)
library(scMINER)
###########Create Seurat Object using Spotiphy-derived iscRNA data##############
# Read h5ad file
ST_decomposition <- readH5AD("ST_decomposition.h5ad")
ST_decomposition.seurat <- as.Seurat(ST_decomposition, counts = "X", data = "X")
ST_decomposition.seurat$samplename <- 'samplename'
# saveRDS(ST_decomposition.seurat, file = "ST_decomposition.seurat.rds")
#Common Seurat pipeline#
ST_decomposition.seurat <- NormalizeData(ST_decomposition.seurat, normalization.method = "LogNormalize", scale.factor=1e6)
ST_decomposition.seurat <- FindVariableFeatures(ST_decomposition.seurat, selection.method = "vst", nfeatures = 2000)
ST_decomposition.seurat <- ScaleData(ST_decomposition.seurat, features = rownames(ST_decomposition.seurat))
ST_decomposition.seurat <- RunPCA(ST_decomposition.seurat, features = VariableFeatures(object = ST_decomposition.seurat))
ST_decomposition.seurat <- FindNeighbors(ST_decomposition.seurat, dims = 1:30)
ST_decomposition.seurat <- FindClusters(ST_decomposition.seurat, resolution = 0.5)
ST_decomposition.seurat <- RunUMAP(ST_decomposition.seurat, dims = 1:30)
DimPlot(ST_decomposition.seurat, reduction = "umap")
# saveRDS(ST_decomposition.seurat, file = "ST_decomposition.seurat.rds")
#Find DE genes of Clusters
markers <- FindAllMarkers(ST_decomposition.seurat, only.pos = TRUE)
# Output cluster information for Spatial Layout in Python
meta.data <- ST_decomposition.seurat@meta.data
meta.data$ID <- rownames(meta.data)
out2excel(meta.data,out.xlsx = "ST_decomposition.seurat_meta.xlsx")
####################Create SparseEset Object using Seurat Object################
#build meta.data
meta.data <- ST_decomposition.seurat@meta.data
#create SparseEset
ST_decomposition.eset <- createSparseEset(input_matrix=ST_decomposition.seurat@assays$originalexp@counts,cellData=meta.data)
colnames(ST_decomposition.eset)[1:(ncol(ST_decomposition.eset))] <- sapply(1:(ncol(ST_decomposition.eset)), function(l) paste0("Cell_", colnames(ST_decomposition.eset)[l]))
# saveRDS(ST_decomposition.eset, file = "ST_decomposition.eset.rds")
#Common scMINER pipeline# For more details: https://jyyulab.github.io/scMINER/index.html
#QC
drawSparseEsetQC(input_eset = ST_decomposition.eset, output_html_file = "ST_decomposition.eset_rawCount.html", overwrite = FALSE)
ST_decomposition.eset <- filterSparseEset(ST_decomposition.eset, filter_mode = "auto", filter_type = "both")
#Normalization
ST_decomposition.eset.log2 <- normalizeSparseEset(ST_decomposition.eset, scale_factor = 1000000, log_base = 2, log_pseudoCount = 1)
# saveRDS(ST_decomposition.eset.log2, file = "ST_decomposition.eset.log2.rds")
#Generate MICA input
generateMICAinput(input_eset= ST_decomposition.eset.log2, output_file = "ST_decomposition.eset.log2_MICA_input.txt", overwrite = TRUE)
#Read in MICA output for annotation
ST_decomposition.eset.log2 <- addMICAoutput(input_eset = ST_decomposition.eset.log2, visual_method = "umap",
mica_output_file = "ST_decomposition.eset.log2/clustering_UMAP_euclidean_20_2.0.txt")
MICAplot(input_eset = ST_decomposition.eset.log2, X = "UMAP_1", Y = "UMAP_2", color_by = "clusterID", point.size = 0.5, fontsize.cluster_label = 4)
# saveRDS(ST_decomposition.eset.log2, file = "ST_decomposition.eset.log2.rds")
#Get DE genes of Clusters
de_res <- getDE(input_eset = ST_decomposition.eset.log2, group_by = "clusterID", use_method = "limma")
# Output cluster information for Spatial Layout in Python
meta.data <- Biobase::pData(ST_decomposition.eset.log2)
meta.data$ID <- rownames(meta.data)
out2excel(meta.data,out.xlsx = "ST_decomposition.eset.log2_meta.xlsx")