使用MuSiC以及MuSiC2来根据单细胞转录组结果推断bulk转录组细胞比例

# BiocManager::install('TOAST')
# devtools::install_github( repo = "crhisto/xbioc")

# install devtools if necessary
if (!"devtools" %in% rownames(installed.packages())) {
  install.packages('devtools')
}
# install MuSiC package
if (!"MuSiC" %in% rownames(installed.packages())) {
  devtools::install_github('xuranw/MuSiC')
}
# install the MuSiC2 package
if (!"MuSiC2" %in% rownames(installed.packages())) {
  devtools::install_github('Jiaxin-Fan/MuSiC2')
}

## Bulk RNA-seq data
#benchmark.eset = readRDS("https://xuranw.github.io/MuSiC/data/bulk-eset.rds")
benchmark.eset = readRDS("bulk-eset.rds")
benchmark.eset 
table(benchmark.eset$group)
bulk.control.mtx = exprs(benchmark.eset)[, benchmark.eset$group == 'healthy']
bulk.case.mtx = exprs(benchmark.eset)[, benchmark.eset$group == 't2d']
bulk.case.mtx[1:4,1:4]

> benchmark.eset 
ExpressionSet (storageMode: lockedEnvironment)
assayData: 10000 features, 200 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: 1 10 ... 200 (200 total)
  varLabels: sampleID group
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation:  

> table(benchmark.eset$group)

healthy     t2d 
    100     100  
    
> bulk.case.mtx[1:4,1:4]
            101    102    103    104
SFT2D1    16791   8362  17888  12894
TRAPPC13   8583   7823   6674   6626
C9orf89    6745   4659   8482   4790
RPS6     218932 142579 190503 160069

## scRNA-seq data
#seger.sce = readRDS("https://xuranw.github.io/MuSiC/data/EMATBsce_healthy.rds")
seger.sce = readRDS("EMTABesethealthy.rds")
seger.sce
# sc.sce  SingleCellExperiment for single cell data
tail(sort(table(seger.sce$cellType)))
# 需要熟练掌握它们的对象，：一些单细胞转录组R包的对象 
# https://mp.weixin.qq.com/s/lpoHhZqi-_ASUaIfpnX96w
seger.sce <- SingleCellExperiment(
  assays = list(counts = exprs(seger.sce) ), 
  colData =   seger.sce@phenoData@data
)
colData(seger.sce)

> seger.sce
ExpressionSet (storageMode: lockedEnvironment)
assayData: 25453 features, 1097 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: AZ_A10 AZ_A11 ... HP1509101_P9 (1097 total)
  varLabels: sampleID SubjectName cellTypeID cellType
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation:  

> tail(sort(table(seger.sce$cellType)))

 delta  gamma acinar ductal   beta  alpha 
    59     75    112    135    171    443 
    
> seger.sce
class: SingleCellExperiment 
dim: 25453 1097 
metadata(0):
assays(1): counts
rownames(25453): SGIP1 AZIN2 ... KIR2DL2 KIR2DS3
rowData names(0):
colnames(1097): AZ_A10 AZ_A11 ... HP1509101_P8 HP1509101_P9
colData names(4): sampleID SubjectName cellTypeID cellType
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):

# MuSiC
bulk.mtx = exprs(benchmark.eset)
prop_music=music_prop(bulk.mtx = bulk.mtx, sc.sce = seger.sce, 
                      clusters = 'cellType', samples = 'sampleID',
                      select.ct = c('acinar','alpha', 'beta', 'delta', 'ductal','gamma'),
                      verbose = F)$Est.prop.weighted

prop_all = cbind('proportion'=c(prop_music),
                 'celltype'=rep(colnames(prop_music), 
                                each=nrow(prop_music)), 
                 'sampleID'=rep(rownames(prop_music),times=ncol(prop_music)),
                 'Method'='MuSiC')
prop_all=as.data.frame(prop_all)
prop_all$proportion=as.numeric(as.character(prop_all$proportion))

# music2 deconvolution
set.seed(1234)
library(scater)
est = music2_prop_t_statistics(bulk.control.mtx = bulk.control.mtx, 
                               bulk.case.mtx = bulk.case.mtx, 
                               sc.sce = seger.sce, 
                               clusters = 'cellType', 
                               samples = 'sampleID', 
                               select.ct = c('acinar','alpha','beta','delta','ductal','gamma'), 
                               n_resample=20, sample_prop=0.5,cutoff_c=0.05,cutoff_r=0.01)

est.prop = est$Est.prop

# plot estimated cell type proportions
prop_all = cbind('proportion'=c(est.prop), 'sampleID'=rep(rownames(est.prop),times=ncol(est.prop)), 'celltype'=rep(colnames(est.prop), each=nrow(est.prop)))
prop_all = as.data.frame(prop_all)
prop_all$proportion = as.numeric(as.character(prop_all$proportion))
prop_all$group = ifelse(prop_all$sampleID %in% seq(from=1, to=100, by=1), 'Healthy', 'T2D')