Scissor：聯(lián)合表型數(shù)據(jù)，Bulk-seq和scRNA的分析軟件

前言

作者開(kāi)發(fā)Scissor的目的是結(jié)合bulk-seq的數(shù)據(jù)，尋找與某一性狀顯著相關(guān)的單細(xì)胞亞群，然后從表型的角度解釋這些細(xì)胞亞群的生物學(xué)意義。作者開(kāi)發(fā)Scissor的動(dòng)機(jī)是由于目前對(duì)細(xì)胞亞群的分群大多基于scRNA的表達(dá)量進(jìn)行無(wú)監(jiān)督聚類(lèi)，卻鮮有人從表型的角度解釋這些單細(xì)胞亞群。作者認(rèn)為，相同的細(xì)胞亞群可能會(huì)導(dǎo)致相同表型的發(fā)生

附上文獻(xiàn)地址：《Identifying phenotype-associated subpopulations by integrating bulk and single-cell sequencing data》

完整流程

# 加載數(shù)據(jù)
## 單細(xì)胞數(shù)據(jù)
location <- "https://xialab.s3-us-west-2.amazonaws.com/Duanchen/Scissor_data/"
load(url(paste0(location, 'scRNA-seq.RData')))
## bulk-seq data
load(url(paste0(location, 'TCGA_LUAD_exp1.RData')))
## 表型 data
load(url(paste0(location, 'TCGA_LUAD_survival.RData')))

# 預(yù)處理單細(xì)胞數(shù)據(jù)
sc_dataset <- Seurat_preprocessing(sc_dataset, verbose = F)

# 預(yù)處理表型數(shù)據(jù)
phenotype <- bulk_survival[,2:3]
colnames(phenotype) <- c("time", "status")

# 運(yùn)行Scissor
infos1 <- Scissor(bulk_dataset, sc_dataset, phenotype, alpha = 0.05, 
                 family = "cox", Save_file = 'Scissor_LUAD_survival.RData')


Scissor_select <- rep(0, ncol(sc_dataset))
names(Scissor_select) <- colnames(sc_dataset)
## infos1$Scissor_pos代表 β > 0 的細(xì)胞
## infos1$Scissor_neg代表 β < 0 的細(xì)胞
Scissor_select[infos1$Scissor_pos] <- 1
Scissor_select[infos1$Scissor_neg] <- 2
sc_dataset <- AddMetaData(sc_dataset, metadata = Scissor_select, col.name = "scissor")
DimPlot(sc_dataset, reduction = 'umap', group.by = 'scissor', cols = c('grey','indianred1','royalblue'), pt.size = 1.2, order = c(2,1))


如上圖，每一個(gè)點(diǎn)代表一個(gè)細(xì)胞，紅色代表與表型呈正相關(guān)的細(xì)胞；藍(lán)色代表與表型呈負(fù)相關(guān)的細(xì)胞
因此，經(jīng)過(guò)計(jì)算獲得的回歸系數(shù)β表征每一個(gè)細(xì)胞與表型的相關(guān)程度，若回歸系數(shù)β > 0代表該細(xì)胞與某種表型呈現(xiàn)正相關(guān)；若回歸系數(shù)β < 0代表該細(xì)胞與某種表型呈現(xiàn)負(fù)相關(guān)

原理


如上圖所示，作者需要的 input 文件有三種，單細(xì)胞數(shù)據(jù)，表型數(shù)據(jù)（可以說(shuō)分類(lèi)表型數(shù)據(jù)，也可以是連續(xù)型表型數(shù)據(jù)）和bulk-seq的表達(dá)矩陣


第一步，軟件利用分位數(shù)回歸去除了bulk-seq和scRNA的批次效應(yīng)；第二步基于單細(xì)胞數(shù)據(jù)構(gòu)建cell與cell間的similarity network（G）；第三步計(jì)算單細(xì)胞表達(dá)矩陣對(duì)bulk-seq表達(dá)矩陣的皮爾斯相關(guān)系數(shù)，記作S={sij}n×m，n為sample的總數(shù)目，m為細(xì)胞的總數(shù)目；第四步，利用相關(guān)性矩陣S作為決策變量，表型數(shù)據(jù)作為響應(yīng)變量建立回歸關(guān)系，設(shè)回歸系數(shù)為β，在計(jì)算β的過(guò)程中將以及cell與cell間的similarity network（G）的部分信息（利用度矩陣和鄰接矩陣構(gòu)建拉普拉斯矩陣）作為估計(jì)的正則項(xiàng)

