官网链接：https://github.com/sunduanchen/Scissor
官方教程：https://sunduanchen.github.io/Scissor/vignettes/Scissor_Tutorial.html
Conda环境设置：
mamba install -c conda-forge -c bioconda r-base==4.1.2 r-qs r-devtools r-rlang==1.1.0 bioconductor-preprocesscore r-progress -y运行代码：
library(Scissor)
library(Seurat)
library(qs)
library(readxl)
library(dplyr)
library(tidyr)
library(Matrix)
## 57.0 设置参数 ####
proteomics_path <- "/Users/zhoukaiwen/Desktop/Breast_Phyllodes_Tumor/Bioinfo/Proteomics/pg_mtx_minhalf_drop.txt"
group_path <- "/Users/zhoukaiwen/Desktop/Breast_Phyllodes_Tumor/Bioinfo/sample_group.xlsx"
sc_path <- "/Users/zhoukaiwen/Desktop/Breast_Phyllodes_Tumor/Bioinfo/R_base/26samples_Merged_AllCells_Annotated_Final_CellType_Minor.qs"
merge_downgraded <- F # TRUE：Stable + Downgraded；FALSE：仅 Stable，排除 Downgraded
proteomics_tissue <- "M" # 选择 "E" 或 "M"；两者都用时改为 c("E", "M")
stopifnot(length(proteomics_tissue) > 0, all(proteomics_tissue %in% c("E", "M")))
scissor_celltypes <- "Fibroblasts"
# 服务器分析原先指定的全部大类时，将上一行改为：
# scissor_celltypes <- c("Fibroblasts", "Pericytes", "VSMCs", "Endothelials", "BCells", "TCells", "Myeloids")
scissor_alpha <- 0.05 # 预先固定，不根据 MP2 富集的 P 值挑选参数
## 57.1 提取蛋白矩阵 ####
# 蛋白矩阵已经 log2 转换并填补缺失，不再重复 log2；按所选 E/M 筛选原发灶。
proteomics <- as.matrix(read.delim(proteomics_path, row.names = 1, check.names = FALSE))
sample_group <- readxl::read_xlsx(group_path)
nonprogression_groups <- if (merge_downgraded) c("Stable", "Downgraded") else "Stable"
primary_group <- sample_group %>%
dplyr::filter(SampleTime == "Primary",
TissueType %in% proteomics_tissue,
Group_Progression %in% c("Progression", nonprogression_groups)) %>%
dplyr::mutate(Scissor_group = ifelse(Group_Progression == "Progression", "Progression", "Non-progression"))
stopifnot(!anyDuplicated(primary_group$SampleName),
all(primary_group$SampleName %in% colnames(proteomics)))
bulk_dataset <- proteomics[, primary_group$SampleName, drop = FALSE]
phenotype <- as.integer(primary_group$Scissor_group == "Progression") # 0=不进展，1=进展
stopifnot(identical(colnames(bulk_dataset), primary_group$SampleName),
setequal(phenotype, 0:1), all(is.finite(bulk_dataset)))
print(table(primary_group$Group_Progression, primary_group$TissueType))
## 57.2 提取单细胞矩阵 ####
scissor_seu <- qs::qread(sc_path)
scissor_seu <- subset(scissor_seu, subset = CellType_Major %in% scissor_celltypes)
print(table(scissor_seu$CellType_Major))
sc_dataset <- GetAssayData(scissor_seu, assay = "RNA", layer = "counts")
sc_dataset <- Seurat_preprocessing(sc_dataset, verbose = T)
class(sc_dataset)
names(sc_dataset)
DimPlot(sc_dataset,
reduction = 'umap',
label = T,
label.size = 10)
Scissor_V5 <- function (bulk_dataset, sc_dataset, phenotype, tag = NULL, alpha = NULL,
cutoff = 0.2, family = c("gaussian", "binomial", "cox"),
Save_file = "Scissor_inputs.RData", Load_file = NULL)
{
library(Seurat)
library(Matrix)
library(preprocessCore)
if (is.null(Load_file)) {
common <- intersect(rownames(bulk_dataset), rownames(sc_dataset))
if (length(common) == 0) {
stop("There is no common genes between the given single-cell and bulk samples.")
}
if (class(sc_dataset) == "Seurat") {
sc_exprs <- as.matrix(sc_dataset@assays$RNA@layers$data)
network <- as.matrix(sc_dataset@graphs$RNA_snn)
}
else {
sc_exprs <- as.matrix(sc_dataset)
Seurat_tmp <- CreateSeuratObject(sc_dataset)
Seurat_tmp <- FindVariableFeatures(Seurat_tmp, selection.method = "vst",
verbose = F)
Seurat_tmp <- ScaleData(Seurat_tmp, verbose = F)
Seurat_tmp <- RunPCA(Seurat_tmp, features = VariableFeatures(Seurat_tmp),
verbose = F)
Seurat_tmp <- FindNeighbors(Seurat_tmp, dims = 1:10,
verbose = F)
network <- as.matrix(Seurat_tmp@graphs$RNA_snn)
}
diag(network) <- 0
network[which(network != 0)] <- 1
dataset0 <- cbind(bulk_dataset[common, ], sc_exprs[common,
])
dataset1 <- normalize.quantiles(dataset0)
rownames(dataset1) <- rownames(dataset0)
colnames(dataset1) <- colnames(dataset0)
Expression_bulk <- dataset1[, 1:ncol(bulk_dataset)]
Expression_cell <- dataset1[, (ncol(bulk_dataset) +
1):ncol(dataset1)]
X <- cor(Expression_bulk, Expression_cell)
quality_check <- quantile(X)
print("|**************************************************|")
print("Performing quality-check for the correlations")
print("The five-number summary of correlations:")
print(quality_check)
print("|**************************************************|")
if (quality_check[3] < 0.01) {
warning("The median correlation between the single-cell and bulk samples is relatively low.")
}
if (family == "binomial") {
Y <- as.numeric(phenotype)
z <- table(Y)
if (length(z) != length(tag)) {
stop("The length differs between tags and phenotypes. Please check Scissor inputs and selected regression type.")
}
else {
print(sprintf("Current phenotype contains %d %s and %d %s samples.",
z[1], tag[1], z[2], tag[2]))
print("Perform logistic regression on the given phenotypes:")
}
}
if (family == "gaussian") {
Y <- as.numeric(phenotype)
z <- table(Y)
if (length(z) != length(tag)) {
stop("The length differs between tags and phenotypes. Please check Scissor inputs and selected regression type.")
}
else {
tmp <- paste(z, tag)
print(paste0("Current phenotype contains ",
paste(tmp[1:(length(z) - 1)], collapse = ", "),
", and ", tmp[length(z)], " samples."))
print("Perform linear regression on the given phenotypes:")
}
}
if (family == "cox") {
Y <- as.matrix(phenotype)
if (ncol(Y) != 2) {
stop("The size of survival data is wrong. Please check Scissor inputs and selected regression type.")
}
else {
print("Perform cox regression on the given clinical outcomes:")
}
}
save(X, Y, network, Expression_bulk, Expression_cell,
file = Save_file)
}
else {
load(Load_file)
}
if (is.null(alpha)) {
alpha <- c(0.005, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5,
0.6, 0.7, 0.8, 0.9)
}
for (i in 1:length(alpha)) {
set.seed(123)
fit0 <- APML1(X, Y, family = family, penalty = "Net",
alpha = alpha[i], Omega = network, nlambda = 100,
nfolds = min(10, nrow(X)))
fit1 <- APML1(X, Y, family = family, penalty = "Net",
alpha = alpha[i], Omega = network, lambda = fit0$lambda.min)
if (family == "binomial") {
Coefs <- as.numeric(fit1$Beta[2:(ncol(X) + 1)])
}
else {
Coefs <- as.numeric(fit1$Beta)
}
Cell1 <- colnames(X)[which(Coefs > 0)]
Cell2 <- colnames(X)[which(Coefs < 0)]
percentage <- (length(Cell1) + length(Cell2))/ncol(X)
print(sprintf("alpha = %s", alpha[i]))
print(sprintf("Scissor identified %d Scissor+ cells and %d Scissor- cells.",
length(Cell1), length(Cell2)))
print(sprintf("The percentage of selected cell is: %s%%",
formatC(percentage * 100, format = "f", digits = 3)))
if (percentage < cutoff) {
break
}
cat("\n")
}
print("|**************************************************|")
return(list(para = list(alpha = alpha[i], lambda = fit0$lambda.min,
family = family), Coefs = Coefs, Scissor_pos = Cell1,
Scissor_neg = Cell2))
}
tag <- c('Stable', 'Progression-related')
infos1 <- Scissor_V5(bulk_dataset,
sc_dataset,
phenotype,
tag = tag,
alpha = scissor_alpha,
family = "binomial",
Save_file = 'Scissor_Fibro_PvsS.RData')