單細(xì)胞加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析scWGCNA-01

scWGCNA是一個(gè)生物信息學(xué)工作流程糖儡,是R包WGCNA的附加組件仅政,用于在單細(xì)胞或單核RNA-seq數(shù)據(jù)集中進(jìn)行加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析。WGCNA最初是為分析大量基因表達(dá)數(shù)據(jù)集而構(gòu)建的,由于scRNA-seq數(shù)據(jù)的固有稀疏性,vanilla WGCNA在單細(xì)胞數(shù)據(jù)上的性能受到限制。為了解決這個(gè)問題津坑,scWGCNA在運(yùn)行WGCNA分析之前將轉(zhuǎn)錄表達(dá)譜相似的細(xì)胞聚集成pseudo-bulk元細(xì)胞。

scWGCNA代碼比較簡短傲霸,對WGCNA的原有算法的一種改進(jìn),以適應(yīng)單細(xì)胞稀疏矩陣(很多表達(dá)為0);

該包只包含兩個(gè)函數(shù)metacells_by_groupsconstruct_metacells昙啄,construct_metacells()函數(shù)將 Seurat 對象穆役,基于knn算法,將相鄰細(xì)胞構(gòu)建平均的“元細(xì)胞”(averaged metacells)梳凛。

github網(wǎng)站:https://github.com/smorabit/scWGCNA
在線教程:https://smorabit.github.io/tutorials/9_scWGCNA_tutorial

作者在2021年nature genetics發(fā)表的文章Single-nucleus chromatin accessibility and transcriptomic characterization of Alzheimer's disease中應(yīng)用了此算法耿币。

1. 先決條件

運(yùn)行 scWGCNA,需要有一個(gè) Seurat 格式的單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)集韧拒,且已經(jīng)進(jìn)行聚類和降維淹接。

2. scWGCNA安裝

# 要運(yùn)行scWGCNA,需要安裝其他一些R包
install.packages('WGCNA') 
install.packages('igraph') 
install.packages('devtools') 
# install Seurat, check their website for the most up-to-date instructions 
install.packages('Seurat')

devtools::install_github('smorabit/scWGCNA')

library(scWGCNA)

3. scWGCNA實(shí)例教程

官方教程介紹從操作Seurat 對象到構(gòu)建metacells逐步展開分析叛溢。 如果你已經(jīng)有一個(gè)已聚類的數(shù)據(jù)集塑悼,你可以跳過接下來的前兩部分,從“構(gòu)建metacells”部分開始分析楷掉。

3.1 聚類和降維

首先厢蒜,從已發(fā)表的人類大腦健康和疾病的 snRNA-seq 研究中收集了幾個(gè)數(shù)據(jù)集,總計(jì)超過 50 萬個(gè)單核轉(zhuǎn)錄組細(xì)胞烹植。這些數(shù)據(jù)集都滿足scWGCNA 的先決條件斑鸦,并且可以通過許多不同的方式進(jìn)行聚類/降維。由于處理來自多個(gè)來源的數(shù)據(jù)草雕,因此使用online iNMF來構(gòu)建不受批次或數(shù)據(jù)來源影響的降維巷屿,以反映潛在的細(xì)胞異質(zhì)性。

代碼中總共包含9個(gè)文獻(xiàn)數(shù)據(jù)集墩虹,數(shù)據(jù)非常不容易下載攒庵,先放棄實(shí)操。

step1: 將9個(gè)seurat對象merge在一起败晴,生成seurat_obj對象浓冒,進(jìn)行Normalize歸一化和scale處理;
step2: 提取9個(gè)數(shù)據(jù)集的基因表達(dá)矩陣尖坤,生成expression_list稳懒,構(gòu)建Liger對象,Liger是一種去批次慢味,多數(shù)據(jù)集整合方法场梆,跟SCTransform-CCA,Harmony方法類似纯路。對Liger對象進(jìn)行如下操作:1)normalize或油,scale等預(yù)處理;2)執(zhí)行online iNMF驰唬;3) quantile歸一化顶岸。將liger_obj保存為9_datasets_online_liger.rds腔彰。
step3:將iNMF坐標(biāo)轉(zhuǎn)移至seurat_obj,我們之前的seurat_obj降維方式通常只有pca辖佣,umap霹抛,tsne,現(xiàn)在將增加iNMF降維卷谈,另外杯拐,用liger_obj的umap坐標(biāo)替換掉seurat_obj的坐標(biāo)。將修改后的seurat_obj保存為9_datasets_online_seurat.rds世蔗。

