稳健和准确top1的空转反卷积算法RCTD实战：GSE242249

###
### Create: Jianming Zeng
### Date:  2023-12-31  
### Email: jmzeng1314@163.com
### Blog: http://www.bio-info-trainee.com/
### Forum:  http://www.biotrainee.com/thread-1376-1-1.html
### CAFS/SUSTC/Eli Lilly/University of Macau
### Update Log: 2023-12-31   First version 
### Update Log: 2024-12-09   by juan zhang (492482942@qq.com)
### 
rm(list=ls())
library(ggsci)
library(dplyr) 
library(future)
library(Seurat)
library(clustree)
library(cowplot)
library(data.table)
library(ggplot2)
library(patchwork)
library(stringr)
library(qs)
library(Matrix)

# 创建目录
getwd()
gse <- "GSE172357"
dir.create(gse)

# 方式一：标准文件夹
###### step1: 导入数据 ######   
samples <- list.files("GSE172357/", pattern = "h5$",recursive = F, full.names = F)
samples

# 挑选文献中的12个
files <- c(samples[1:9],samples[15],samples[17],samples[19])
files

scRNAlist <- lapply(files, function(pro){
#pro <- files[1]
print(pro)
  counts <- Read10X_h5(paste0("GSE172357/",pro))
  name <- str_split(pro, pattern = "_Via_|_Fcol_|_filtered",simplify = T,n = 2)[,1]
print(name)
  sce <- CreateSeuratObject(counts, project=name, min.cells = 3)
return(sce)
})
names(scRNAlist) <-  str_split(files, pattern = "_Via_|_Fcol_|_filtered",simplify = T,n = 2)[,1]
scRNAlist

# merge
sce.all <- merge(scRNAlist[[1]], y=scRNAlist[-1], add.cell.ids=names(scRNAlist))
sce.all <- JoinLayers(sce.all) # seurat v5
sce.all

# 查看特征
as.data.frame(sce.all@assays$RNA$counts[1:10, 1:2])
head(sce.all@meta.data, 10)
table(sce.all$orig.ident) 

library(qs)
qsave(sce.all, file="GSE172357/sce.all.qs")

rm(list=ls())
library(spacexr)
library(Matrix)
library(doParallel)
library(ggplot2)

# 输出文件夹，可以换成自己指定的
savedir <- 'RCTD_results_Run'
if(!dir.exists(savedir))
  dir.create(savedir)
knitr::opts_chunk$set(collapse = TRUE,comment = "#>",cache = TRUE,out.width = "100%" )

load("../2024-GSE242249-空间单细胞-良性前列腺增生/spatial.Rdata")
counts <- as.matrix(spatial@assays$Spatial$counts)
counts[1:10,1:5]

# 读取坐标
centroids <- spatial[["slice1"]]@boundaries$centroids
coords <- setNames(as.data.frame(centroids@coords), c("x", "y"))
rownames(coords) <- centroids@cells
head(coords)

# 计算每个位置的umi
nUMI <- colSums(counts) # In this case, total counts per pixel is nUMI
head(nUMI)

# 构建一个对象
puck <- SpatialRNA(coords, counts, nUMI)
puck

barcodes <- colnames(puck@counts) # pixels to be used (a list of barcode names). 
plot_puck_continuous(puck, barcodes, puck@nUMI, 
                     ylimit = c(0,round(quantile(puck@nUMI,0.9))), 
                     size = 2,
                     title ='plot of nUMI') 

# 看一下构建好的SpatialRNA对象
str(puck)

# count值
sce <- readRDS("../../03-scRNA/2021-GSE172357-12-BPH/3-check-by-0.5/sce.all.int.rds")
counts <- as.matrix(sce@assays$RNA$counts)
counts[1:5,1:5]

# 细胞类型：
cell_types <- sce$celltype
levels(cell_types)[-9]
cell_types <- factor(cell_types, levels = levels(cell_types)[-9])
head(cell_types)
table(cell_types)

# umi
nUMI <- sce$nCount_RNA
head(nUMI)

# 构建对象
reference <- Reference(counts, cell_types, nUMI)
str(reference)

# 运行
myRCTD <- create.RCTD(puck, reference, max_cores = 60)
myRCTD <- run.RCTD(myRCTD, doublet_mode = 'full')
saveRDS(myRCTD,'myRCTD_visium_full.rds')

Begin: process_cell_type_info
process_cell_type_info: number of cells in reference: 64142
process_cell_type_info: number of genes in reference: 27508