那么接下來(lái)的任務(wù)就是估計(jì)出回歸系數(shù)β，轉(zhuǎn)換為最優(yōu)化問(wèn)題：
β等于

其中：

  1. L(拉普拉斯矩陣） = I - D-1/2AD-1/2
  2. 上式的A={aij}m×m代表cell與cell間的similarity network（G）的鄰接矩陣（aij介于0—1之間）
  3. 上式的D={dij}m×m是G的度矩陣，dij=∑j=1aij，其中 j = 1—m
  4. λ 為 penalty

其中：

  1. 設(shè)
    表示sample i 對(duì)每一個(gè)細(xì)胞的皮爾斯相關(guān)系數(shù)，si2表示sample i 與第2個(gè)細(xì)胞的皮爾斯相關(guān)系數(shù)
  2. 若表型數(shù)據(jù)Y是連續(xù)型變量
    l(β)定義為：

    3.若表型數(shù)據(jù)Y是分類(lèi)型變量
    l(β)定義為：

那么計(jì)算出來(lái)的回歸系數(shù)為β越高代表某細(xì)胞亞群與某表型的相關(guān)性比較高，反之比較低

事實(shí)上bulk-seq的sample數(shù)量與表型數(shù)據(jù)的數(shù)量是一致的。而決策變量Si表征每個(gè)細(xì)胞與sample i的相關(guān)性，相關(guān)性高即代表該細(xì)胞與該sample的表達(dá)模式相同，也就是sample i 中這個(gè)細(xì)胞的含量較多（該細(xì)胞含量多才會(huì)使得該細(xì)胞與該sample之間表達(dá)模式相同），因此可以等量代換為某個(gè)細(xì)胞的含量與表型之間的關(guān)系，因此β值為正且越大，則說(shuō)明該細(xì)胞含量對(duì)表型影響呈正相關(guān)且影響大。反之β為負(fù)且越小，則說(shuō)明該細(xì)胞含量對(duì)表型影響呈正負(fù)相關(guān)且影響大

代碼

第一部分：

# single cell
location <- "https://xialab.s3-us-west-2.amazonaws.com/Duanchen/Scissor_data/"
load(url(paste0(location, 'scRNA-seq.RData')))
sc_dataset <- Seurat_preprocessing(sc_dataset, verbose = F)

# bulk seq
load(url(paste0(location, 'TCGA_LUAD_exp1.RData')))
load(url(paste0(location, 'TCGA_LUAD_survival.RData')))

# 提取表型數(shù)據(jù)
all(colnames(bulk_dataset) == bulk_survival$TCGA_patient_barcode)
phenotype <- bulk_survival[,2:3]
colnames(phenotype) <- c("time", "status")

其中bulk-seq的列表示不同sample，行表示不同基因；表型數(shù)據(jù)：

phenotype

第二部分：

# 讀取數(shù)據(jù)
bulk_dataset = bulk_dataset
sc_dataset = sc_dataset
phenotype = phenotype
tag = NULL
alpha = 0.05
cutoff = 0.2
family = "cox"
Save_file = "Scissor_inputs.RData"
Load_file = NULL

# 單細(xì)胞基因和bulk-seq基因取交集
common <- intersect(rownames(bulk_dataset), rownames(sc_dataset))

# 構(gòu)建cell與cell間的similarity network（G）
sc_exprs <- as.matrix(sc_dataset@assays$RNA@data)
network  <- as.matrix(sc_dataset@graphs$RNA_snn)
diag(network) <- 0
network[which(network != 0)] <- 1

# Dataset before quantile normalization.
dataset0 <- cbind(bulk_dataset[common,], sc_exprs[common,])        
# Dataset after  quantile normalization.
dataset1 <- normalize.quantiles(dataset0)                           
rownames(dataset1) <- rownames(dataset0)
colnames(dataset1) <- colnames(dataset0)

# 得到單細(xì)胞表達(dá)矩陣和bulk-seq表達(dá)矩陣
Expression_bulk <- dataset1[,1:ncol(bulk_dataset)]
Expression_cell <- dataset1[,(ncol(bulk_dataset) + 1):ncol(dataset1)]
# 單細(xì)胞矩陣和bulk-seq表達(dá)量矩陣求皮爾斯相關(guān)系數(shù)所構(gòu)成的相關(guān)系數(shù)矩陣，作為決策變量
X <- cor(Expression_bulk, Expression_cell)
quality_check <- quantile(X)

# 定義表型的響應(yīng)變量
Y <- as.matrix(phenotype)

# 計(jì)算回歸系數(shù)β
set.seed(123)
# 初次計(jì)算
fit0 <- APML1(X, Y, family = family, penalty = "Net", alpha = alpha[i], Omega = network, nlambda = 100, nfolds = 10)
# 再次計(jì)算，將fit0$lambda.min作為新模型的 lambda
fit1 <- APML1(X, Y, family = family, penalty = "Net", alpha = alpha[i], Omega = network, lambda = fit0$lambda.min)
# 獲取回歸系數(shù)β
Coefs <- as.numeric(fit1$Beta)
Cell1 <- colnames(X)[which(Coefs > 0)]
Cell2 <- colnames(X)[which(Coefs < 0)]
percentage <- (length(Cell1) + length(Cell2)) / ncol(X)

因此，經(jīng)過(guò)計(jì)算獲得的回歸系數(shù)β表征每一個(gè)細(xì)胞與表型的相關(guān)程度，若回歸系數(shù)β > 0代表該細(xì)胞與某種表型呈現(xiàn)正相關(guān)；若回歸系數(shù)β < 0代表該細(xì)胞與某種表型呈現(xiàn)負(fù)相關(guān)

細(xì)胞水平評(píng)估

此外，作者還提供了一個(gè)估計(jì)細(xì)胞水平的函數(shù)evaluate.cell()，其核心代碼如下：

Scissor_result = infos1

selected_cell <- c(Scissor_result$Scissor_pos, Scissor_result$Scissor_neg)

# Expression_bulk意義參見(jiàn)第二部分代碼
m <- ncol(Expression_bulk)
# selected_cell意義參見(jiàn)第二部分代碼
n <- length(selected_cell)
evaluate_summary <- as.data.frame(matrix(0, n, 11, dimnames = list(selected_cell,
                        c("Mean correlation","Correlation > 0","Correlation < 0","Significant Correlation","Coefficient",
                          "Beta 0%","Beta 25%","Beta 50%","Beta 75%","Beta 100%","Probability of zero"))))

   

cor_test_p <- matrix(0, m, n)
for (i in 1:m){
     Sys.sleep(1 / 100)
     for (j in 1:n){
          #相關(guān)性的顯著性檢驗(yàn)
          cor_test_p[i,j] <- cor.test(Expression_bulk[,i], Expression_cell[,selected_cell[j]])$p.value
    }
}
    
cor_test_FDR <- matrix(p.adjust(as.numeric(cor_test_p), method = "fdr"), m)
for (j in 1:n){
     # X代表單細(xì)胞矩陣和bulk-seq表達(dá)量矩陣求皮爾斯相關(guān)系數(shù)所構(gòu)成的相關(guān)系數(shù)矩陣，參考第二部分
      evaluate_summary[j,1] <- mean(X[,selected_cell[j]])
      evaluate_summary[j,2] <- percent(sum(X[,selected_cell[j]] > 0)/m)
      evaluate_summary[j,3] <- percent(sum(X[,selected_cell[j]] < 0)/m)
      evaluate_summary[j,4] <- percent(sum(cor_test_FDR[,j] < FDR_cutoff)/m)
}

其中：X代表單細(xì)胞矩陣和bulk-seq表達(dá)量矩陣求皮爾斯相關(guān)系數(shù)所構(gòu)成的相關(guān)系數(shù)矩陣



因此結(jié)果的第一列表示的是某一個(gè)細(xì)胞與所有bulk-seq樣本相關(guān)性的均值；第二列是統(tǒng)計(jì)X矩陣相關(guān)系數(shù)大于0的比例；第三列是統(tǒng)計(jì)X矩陣相關(guān)系數(shù)小于0的比例

小tips：關(guān)于APML1函數(shù)

以cox為參數(shù)作為例子

# 讀取數(shù)據(jù)
x=X
y=Y
Omega = network
alpha = alpha[i]
lambda=NULL
nlambda = 100
rlambda=NULL
nfolds = 10
foldid=NULL
inzero=TRUE
wbeta=rep(1,ncol(x))
sgn=rep(1,ncol(x))
isd=FALSE
iysd=FALSE
keep.beta=FALSE
ifast=TRUE
thresh=1e-7
maxit=1e+5
penalty="Net"
wbeta=abs(wbeta)
family = 'cox'

fit=CoxL0(x,y,Omega,alpha,lambda,nlambda,rlambda,nfolds,foldid,inzero,wbeta,sgn,isd,keep.beta,ifast,thresh,maxit)

class(fit)="APML1"
# 最終返回 fit
return(fit)

其中CoxL0函數(shù)的功能是估計(jì)β的算式，是用C++寫(xiě)的

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書(shū)系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容