library(Seurat) 
library(Matrix) 
library(tidyverse) 
library(rliger) 
data_dir <- "data/" 

# load individual seurat objects: 
seurat_trem2 <- readRDS(file=paste0(data_dir, 'seurat_objects/trem2_nmed_unprocessed.rds')) 
seurat_allen <- readRDS(file=paste0(data_dir, 'seurat_objects/allen_brain_map_2020_unprocessed.rds')) 
seurat_jakel <- readRDS(file=paste0(data_dir, 'seurat_objects/jakel_2019_unprocessed.rds')) 
seurat_grubman <- readRDS(file=paste0(data_dir, 'seurat_objects/grubman_2019_unprocessed.rds')) 
seurat_schirmer <- readRDS(file=paste0(data_dir, 'seurat_objects/schirmer_2019_unprocessed.rds')) 
seurat_velmeshev <- readRDS(file=paste0(data_dir, 'seurat_objects/velmeshev_2019_unprocessed.rds')) 
seurat_swarup_splitseq <- readRDS(file=paste0(data_dir, 'seurat_objects/swarup_AD_splitseq_2019_unprocessed.rds')) 
seurat_swarup <- readRDS(file=paste0(data_dir, 'seurat_objects/swarup_AD_2019_unprocessed.rds')) 
seurat_tsai <- readRDS(file=paste0(data_dir, 'seurat_objects/tsai_AD_2019_unprocessed.rds')) 
# make one big merged Seurat object: 
seurat_obj <- merge(c(seurat_trem2, seurat_allen, seurat_jakel, seurat_grubman, seurat_schirmer, seurat_velmeshev, seurat_swarup_splitseq, seurat_swarup, seurat_tsai)) 
# Normalize and scale data: 
seurat_obj <- NormalizeData(seurat_obj) 
seurat_obj <- ScaleData(seurat_obj) 
# make a list of expression matrices: 
expression_list <- list( 
  'trem2' = GetAssayData(seurat_trem2, slot='counts'), 
  'allen'= GetAssayData(seurat_allen, slot='counts'), 
  'jakel'= GetAssayData(seurat_jakel, slot='counts'), 
  'grubman'= GetAssayData(seurat_grubman, slot='counts'), 
  'schirmer'= GetAssayData(seurat_schirmer, slot='counts'), 
  'velmeshev'= GetAssayData(seurat_velmeshev, slot='counts'), 
  'swarup_splitseq'= GetAssayData(seurat_swarup_splitseq, slot='counts'), 
  'swarup_AD'= GetAssayData(seurat_swarup, slot='counts'), 
  'mathys'= GetAssayData(seurat_tsai, slot='counts') 
) 
# append dataset names to barcodes: 
for(dataset in names(expression_list)){ 
  colnames(expression_list[[dataset]]) <- paste0(dataset, '_', colnames(expression_list[[dataset]])) 
} 
# combined metadata: 
meta_columns <- c('orig.ident', 'nCount_RNA', 'nFeature_RNA', 'Age', 'Sex', 'Condition', 'Condition_specific', 'Batch', 'SampleID', 'DonorID', 'Original_cluster', 'Dataset', 'Tissue', 'Technology' ) 
seurat_meta <- rbind( 
  seurat_trem2@meta.data[,meta_columns], 
  seurat_allen@meta.data[,meta_columns], 
  seurat_jakel@meta.data[,meta_columns], 
  seurat_grubman@meta.data[,meta_columns], 
  seurat_schirmer@meta.data[,meta_columns], 
  seurat_velmeshev@meta.data[,meta_columns], 
  seurat_swarup_splitseq@meta.data[,meta_columns], 
  seurat_swarup@meta.data[,meta_columns], 
  seurat_tsai@meta.data[,meta_columns] 
) 
rownames(seurat_meta) <- paste0(seurat_meta$Dataset, '_', rownames(seurat_meta)) 
# remove individual seurat objects to save space: 
rm(seurat_trem2, seurat_allen, seurat_jakel, seurat_grubman, seurat_schirmer, seurat_velmeshev, seurat_swarup, seurat_swarup_splitseq, seurat_tsai); gc(); 
# create liger object: 
liger_obj <- createLiger(expression_list) 
# pre-processing 
liger_obj <- normalize(liger_obj) 
liger_obj <- selectGenes(liger_obj,) 
liger_obj <- scaleNotCenter(liger_obj) 
# perform online iNMF 
liger_obj = online_iNMF(liger_obj, k=50, max.epochs=5) 
# quantile normalization 
liger_obj  = quantile_norm(liger_obj) 
liger_obj  = runUMAP(liger_obj) 
saveRDS(liger_obj, file=paste0(data_dir, '9_datasets_online_liger.rds')) 
# transfer iNMF matrix to seurat obj: 
seurat_obj@reductions$iNMF <- CreateDimReducObject( 
    loadings=t(liger_obj@W), 
    embeddings=liger_obj@H.norm[colnames(seurat_obj),], 
    key="iNMF_", 
    assay="RNA" 
) 
# add UMAP: 
seurat_obj@reductions$UMAP <- CreateDimReducObject( 
  embeddings = liger_obj@tsne.coords[colnames(seurat_obj),], 
  key='umap', 
  assay='RNA' 
) 
# add clusters from liger 
seurat_obj@meta.data$liger_clusters <- liger_obj@clusters 
saveRDS(seurat_obj, file=paste0(data_dir, '9_datasets_online_seurat.rds'))

