在之前的文章裡，我主要講了如下兩個內容：(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_tag7、tag6、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("/home/ruibinxi_pkuhpc/lustre1/huangsiyuan/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返回的標籤不容易識別，比如1、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資料的處理中可以說是相當重要了，如果拆分質量不過關，錯誤地將不同來源的細胞劃分到一個矩陣中，對後續分析結果影響很大。
上述程式碼只呈現了拆分的關鍵步驟，詳細的畫圖程式碼沒有放上來，如果需要可以在微信後臺私信我。
因水平有限，有錯誤的地方，歡迎批評指正！