最近真的是懶得要死，很多東西真的是做了沒有及時整理出來，像這篇文獻(xiàn)已經(jīng)讀過了，方法也去嘗試過了，但是就是沒有整理出來，正所謂好記性不如爛筆頭。因為自己一直在分析單細(xì)胞測序的數(shù)據(jù)，總是畫了聚類圖后，想著怎么鑒定出細(xì)胞類型來，然后剛好讀到這篇文獻(xiàn)，講了用CHETAH，基于已發(fā)表的文獻(xiàn)的鑒定出來的細(xì)胞類型來鑒定，這種感覺就像是以前用已知數(shù)據(jù)來訓(xùn)練未知數(shù)據(jù)一樣，我感覺隨著單細(xì)胞測序數(shù)據(jù)越來越多，這種鑒定方法將成為一種趨勢，以后應(yīng)該會用更好的算法來實現(xiàn)，我想想自己就是個菜鳥哈哈哈，能用能看懂就很開心呢了。
一、文章總體思路
文章的思路：
A：利用已知的數(shù)據(jù)，構(gòu)建一棵分類樹。比如一套數(shù)據(jù)中有四種細(xì)胞類型 RP1 RP2 RP3 RP4,文章首先是計算這四種細(xì)胞類型的平均值，然后算著四種類型細(xì)胞的相關(guān)性，然后進(jìn)行層次聚類，構(gòu)建出分類樹。
B: 拿到一個cell j，比如假定說屬于RP1 類，那么久計算這個細(xì)胞和RP1的相關(guān)性，也要計算cell j與除RP1 以外的 【RP3 RP4】（將這兩個當(dāng)作一個整體）的相關(guān)性，為什么不計算RP2呢，因為我是這么想的，RP2 和RP1屬于都一個根下，所以都加進(jìn)去計算會不太好吧，具體原理不是很清楚。這里計算完之后，因為這個算法是對每個細(xì)胞選定前200個差異基因（RP1 對RP3 RP4的差異基因，文章說效果比較好）來計算相關(guān)性，因此每個細(xì)胞的200個基因是不一樣，所以作者這邊就計算RP1 【RP3 RP4】的理論分布，然后把相關(guān)性轉(zhuǎn)化為落在理論分布的profile score，之后還計算了像p值一樣的confidence score。
C：就是對每個細(xì)胞做上述過程，然后得到cell type.
細(xì)胞類型的鑒定過程取決于你所選的reference data的可靠性咯。
作者還比較這個算法與其他方法的結(jié)果，與之前的SingleR的結(jié)果不相上下，但是SingleR是用bulk RNA數(shù)據(jù)。
二、代碼實現(xiàn)
1. 網(wǎng)址
vignette("CHETAH_introduction")
2. 輸入輸出
1）輸入
If you have your data stored as SingleCellExperiments, continue to the next step. Otherwise, you need the following data before you begin:
- 輸入你要分析的單細(xì)胞數(shù)據(jù)矩陣【 input scRNA-seq count data of the cells to be classified a data.frame or matrix, with cells in the columns and genes in the rows】
- 標(biāo)準(zhǔn)化后的參考的單細(xì)胞數(shù)據(jù)矩陣(!) 【normalized scRNA-seq count data of reference cells in similar format as the input】
- 參考的單細(xì)胞數(shù)據(jù)的細(xì)胞類型【the cell types of the reference cells
a named character vector (names corresponding to the colnames of the reference counts)】
- 可選用來可視化的要分析的單細(xì)胞數(shù)據(jù)的2維的TSNE PCA數(shù)值【(optional) a 2D reduced dimensional representation of your input data for visualization: e.g. t-SNE1, PCA，a two-column matrix/data.frame, with the cells in the rows and the two dimensions in the columns】
2）輸出
CHETAH constructs a classification tree by hierarchically clustering the reference data. The classification is guided by this tree. In each node of the tree, input cells are either assigned to the right, or the left branch. A confidence score is calculated for each of these assignments. When the confidence score for an assignment is lower than the threshold (default = 0.1), the classification for that cell stops in that node.
This results in two types of classifications:
- 1.最終的細(xì)胞類型【 final types: Cells that are classified to one of the leaf nodes of the tree (i.e. a cell type of the reference)】
- 中間細(xì)胞類型，其實就沒辦法細(xì)分的細(xì)胞類型，也就是處于根的時候，左右兩邊的細(xì)胞類型的得分相似?！緄ntermediate types: Cells that had a confidence score lower than the threshold in a certain node are assigned to that intermediate node of the tree. This happens when a cell has approximately the same similarity to the cell types in right and the left branch of that node.】
- 3.中間細(xì)胞類型在圖中的顯示就是Unassigned “Node1”, “Node2”等
【CHETAH generates generic names for the intermediate types: Unassigned for cells that are classified to the very first node, and “Node1”, “Node2”, etc for the additional nodes】
2.代碼實例
（1）安裝加載包
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("CHETAH")
vignette("CHETAH_introduction")
library(CHETAH)
(2) 構(gòu)建reference data
1) 首先你需要下載你要參考的單細(xì)胞的數(shù)據(jù)，我是從Mouse Cell Atlas中下載小鼠腦的數(shù)據(jù)，然后導(dǎo)入R中分析
setwd("數(shù)據(jù)存的路徑/Mouse Cell Atls")
adult_brain_exp <- read.csv("adult_brain_cell_exp.csv",header=T)
adult_cell_type <- read.csv("adult_brain_cell.csv",header=T)
## filter gene by a gene is expressed at least in 50 cells 因為基因數(shù)目比較多，我就把在少于50個細(xì)胞中表達(dá)的基因過濾了
row_sum <- apply(adult_brain_exp[,-1],1,sum)
filter_exp <- adult_brain_exp[row_sum>50,]#一個基因至少在50個細(xì)胞中表達(dá)
genename <- as.character(filter_exp[,1])#genename中有兩個2-Mar 基因名字
final_exp <- as.matrix(filter_exp[,-1])
index <- which(genename=="2-Mar")
genename[index[1]]<- "2-Mar.1"
rownames(final_exp)<-genename
final_exp <- na.omit(final_exp) # 我數(shù)據(jù)中不知道怎么有缺失值，然后之前跑一直不成功，我試了那么久才發(fā)現(xiàn)
2)我把數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化和中心化，我使用seurat來做的，可以自己寫代碼
#Normalization and Scale data and use the Top 3000 variable gene
library(Seurat)
exp_1 <- CreateSeuratObject(counts = final_exp)
exp_1 <- SCTransform(exp_1)
exp_2 <- exp_1@assays$SCT@scale.data
3) 因為參考的數(shù)據(jù)里面可能不同細(xì)胞類型的數(shù)目存在差異，作者認(rèn)為細(xì)胞數(shù)目在10-200之間是比較好，計算量不會太大，所以有一步選擇細(xì)胞的過程。我原始數(shù)據(jù)中有一類細(xì)胞有1000多，如果不選的話，會導(dǎo)致構(gòu)建分類樹的時候，結(jié)果不準(zhǔn)確。
# Cell selction
cell_selection <- unlist(lapply(unique(adult_cell_type$cell.type), function(type) {
type_cells <- which(adult_cell_type$cell.type == type)
if (length(type_cells) > 200) {
type_cells[sample(length(type_cells), 200)]
} else type_cells
}))
4） 構(gòu)建reference data 并用自身數(shù)據(jù)進(jìn)行分類看看結(jié)果
## make reference data
ref_count <- exp_2[,cell_selection]
ref_cellid <- adult_cell_type$cell.id[cell_selection]
ref_celltype <- adult_cell_type$cell.type[cell_selection]
all.equal(as.character(ref_cellid), colnames(ref_count))
adult_brain_ref <- SingleCellExperiment(assays = list(counts = ref_count),
colData = DataFrame(celltypes = ref_celltype))
#Reference QC
CorrelateReference(ref_cells = adult_brain_ref)
ClassifyReference(ref_cells = adult_brain_ref)
(3) 導(dǎo)入自己的數(shù)據(jù)
雖然文章沒有說自己的數(shù)據(jù)需不需要標(biāo)準(zhǔn)化，但是我還是做了，我感覺沒有多大影響
## input data
input_count <- readsCountSM_TSNE@assays$RNA@counts
## Normalization
input_count <- apply(input_count,2,function(column) log2((column/sum(column) * 100000) + 1))
reduced_tsne <- Embeddings(readsCountSM_TSNE, reduction = "tsne")#使用embeddings 函數(shù)來調(diào)用tsne 值
all.equal(rownames(reduced_tsne), colnames(input_count))
input_drug <- SingleCellExperiment(assays = list(counts = input_count),
reducedDims = SimpleList(TSNE = reduced_tsne))
（4）分類，然后把分類結(jié)果導(dǎo)入從Seurat得到的object中， 就可以了
## Classfication
input_drug <- CHETAHclassifier(input = input_drug,
ref_cells = adult_brain_ref)
PlotCHETAH(input_drug)
input_drug <- Classify(input_drug, 0)
PlotCHETAH(input = input_drug, tree = FALSE)
pdf("cell_type_analyzed_by_adult_brain_type_from_Mouse_cell_atlas.pdf",8,8)
DimPlot(object = readsCountSM_TSNE, reduction = 'tsne',label = TRUE,group.by = 'ident',shape.by="orig.ident",repel=TRUE,label.size=5,pt.size=1.5)
dev.off()
levels(readsCountSM_TSNE$orig.ident)[1]="3-HB"
levels(readsCountSM_TSNE$orig.ident)[2]="Control"
Sample_cluster(readsCountSM_TSNE)
3. 總結(jié)
我感覺這個方法還是不錯的，結(jié)果還是挺準(zhǔn)的。但是還是取決于所參考的數(shù)據(jù)。