問題1:liger_obj進(jìn)行runUMAP降維端逼,但是在后面的代碼中,出現(xiàn)liger_obj@tsne.coords污淋,有點(diǎn)看不懂顶滩,是不是誤寫?
seurat_obj@reductions$UMAP <- CreateDimReducObject( 
  embeddings = liger_obj@umap.coords[colnames(seurat_obj),],  
  key='umap', 
  assay='RNA' 
)

找同行確認(rèn)了下芙沥,這應(yīng)該是作者誤寫诲祸。


image.png
問題2:seurat_meta整合9個(gè)樣本的meta.data信息(orig.ident,nCount_RNA等)而昨,但是代碼中沒有將seurat_meta賦值給liger_obj@meta.data或seurat_obj@meta.data救氯, 也沒有保存?
liger_obj@meta.data <- seurat_meta
seurat_obj@meta.data <- seurat_meta

下圖是用online iNMF整合方法得到的UMAP聚類圖

image
3.2 獲取感興趣的特定細(xì)胞類型子集

為了識別給定細(xì)胞類型內(nèi)的共表達(dá)模塊歌憨,作者僅使用數(shù)據(jù)的一個(gè)子集重復(fù)上述分析着憨。在這里,我們對感興趣的少突膠質(zhì)細(xì)胞譜系細(xì)胞進(jìn)行這種分析务嫡,包括少突膠質(zhì)祖細(xì)胞和成熟的少突膠質(zhì)細(xì)胞甲抖。
step1:從seurat_obj提取細(xì)胞類型為ODC,OPC的子集seurat_odc心铃;同上述操作准谚,對seurat_odc數(shù)據(jù)集進(jìn)行iNMF整合,降維去扣,將liger_obj保存為ODC_9_datasets_online_liger.rds柱衔。
step2:將iNMF坐標(biāo)轉(zhuǎn)移至seurat_odc,命名為ctiNMF愉棱,這將新增ctiNMF降維方式(類似于pca)唆铐,seurat_odc對象scale處理后,基于ctiNMF進(jìn)行umap降維奔滑,將生成的seurat_odc保存為9_datasets_ODC_seurat.rds艾岂。

