我導(dǎo)再也不用擔(dān)心我不認(rèn)識marker啦

我們在進(jìn)行單細(xì)胞測序的時候，通常情況下是通過高變genes來辨別細(xì)胞類型（于是一大堆不認(rèn)識的），除了免疫細(xì)胞可能容易分析出來，其他的細(xì)胞我是兩眼一抹黑。雖然具有single cell marker基因庫，但在判斷細(xì)胞亞型的時候有時并不是那么solid，有一些原則上比較特異的gene不知道為什么在cell marker上出現(xiàn)在不同的細(xì)胞上，比如S100A8,S100A9，在cell marker基因庫中多為中性粒細(xì)胞的marker，但結(jié)合生理學(xué)意義，其實(shí)在單細(xì)胞中由于中性粒細(xì)胞本身含量較少且較為脆弱等原因?qū)е虏东@中性粒細(xì)胞是極為困難的。此時如果貿(mào)然下結(jié)論為中性粒細(xì)胞其實(shí)不利于后期的分析。

單細(xì)胞分群后，怎么找到Marker基因定義每一類群？

celaref R包通過與已知細(xì)胞類型的參考數(shù)據(jù)集的相似度進(jìn)行比較。輸入每個細(xì)胞中每個基因的reads counts數(shù)（gene-cell matrix）和每個細(xì)胞所屬的簇（cluster）信息，和每個查詢組中最明顯富集的基因的參考樣本比較，通過排名來匹配細(xì)胞類型。

Workflow

下面是典型的celaref工作流程，根據(jù)與參考數(shù)據(jù)集的轉(zhuǎn)錄組相似性來表征查詢數(shù)據(jù)集的細(xì)胞簇。

數(shù)據(jù)及其格式

gene-cell matrix

每個細(xì)胞所屬的cluster信息

輸入數(shù)據(jù)之后利用MAST包進(jìn)行基因表達(dá)差異分析，每個基因被指定為從0（最富集）到1（最不存在）的rescaled rank。

比較查詢數(shù)據(jù)和參考數(shù)據(jù)

得到每個查詢細(xì)胞簇的Up基因列表 — 在該簇中具有顯著更高表達(dá)的基因。在每個參考細(xì)胞簇的基因排名中查找這些基因，比較并繪制相似性。

輸出結(jié)果

通常，查詢數(shù)據(jù)中的每個細(xì)胞簇都針對參考數(shù)據(jù)（X軸）中的所有內(nèi)容繪制。刻度線表示up基因，并且**中位基因（middle?generank）**顯示為粗條。頂部附近的偏差分布表示各組的相似性 - 基本上相同的基因代表它們各自樣品中的cluster。也就是說查詢組聚類分析后的代表基因如果和reference具有一定的相似性，則可以通過其相同的基因代表與其對應(yīng)，也就是細(xì)胞類型的對應(yīng)。

R語言 - 箱線圖（小提琴圖、抖動圖、區(qū)域散點(diǎn)圖）

為clusters分配標(biāo)簽

通過上圖第二列可以發(fā)現(xiàn)其實(shí)可能存在兩種以上的標(biāo)簽，這時候就需要進(jìn)行綜合判斷了。

實(shí)例分析

# Installing BiocManager if necessary: #通過BiocManager進(jìn)行R包加載 # install.packages("BiocManager") BiocManager::install("celaref")

我們以系統(tǒng)包中的dataset為例進(jìn)行演示。

假設(shè)有一個新的scRNAseq數(shù)據(jù)集（查詢），其聚類已經(jīng)聚集成4組：組1—4。但是我們還不知道哪個組對應(yīng)于哪種細(xì)胞類型。

現(xiàn)在有一個相同組織類型的舊數(shù)據(jù)集（參考），其他人已經(jīng)確定了細(xì)胞類型：Weird subtype,?Exciting,?Mystery cell typeand?Dunno。

library(celaref) # Paths to data files. counts_filepath.query <- system.file("extdata", "sim_query_counts.tab", package = "celaref") cell_info_filepath.query <- system.file("extdata", "sim_query_cell_info.tab", package = "celaref") counts_filepath.ref <- system.file("extdata", "sim_ref_counts.tab", package = "celaref") cell_info_filepath.ref <- system.file("extdata", "sim_ref_cell_info.tab", package = "celaref") # 加載數(shù)據(jù) toy_ref_se <- load_se_from_files(counts_file=counts_filepath.ref, cell_info_file=cell_info_filepath.ref) toy_query_se <- load_se_from_files(counts_file=counts_filepath.query, cell_info_file=cell_info_filepath.query) head(de_table.demo_ref) # 參考數(shù)據(jù)格式

