單細胞分析實錄(3): Cell Hashing數據拆分

在以前的文章裏，我主要講了以下兩個內容：(1) 認識Cell Hashing；(2): 使用Cell Ranger獲得表達矩陣。相信你們已經知道了cell hashing與普通10X轉錄組的差別，以及使用cellranger獲得表達矩陣。

這一篇講如何使用Seurat的HTODemux函數，CiteFuse的crossSampleDoublets函數兩種方法拆分表達矩陣（混了不一樣來源的細胞），最後還會略微比較一下兩種方法獲得的結果的差別。

HTODemux

這種方法的原理我在第一篇筆記中已經講過，感興趣的小夥伴能夠看以前的文章。主要R代碼以下：

library(Seurat)
library(ggplot2)
library(tidyverse)
args <- commandArgs(TRUE)

加載R包，導入外部參數，args[1]表示樣本名稱，args[2]表示ensembl_ID和基因symbol對應關係的文本文件，前面獲得的表達矩陣行名是ensembl_ID，爲了在後續可視化的時候更省事，建議在這一步更換基因名稱。

df <- read.table(paste(args[1],".mat.count.txt",sep = ""),header = T,row.names = 1) #df的行數包括基因和tag
colnames(df) <- str_replace(colnames(df),"\\.1","")
ensembl_symbol <- read.table(args[2],header = F,row.names = 1,stringsAsFactors = F)
df1 <- df[intersect(rownames(ensembl_symbol),rownames(df)),] #提取出基因表達矩陣
df2 <- df[setdiff(rownames(df),rownames(ensembl_symbol)),] #提取出tag表達矩陣
rownames(df1) <- ensembl_symbol[rownames(df1),] #更換基因表達矩陣的行名

接下來利用df2數據框的信息拆分，df2行爲tag，列爲cellular barcode

cellhash <- CreateSeuratObject(counts = df2,project = "cell_hashing", assay = "HTO")
cellhash <- NormalizeData(cellhash, assay = "HTO", normalization.method = "CLR")
cellhash <- HTODemux(cellhash, assay = "HTO", positive.quantile = 0.85)

最後一步就是拆分，第一篇筆記說過，positive.quantile參數表示在擬合負二項分佈以後使用什麼分位數來判斷UMI是相對大仍是相對小，默認值是0.99，實際處理時，發現這個值可能並不合理，好比最終拆分出來的有效細胞數、不一樣來源細胞數比例和預期差異很大，再好比從圖形上看，明顯不對（下文有圖形說明）。
這一步以後，每個CB都會帶上一個標籤，好比個人數據只有兩個樣原本源，標籤會有這4種：Negative、tag6_tag七、tag六、tag7，前面兩個表示空液滴、(跨樣本的)doublet。

Idents(cellhash) <- "HTO_classification"
FeatureScatter(cellhash, feature1 = paste("hto_",rownames(cellhash)[1],sep=""),
               feature2 = paste("hto_",rownames(cellhash)[2],sep = ""),slot = "counts")

HTOHeatmap(cellhash, assay = "HTO")

上面兩個圖，能夠用來檢驗拆分的質量，第一張每一個點的橫縱座標表示每一個CB兩個tag的UMI，第二張圖的每一列表示每一個CB兩個tag的標準化以後的表達量。
而後根據每一個CB的標籤提取出有效的singlet就能夠了。

small_df1 <- df1[,colnames(cellhash)[cellhash$HTO_classification==rownames(cellhash)[1]]]
write.table(small_df1,paste(args[1],"_",rownames(cellhash)[1],".mat.count.txt",sep = ""),quote = F,row.names = T,col.names = T,sep="\t")
small_df2 <- df1[,colnames(cellhash)[cellhash$HTO_classification==rownames(cellhash)[2]]]
write.table(small_df2,paste(args[1],"_",rownames(cellhash)[2],".mat.count.txt",sep = ""),quote = F,row.names = T,col.names = T,sep="\t")