# subset by neuronal clusters
seurat_odc <- subset(seurat_obj, celltype %in% c('ODC', 'OPC'))

# iNMF for just ODCs:
expr_matrix <- GetAssayData(seurat_odc, slot='counts')

expression_list <- list(
  'zhou' = expr_matrix[,colnames(seurat_odc)[seurat_odc$Dataset == 'zhou']],
  'allen'= expr_matrix[,colnames(seurat_odc)[seurat_odc$Dataset == 'allen']],
  'jakel'= expr_matrix[,colnames(seurat_odc)[seurat_odc$Dataset == 'jakel']],
  'grubman'= expr_matrix[,colnames(seurat_odc)[seurat_odc$Dataset == 'grubman']],
  'schirmer'= expr_matrix[,colnames(seurat_odc)[seurat_odc$Dataset == 'schirmer']],
  'velmeshev'= expr_matrix[,colnames(seurat_odc)[seurat_odc$Dataset == 'velmeshev']],
  'swarup_splitseq'= expr_matrix[,colnames(seurat_odc)[seurat_odc$Dataset == 'swarup_splitseq']],
  'swarup_AD'= expr_matrix[,colnames(seurat_odc)[seurat_odc$Dataset == 'swarup_AD_2020']],
  'mathys'= expr_matrix[,colnames(seurat_odc)[seurat_odc$Dataset == 'mathys']]
)

seurat_meta <- rbind(
  subset(seurat_odc@meta.data, Dataset=='zhou'),
  subset(seurat_odc@meta.data, Dataset=='allen'),
  subset(seurat_odc@meta.data, Dataset=='jakel'),
  subset(seurat_odc@meta.data, Dataset=='grubman'),
  subset(seurat_odc@meta.data, Dataset=='schirmer'),
  subset(seurat_odc@meta.data, Dataset=='velmeshev'),
  subset(seurat_odc@meta.data, Dataset=='swarup_splitseq'),
  subset(seurat_odc@meta.data, Dataset=='swarup_AD_2020'),
  subset(seurat_odc@meta.data, Dataset=='mathys')
)

# create liger object:
liger_obj <- createLiger(expression_list)

liger_obj <- normalize(liger_obj)
pdf("liger_variable_genes.pdf", width=8, height=8)
liger_obj <- selectGenes(
  liger_obj,
  var.thresh =c(0.05, 0.15, 0.025, 0.025, 0.075, 0.075, 0.05, 0.20, 0.025),
  do.plot=T
)
dev.off()
liger_obj@var.genes %>% length
liger_obj <- scaleNotCenter(liger_obj)

# perform online iNMF
liger_obj = online_iNMF(liger_obj, k=20, max.epochs=5)

# quantile normalization
liger_obj  = quantile_norm(liger_obj)
liger_obj  = runUMAP(liger_obj)

saveRDS(liger_obj, file='ODC_9_datasets_online_liger.rds')

# transfer iNMF matrix to seurat obj:
seurat_odc@reductions$ctiNMF <- CreateDimReducObject(
    loadings=t(liger_obj@W),
    embeddings=liger_obj@H.norm[colnames(seurat_odc),],
    key="ctiNMF_",
    assay="RNA"
  )
VariableFeatures(seurat_odc) <- liger_obj@var.genes

# scale expression data:
seurat_odc <- ScaleData(seurat_odc, features=VariableFeatures(seurat_odc))

# UMAP + clustering:
seurat_odc <- RunUMAP(seurat_odc, reduction='ctiNMF', dims=1:20)

# save the results
saveRDS(seurat_odc, file ='9_datasets_ODC_seurat.rds')

# optional: format for use in Scanpy
library(SeuratDisk)
SaveH5Seurat(cur_save, filename = "9_datasets_ODC_seurat.h5Seurat", overwrite=TRUE)
Convert("9_datasets_ODC_seurat.h5Seurat", dest='h5ad', overwrite=TRUE)