head(de_table.toy_query ) # 查詢數(shù)據(jù)格式

# 過濾 toy_ref_se <- trim_small_groups_and_low_expression_genes(toy_ref_se) toy_query_se <- trim_small_groups_and_low_expression_genes(toy_query_se) # 分別得到各自的差異表達(dá)基因 de_table.toy_ref <- contrast_each_group_to_the_rest(toy_ref_se, dataset_name="ref") de_table.toy_query <- contrast_each_group_to_the_rest(toy_query_se, dataset_name="query")

以上兩步的運(yùn)行均需要較長的時間。

# 展示查詢組top基因在參考組中的分布 make_ranking_violin_plot(de_table.test=de_table.toy_query, de_table.ref=de_table.toy_ref)

# 獲得各個group的標(biāo)簽 make_ref_similarity_names(de_table.toy_query, de_table.toy_ref)

通過以上方法我們基本可以辨別細(xì)胞類型，其中我們可以看到部分group并沒有給出細(xì)胞類型（如group1,2在reciprocal_matches中沒有找到匹配類型），有的可能會出現(xiàn)多個類型，需要進(jìn)一步判斷。

我們再以10X數(shù)據(jù)為例進(jìn)行演示|PBMCs - 10X vs Microarray Reference

10X genomic****s有幾個數(shù)據(jù)集可供從他們的網(wǎng)站下載，比如pbmc4k datasets（https://support.10xgenomics.com/single-cell-gene-expression/datasets/2.1.0/pbmc4k），其中包含來自健康個體的外周血細(xì)胞(PBMC)數(shù)據(jù)，下面示例是篩選后的文件Gene/cell matrix(filtered)和聚類文件Clustering analysis。

查詢數(shù)據(jù)：10X query dataset

library(celaref) # 首先設(shè)置路徑并且加載數(shù)據(jù)，我們以kmeans_7_clusters 的聚類結(jié)果來進(jìn)行細(xì)胞定義 datasets_dir <- "~/celaref_extra_vignette_data/datasets"dataset_se.10X_pbmc4k_k7 <- load_dataset_10Xdata(dataset_path = file.path(datasets_dir,'10X_pbmc4k'),dataset_genome = "GRCh38",clustering_set = "kmeans_7_clusters",id_to_use = "GeneSymbol") dataset_se.10X_pbmc4k_k7 <- trim_small_groups_and_low_expression_genes(dataset_se.10X_pbmc4k_k7) # 進(jìn)行部分?jǐn)?shù)據(jù)篩選并轉(zhuǎn)換格式 ？trim_small_groups_and_low_expression_genes # 查詢函數(shù)幫助

de_table.10X_pbmc4k_k7 <- contrast_each_group_to_the_rest(dataset_se.10X_pbmc4k_k7, dataset_name="10X_pbmc4k_k7", num_cores=7) ##進(jìn)行差異比較

參考數(shù)據(jù)：microarray dataset

查詢數(shù)據(jù)有了，那么參考數(shù)據(jù)（reference）呢？

Watkins?等人在2009年的一篇文獻(xiàn)已經(jīng)發(fā)表過合適的? PBMC reference (a HaemAtlas)，這里使用的數(shù)據(jù)是從haemosphere網(wǎng)站（Graaf et al.2016，https://www.haemosphere.org/datasets/show）下載得到。

因為這是微陣列數(shù)據(jù)，所以還需要一些處理步驟：

• Logged, normalised expression values. Any low expression or poor quality ?measurements should have already been removed.

• Sample information (see ?contrast_each_group_to_the_rest_for_norm_ma_with_limma?for details)

this_dataset_dir<-file.path(datasets_dir, 'haemosphere_datasets','watkins') norm_expression_file<-file.path(this_dataset_dir,"watkins_expression.txt") samples_file<- file.path(this_dataset_dir, "watkins_samples.txt") norm_expression_table.full<-ead.table(norm_expression_file, sep="\t", header=TRUE,quote="",comment.char="",row.names=1,check.names=FALSE) samples_table<- read_tsv(samples_file, col_types = cols()) samples_table$description<- make.names( samples_table$description) # Avoid group or extra_factor names starting with numbers, for microarrays