除了上面兩種Seurat自帶圖形，下面兩種圖形也頗有參考意義，代碼就先不放了，若有須要能夠在公衆號後臺小窗我。

將UMI取對數以後作圖，能夠從另外一個角度看結果，能夠看到右上角被HTODemux認定爲doublet的CB，像是包含了本應該是singlet的CB。我嘗試過positive.quantile用默認值0.99，這種現象會更明顯，因此我以爲在作這一步的時候，能夠畫畫這個圖，選擇一個適中的positive.quantile值。

crossSampleDoublets

CiteFuse包在作這一步的時候，是從取對數以後的UMI矩陣開始的，分別從兩個維度擬合正態分佈，所以最終獲得的結果在散點圖上，比上一種方法更說得過去。示意圖以下：

具體使用的R代碼以下：

library(tidyverse)
library(ggplot2)
library(SingleCellExperiment)
library(CiteFuse)
args <- commandArgs(TRUE)

allexp <- read.table(paste(args[1],".mat.count.txt",sep = ""),header = T,row.names = 1)
colnames(allexp) <- str_replace(colnames(allexp),"\\.1","")
allexp_sce <- preprocessing(exprsMat = as.matrix(allexp)) #生成特定的對象
is.HTO <- grepl("^tag[123678]", rownames(allexp_sce)) #根據本身的tag命名修改正則表達式
allexp_sce <- splitAltExps(allexp_sce, ifelse(is.HTO, "HTO", "gene")) #給每一行加一個標籤，HTO或者gene

allexp_sce=normaliseExprs(allexp_sce, altExp_name = "HTO", exprs_value = "counts",transform = c("log")) #僅針對HTO行，取對數
allexp_sce=crossSampleDoublets(allexp_sce,altExp_name = "HTO",totalExp_threshold = 10)

最後一行就是拆分關鍵步驟，會給每一個CB一個標籤，totalExp_threshold表示只會保留表達數大於10的CB。

ensembl_symbol <- read.table("/ref/10x/Ensembl_symbol_new.txt",header = F,row.names = 1,stri
ngsAsFactors = F)
df1 <- allexp[intersect(rownames(ensembl_symbol),rownames(allexp)),]
df2 <- allexp[setdiff(rownames(allexp),rownames(ensembl_symbol)),]
rownames(df1) <- ensembl_symbol[rownames(df1),]

tmp1=as.data.frame(t(df2))
tmp2=as.data.frame(allexp_sce$doubletClassify_between_label)
colnames(tmp2)="anno"
tmp2$anno=as.character(tmp2$anno)

crossSampleDoublets返回的標籤不容易識別，好比一、2，還須要從新更換名稱，以下

for (i in seq(1,length(rownames(tmp2)),1)) {
  for (j in seq(1,length(colnames(tmp2)),1)) {
    if (tmp2[i,j] == "1") {
      tmp2[i,j] = colnames(tmp1)[1]
    }
    if (tmp2[i,j] == "2") {
      tmp2[i,j] = colnames(tmp1)[2]
    }
    if (tmp2[i,j] == "doublet/multiplet") {
      tmp2[i,j] = "doublet"
    }
  }
}

df_point=cbind(tmp1,tmp2)
colnames(df_point)=c("taga","tagb","anno")

這一步以後就能根據tag標籤畫散點圖，以及提取想要的矩陣了。

---

實際處理中，上面兩種方法我都用了，最後選了兩者交集的CB來提取矩陣（相對保險的作法）。這一步在cell hashing數據的處理中能夠說是至關重要了，若是拆分質量不過關，錯誤地將不一樣來源的細胞劃分到一個矩陣中，對後續分析結果影響很大。
上述代碼只呈現了拆分的關鍵步驟，詳細的畫圖代碼沒有放上來，若是須要能夠在微信後臺私信我。
因水平有限，有錯誤的地方，歡迎批評指正！