下面可以看到少突膠質(zhì)細(xì)胞譜系 細(xì)胞的UMAP聚類圖,作者根據(jù)已知標(biāo)記基因?qū)⑵浞譃?4 個(gè)主要cluster朋其。


image
3.3 構(gòu)建metacells

現(xiàn)在我們準(zhǔn)備構(gòu)建metacells王浴。為了節(jié)省內(nèi)存脆炎,我們只使用高可變基因來構(gòu)建metacells并執(zhí)行下游 WGCNA分析。 可以使用 HVG基因叼耙,也可以用所有基因腕窥,下游結(jié)果將基因數(shù)目的差異而有所不同粒没。因此筛婉,你可以重新運(yùn)行此步驟以嘗試不同的參數(shù)。

step1:這個(gè)數(shù)據(jù)集有幾種不同的疾病處理?xiàng)l件癞松,但是爽撒,我們只想聚合相同條件下的細(xì)胞,即odc_group和Condition屬性一致的細(xì)胞群响蓉。 因此硕勿,我們在 Seurat 對象中創(chuàng)建了一個(gè)新的metadata列metacell_group ,它將 ODC 亞群(odc_group)與處理?xiàng)l件(Condition)組合在一起枫甲。

step2:接下來源武,我們將來自同一 ODC 亞群和相同疾病處理組的細(xì)胞群構(gòu)建metacells。construct_metacells()函數(shù)有一個(gè)關(guān)鍵參數(shù)是 k想幻,有多少個(gè)細(xì)胞合并成一個(gè)metacells粱栖。 這里我們使用了一個(gè)較大數(shù)值 100,這僅僅是因?yàn)槲覀冞@里的數(shù)據(jù)集非常大脏毯。 你可以根據(jù)數(shù)據(jù)集中的細(xì)胞數(shù)來調(diào)整 k 值大小闹究。我只是通過 for 循環(huán)單線程運(yùn)行它,當(dāng)然食店,您可以更快速渣淤,例如,如果您要將每個(gè)調(diào)用construct_metacells函數(shù)的程序投入到HPC計(jì)算節(jié)點(diǎn)吉嫩,實(shí)現(xiàn)并行運(yùn)算价认。該construct_metacells() 函數(shù)輸出一個(gè)新的Seurat 對象,包含聚合后的基因表達(dá)譜自娩。

step3:將每個(gè)metacell_group構(gòu)建好的metacells(Seurat 對象)進(jìn)行merge用踩,scale縮放,Harmony整合椒功,umap聚類捶箱。

# load scWGCNA package
library(scWGCNA)

# construct metacells by ODC group, condition
seurat_odc$metacell_group <- paste0(
  as.character(seurat_odc$odc_group), '_',
  as.character(seurat_odc$Condition)
)

genes.keep <- VariableFeatures(seurat_odc)

# loop through each group and construct metacells
seurat_list <- list()
for(group in unique(seurat_odc$metacell_group)){
  print(group)
  cur_seurat <- subset(seurat_odc, metacell_group == group)
  cur_seurat <- cur_seurat[genes.keep,]
  cur_metacell_seurat <- scWGCNA::construct_metacells(
    cur_seurat, name=group,
    k=100, reduction='umap',
    assay='RNA', slot='data'
  )
  cur_metacell_seurat$Condition <- as.character(unique(cur_seurat$Condition))
  cur_metacell_seurat$odc_group <- as.character(unique(cur_seurat$odc_group))
  seurat_list[[group]] <- cur_metacell_seurat
}

# merge all of the metacells objects
metacell_seurat <- merge(seurat_list[[1]], seurat_list[2:length(seurat_list)])
saveRDS(metacell_seurat, file='data/metacell_seurat.rds')

Optionally, we can run a dimensionality reduction on the metacell Seurat object to check if some cellular heterogeneity has been retained. The differences in the transcriptional profiles of the different disease groups are quite prominent in the UMAP space, so I am harmonizing the data on the basis of condition merely for visualization purposes.