該數(shù)據(jù)集還包括其他組織，但作為PBMC數(shù)據(jù)的參考，我們只想考慮外周血樣本。

samples_table <- samples_table[samples_table$tissue == "Peripheral Blood",]

這就是大致的流程，最難的部分是數(shù)據(jù)格式。微陣列表達(dá)值應(yīng)使用經(jīng)過標(biāo)準(zhǔn)化，對數(shù)轉(zhuǎn)換并具有相同ID號的數(shù)據(jù)作為query dataset。從haemosphere網(wǎng)站能得到標(biāo)準(zhǔn)化的數(shù)據(jù)?—?但仍需要匹配ID。

該數(shù)據(jù)來自Illumina HumanWG-6 v2 Expression BeadChips，并在探針?biāo)缴辖o出表達(dá)。需要將這些探針轉(zhuǎn)換為gene symbol以匹配PBMC數(shù)據(jù)。

library("tidyverse") library("illuminaHumanv2.db") probes_with_gene_symbol_and_with_data <- intersect(keys(illuminaHumanv2SYMBOL),rownames(norm_expression_table.full))# Get mappings - non NA probe_to_symbol <- select(illuminaHumanv2.db, keys=rownames(norm_expression_table.full), columns=c("SYMBOL"), keytype="PROBEID") probe_to_symbol <- unique(probe_to_symbol[! is.na(probe_to_symbol$SYMBOL),]) # no multimapping probes genes_per_probe <- table(probe_to_symbol$PROBEID) # How many genes a probe is annotated against? multimap_probes <- names(genes_per_probe)[genes_per_probe > 1] probe_to_symbol <- probe_to_symbol[!probe_to_symbol$PROBEID %in% multimap_probes, ]convert_expression_table_ids<- function(expression_table, the_probes_table, old_id_name, new_id_name){the_probes_table <- the_probes_table[,c(old_id_name, new_id_name)]colnames(the_probes_table) <- c("old_id", "new_id")# Before DE, just pick the top expresed probe to represent the gene# Not perfect, but this is a ranking-based analysis.# hybridisation issues aside, would expect higher epressed probes to be more relevant to Single cell data anyway.probe_expression_levels <- rowSums(expression_table)the_probes_table$avgexpr <- probe_expression_levels[as.character(the_probes_table$old_id)]the_genes_table <- the_probes_table %>%group_by(new_id) %>%top_n(1, avgexpr)expression_table <- expression_table[the_genes_table$old_id,]rownames(expression_table) <- the_genes_table$new_idreturn(expression_table) }# Just the most highly expressed probe foreach gene. norm_expression_table.genes <- convert_expression_table_ids(norm_expression_table.full,probe_to_symbol, old_id_name="PROBEID", new_id_name="SYMBOL") # 現(xiàn)在讀取數(shù)據(jù)并用contrast_each_group_to_the_rest_for_norm_ma_with_limma進(jìn)行對比。de_table.Watkins2009PBMCs <- contrast_each_group_to_the_rest_for_norm_ma_with_limma(norm_expression_table = norm_expression_table.genes,sample_sheet_table = samples_table,dataset_name = "Watkins2009PBMCs",extra_factor_name = 'description',sample_name = "sampleId",group_name = 'celltype')</pre>

將10X查詢數(shù)據(jù)與參考數(shù)據(jù)進(jìn)行比較

make_ranking_violin_plot(de_table.test=de_table.10X_pbmc4k_k7,de_table.ref=de_table.Watkins2009PBMCs)

# 提供分組標(biāo)簽 label_table.pbmc4k_k7_vs_Watkins2009PBMCs <- make_ref_similarity_names(de_table.10X_pbmc4k_k7, de_table.Watkins2009PBMCs)

這么熱的單細(xì)胞，中科院的算法開發(fā)博士帶你真正玩轉(zhuǎn)這項平均每個月都有多篇高IF文章的技術(shù)。（點(diǎn)擊藍(lán)色字體了解詳細(xì)）

Reference

Sarah Williams (2019). celaref: Single-cell RNAseq cell cluster labelling by reference. R package version 1.0.1.

https://bioconductor.org/packages/release/bioc/vignettes/celaref/inst/doc/celaref_doco.html

總結(jié)

以上是生活随笔為你收集整理的Celaref | 单细胞测序细胞类型注释工具的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。