metacell_subset <- ScaleData(metacell_subset, features = rownames(metacell_subset))
metacell_subset <- RunPCA(metacell_subset, features=rownames(metacell_subset))
metacell_subset <- RunHarmony(metacell_subset, group.by='Condition', dims=1:15)
metacell_subset <- RunUMAP(metacell_subset, reduction='harmony', dims=1:15)

pdf(paste0(fig_dir, "metacell_umap_group.pdf"), width=6, height=6)
DimPlot(metacell_subset, group.by='odc_group', reduction='umap', label=TRUE) + umap_theme
dev.off()

由下圖可以看出,通過construct_metacells處理后的seurat對象在umap圖中仍有很好地分離动漾。

image.png

前面代碼最關(guān)鍵的一步是construct_metacells丁屎,將cur_seurat轉(zhuǎn)成cur_metacell_seurat ,我們看下源碼旱眯,這步做了什么處理晨川。

construct_metacells <- function(seurat_obj, name='agg', k=50, reduction='umap', assay='RNA', slot='data'){
  reduced_coordinates <- as.data.frame(seurat_obj@reductions[[reduction]]@cell.embeddings)
  nn_map <- FNN::knn.index(reduced_coordinates, k = (k - 1))
  row.names(nn_map) <- row.names(reduced_coordinates)
  nn_map <- cbind(nn_map, seq_len(nrow(nn_map)))
  good_choices <- seq_len(nrow(nn_map))
  choice <- sample(seq_len(length(good_choices)), size = 1,
      replace = FALSE)
  chosen <- good_choices[choice]
  good_choices <- good_choices[good_choices != good_choices[choice]]
  it <- 0
  k2 <- k * 2
  get_shared <- function(other, this_choice) {
      k2 - length(union(cell_sample[other, ], this_choice))
  }
  while (length(good_choices) > 0 & it < 5000) {
      it <- it + 1
      choice <- sample(seq_len(length(good_choices)), size = 1,
          replace = FALSE)
      new_chosen <- c(chosen, good_choices[choice])
      good_choices <- good_choices[good_choices != good_choices[choice]]
      cell_sample <- nn_map[new_chosen, ]
      others <- seq_len(nrow(cell_sample) - 1)
      this_choice <- cell_sample[nrow(cell_sample), ]
      shared <- sapply(others, get_shared, this_choice = this_choice)

      if (max(shared) < 0.9 * k) {
          chosen <- new_chosen
      }
  }

  cell_sample <- nn_map[chosen, ]
  combs <- combn(nrow(cell_sample), 2)
  shared <- apply(combs, 2, function(x) {
      k2 - length(unique(as.vector(cell_sample[x, ])))
  })

  message(paste0("Overlap QC metrics:\nCells per bin: ",
      k, "\nMaximum shared cells bin-bin: ", max(shared),
      "\nMean shared cells bin-bin: ", mean(shared), "\nMedian shared cells bin-bin: ",
      median(shared)))
  if (mean(shared)/k > 0.1)
      warning("On average, more than 10% of cells are shared between paired bins.")

  exprs_old <- GetAssayData(seurat_obj, assay=assay, slot=slot)

  mask <- sapply(seq_len(nrow(cell_sample)), function(x) seq_len(ncol(exprs_old)) %in%
      cell_sample[x, , drop = FALSE])
  mask <- Matrix::Matrix(mask)
  new_exprs <- (exprs_old %*% mask) / k
  colnames(new_exprs) <- paste0(name, '_', 1:ncol(new_exprs))
  rownames(cell_sample) <- paste0(name, '_', 1:ncol(new_exprs))
  colnames(cell_sample) <- paste0('knn_', 1:ncol(cell_sample))

  # make seurat obj:
  seurat_aggr <- CreateSeuratObject(
    counts = new_exprs
  )
  seurat_aggr
}
代碼解析:

1.假設(shè)seurat_obj對象有5000個(gè)細(xì)胞证九,10000個(gè)基因,即5000*10000的基因表達(dá)矩陣共虑。提取seurat_obj對象的umap坐標(biāo)信息(umap1,umap2)愧怜,生成reduced_coordinates(5000*2矩陣)。

2.利用FNN快速最近鄰搜索算法妈拌,為每個(gè)細(xì)胞尋找99個(gè)最近鄰細(xì)胞拥坛,nn_map為5000*99矩陣;第1-99列均為最臨近細(xì)胞的次序號尘分;



添加第100列猜惋,為當(dāng)前細(xì)胞自身的次序號;nn_map為5000*100矩陣培愁;99個(gè)鄰近細(xì)胞加上自身細(xì)胞著摔,共100個(gè)。



3.隨機(jī)從5000個(gè)細(xì)胞中抽取一個(gè)細(xì)胞定续,如編號為2000的細(xì)胞N2000谍咆,然后從剩下的4999個(gè)細(xì)胞中再隨機(jī)抽取一個(gè)編號N100細(xì)胞,計(jì)算N2000和N100的最臨近細(xì)胞群數(shù)目(200個(gè)細(xì)胞編號)私股,對著200個(gè)細(xì)胞編號摹察,取其并集,如滿足max(shared) < 0.9 * k條件庇茫,也就是說港粱,如果這兩細(xì)胞越相似,他們的共享鄰居應(yīng)該越多旦签,把N2000和N100細(xì)胞合并成chosen組查坪,反之亦然。接著從剩下的細(xì)胞群抽取宁炫,與上一步的chosen 組比較偿曙,看相似程度,迭代5000次羔巢,把相似的細(xì)胞編號放入chosen組望忆。
4.通過上一步的找鄰居策略,chosen組為2000個(gè)竿秆,5000個(gè)總細(xì)胞群启摄,有2000個(gè)最相似的最近鄰居群, cell_sample為2000*100矩陣幽钢。

5.2000個(gè)細(xì)胞兩兩比較歉备,查看共有鄰居數(shù),統(tǒng)計(jì)max(shared)匪燕,median(shared)蕾羊,mean(shared)喧笔。
6.提取5000*10000的基因表達(dá)矩陣exprs_old ,(exprs_old %*% mask) / k的計(jì)算是2000個(gè)細(xì)胞鄰近細(xì)胞群的平均表達(dá)量龟再,轉(zhuǎn)換成2000*10000书闸,即5000個(gè)細(xì)胞由2000個(gè)細(xì)胞所代表,數(shù)據(jù)集也實(shí)現(xiàn)了降維利凑。

construct_metacells算法小結(jié):

1.seurat_obj對象有5000個(gè)細(xì)胞浆劲,10000個(gè)基因,是5000*10000的基因表達(dá)矩陣截碴;設(shè)置參數(shù)k=100梳侨,通過KNN算法查找100個(gè)最臨近鄰居群(包括自身細(xì)胞)蛉威,如A細(xì)胞日丹,找99個(gè)鄰居,加上自己蚯嫌,共100個(gè)哲虾。
2.指定一套篩選相似細(xì)胞策略,將細(xì)胞的鄰居群相互比較择示,保留相似的細(xì)胞束凑,5000個(gè)細(xì)胞保留了最相似的2000個(gè)細(xì)胞,用這2000個(gè)細(xì)胞代表整體細(xì)胞栅盲。
3.之前是5000*10000的基因表達(dá)矩陣汪诉,新的基因表達(dá)矩陣是:2000個(gè)細(xì)胞的鄰居群(100個(gè)細(xì)胞)的平均基因表達(dá)量,如A細(xì)胞的gene1表達(dá)量為100個(gè)鄰居群在gene1的表達(dá)平均值谈秫,因此新矩陣為2000*10000扒寄,剔除了噪音數(shù)據(jù),也實(shí)現(xiàn)了數(shù)據(jù)降維拟烫。這應(yīng)該也是處理單細(xì)胞表達(dá)譜的一種數(shù)據(jù)處理策略该编。

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末,一起剝皮案震驚了整個(gè)濱河市硕淑,隨后出現(xiàn)的幾起案子课竣,更是在濱河造成了極大的恐慌,老刑警劉巖置媳,帶你破解...
    沈念sama閱讀 222,104評論 6 515
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件于樟,死亡現(xiàn)場離奇詭異,居然都是意外死亡拇囊,警方通過查閱死者的電腦和手機(jī)迂曲,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 94,816評論 3 399
  • 文/潘曉璐 我一進(jìn)店門,熙熙樓的掌柜王于貴愁眉苦臉地迎上來寂拆,“玉大人奢米,你說我怎么就攤上這事抓韩。” “怎么了鬓长?”我有些...
    開封第一講書人閱讀 168,697評論 0 360
  • 文/不壞的土叔 我叫張陵谒拴,是天一觀的道長。 經(jīng)常有香客問我涉波,道長英上,這世上最難降的妖魔是什么? 我笑而不...
    開封第一講書人閱讀 59,836評論 1 298
  • 正文 為了忘掉前任啤覆,我火速辦了婚禮苍日,結(jié)果婚禮上,老公的妹妹穿的比我還像新娘窗声。我一直安慰自己相恃,他們只是感情好,可當(dāng)我...
    茶點(diǎn)故事閱讀 68,851評論 6 397
  • 文/花漫 我一把揭開白布笨觅。 她就那樣靜靜地躺著拦耐,像睡著了一般。 火紅的嫁衣襯著肌膚如雪见剩。 梳的紋絲不亂的頭發(fā)上杀糯,一...
    開封第一講書人閱讀 52,441評論 1 310
  • 那天,我揣著相機(jī)與錄音苍苞,去河邊找鬼固翰。 笑死,一個(gè)胖子當(dāng)著我的面吹牛羹呵,可吹牛的內(nèi)容都是我干的骂际。 我是一名探鬼主播,決...
    沈念sama閱讀 40,992評論 3 421
  • 文/蒼蘭香墨 我猛地睜開眼担巩,長吁一口氣:“原來是場噩夢啊……” “哼方援!你這毒婦竟也來了?” 一聲冷哼從身側(cè)響起涛癌,我...
    開封第一講書人閱讀 39,899評論 0 276
  • 序言:老撾萬榮一對情侶失蹤犯戏,失蹤者是張志新(化名)和其女友劉穎,沒想到半個(gè)月后拳话,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體先匪,經(jīng)...
    沈念sama閱讀 46,457評論 1 318
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 38,529評論 3 341
  • 正文 我和宋清朗相戀三年弃衍,在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了呀非。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
    茶點(diǎn)故事閱讀 40,664評論 1 352
  • 序言:一個(gè)原本活蹦亂跳的男人離奇死亡,死狀恐怖岸裙,靈堂內(nèi)的尸體忽然破棺而出猖败,到底是詐尸還是另有隱情,我是刑警寧澤降允,帶...
    沈念sama閱讀 36,346評論 5 350
  • 正文 年R本政府宣布恩闻,位于F島的核電站,受9級特大地震影響剧董,放射性物質(zhì)發(fā)生泄漏幢尚。R本人自食惡果不足惜,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 42,025評論 3 334
  • 文/蒙蒙 一翅楼、第九天 我趴在偏房一處隱蔽的房頂上張望尉剩。 院中可真熱鬧,春花似錦毅臊、人聲如沸理茎。這莊子的主人今日做“春日...
    開封第一講書人閱讀 32,511評論 0 24
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽功蜓。三九已至,卻和暖如春宠蚂,著一層夾襖步出監(jiān)牢的瞬間,已是汗流浹背童社。 一陣腳步聲響...
    開封第一講書人閱讀 33,611評論 1 272
  • 我被黑心中介騙來泰國打工求厕, 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留,地道東北人扰楼。 一個(gè)月前我還...
    沈念sama閱讀 49,081評論 3 377
  • 正文 我出身青樓呀癣,卻偏偏與公主長得像,于是被迫代替她去往敵國和親弦赖。 傳聞我的和親對象是個(gè)殘疾皇子项栏,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 45,675評論 2 359

推薦閱讀更多精彩內(nèi)容