單細(xì)胞免疫組庫VDJ——scRepertoire （四）

一、介紹

scRepertoire 旨在從 10x Genomics Cell Ranger 管道中獲取過濾器 contig 輸出，處理該數(shù)據(jù)以根據(jù)兩個 TCR 或 Ig 鏈分配克隆型，并分析克隆型動態(tài)。后者可以分為 1) 僅克隆型分析功能，例如獨(dú)特的克隆型或克隆空間量化和 2) 使用 Seurat、SingleCellExperiment 或 Monocle 3 包與 mRNA 表達(dá)數(shù)據(jù)能夠交互。

1、安裝

> devtools::install_github("ncborcherding/scRepertoire")
> packageVersion('scRepertoire')
[1] ‘1.5.2’

2、加載庫

suppressMessages(library(scRepertoire))
suppressMessages(library(Seurat))

二、數(shù)據(jù)處理

1、加載和處理 Contig Data

使用來自 10x Genomics Cell Ranger的過濾后的_contig_annotations.csv 輸出。此文件位于 VDJ 對齊文件夾的 ./outs/ 目錄中。要生成用于 scRepertoire 的重疊群列表，只需 1) 為每個樣本加載 filters_contig_annotations.csv，然后2) 制作一個列表。

S1 <- read.csv(".../Sample1/outs/filtered_contig_annotations.csv")
S2 <- read.csv(".../Sample2/outs/filtered_contig_annotations.csv")
S3 <- read.csv(".../Sample3/outs/filtered_contig_annotations.csv")
S4 <- read.csv(".../Sample4/outs/filtered_contig_annotations.csv")
contig_list <- list(S1, S2, S3, S4)

2、Combining the Contigs

由于 CellRanger 的輸出是 TCRA 和 TCRB 鏈的量化，下一步是通過細(xì)胞條形碼創(chuàng)建一個包含 TCR 基因（由 VDJC 基因組成）和 CDR3 序列的單個列表對象。這是使用 執(zhí)行的combineTCR()，其中輸入是剝離的 contig_list。還可以通過樣本和 ID 信息重新標(biāo)記條形碼，以防止重復(fù)。

#若有分組信息，則填入samples中，ID可填入樣本名稱，cells：T-AB 表示T 細(xì)胞的α-β；
T-GD 表示T細(xì)胞的γ-δ ，一般此處選擇T-AB。
combined <- combineTCR(contig_list, 
                samples = c("group1", "group1", "group2", "group2"), 
                ID = c("S1", "S2", "S3", "S4"), cells ="T-AB")

處理BCR時，對應(yīng)函數(shù)combineBCR()，但有兩個主要注意事項：1）每個條形碼最多只能有 2 個序列，如果存在更多，則選擇具有最高讀數(shù)的 2 個。2) 克隆型的嚴(yán)格定義 (CTstrict) 基于 V 基因和核苷酸序列的 >85% 標(biāo)準(zhǔn)化 Levenshtein 距離。Levenshtein 距離是在所有恢復(fù)的 BCR 序列中計算的，無論運(yùn)行如何。

3、添加附加變量

如果要添加的變量不僅僅是樣本和 ID 怎么辦？我們可以使用addVariable()函數(shù)添加它們。我們需要的只是您要添加的變量的名稱以及特定的字符或數(shù)值（變量）。例如，在這里我們添加了處理和測序樣品的批次。

example <- addVariable(combined, name = "batch",  variables = c("b1", "b2", "b1", "b2"))
example[[1]][1:5,ncol(example[[1]])] # This is showing the first 5 values of the new column added

同樣，我們可以在使用該subsetContig()函數(shù)后刪除特定的列表元素。為了進(jìn)行子集化，我們需要確定要用于子集化的向量（名稱）以及子集的變量值（變量）。

subset <- subsetContig(combined, name = "sample", variables = c("group1", "group2"))

4、可視化及高級分析

4.1、Quantify Clonotypes

探索克隆型的第一個函數(shù)是quantContig()返回唯一克隆型的總數(shù)或相對數(shù)量。先來個圖我們再說參數(shù)：

quantContig(combined, cloneCall="gene+nt", scale = T)+theme(axis.text.x = element_text(angle = 30,vjust = 0.85,hjust = 0.75))#X坐標(biāo)加點(diǎn)斜體，出圖好看
quantContig(combined, cloneCall="gene+nt", scale = F)+theme(axis.text.x = element_text(angle = 30,vjust = 0.85,hjust = 0.75))

image.png

兩張圖的唯一區(qū)別就是scale參數(shù)。其余參數(shù)的對比不一一展示。
其余參數(shù)含義如下：

  • cloneCall :這里有四個選項
    (1) “gene” ：使用包含 TCR/Ig 的 VDJC 基因
    (2) “nt”：使用 CDR3 區(qū)域的核苷酸序列
    (3) “aa” ：使用 CDR3 區(qū)域的氨基酸序列
    (4) “gene+nt” 使用包含 TCR/Ig + CDR3 區(qū)域的核苷酸序列的 VDJC 基因。這是克隆型的正確定義
    需要注意的重要一點(diǎn)是，克隆型基本上是使用兩個基因座的基因或 nt/aa CDR3 序列的組合來調(diào)用的。在 scRepertoire 的這種實施中，克隆型調(diào)用并未在 CDR3 序列中包含小的變化。因此，基因方法將是最敏感的，而nt或aa的使用則較為敏感，而對克隆型最特異的是gene+nt。此外，克隆型調(diào)用試圖合并兩個基因座，即TCRA和TCRB鏈，并且如果單個細(xì)胞條形碼具有識別的多個序列（即, 2 個 TCRA 鏈在一個細(xì)胞中表達(dá)）。使用 10x 方法，有一個條形碼子集只返回一個免疫受體鏈，未返回的鏈被分配一個NA值。

  • scale：將圖表轉(zhuǎn)換為獨(dú)特克隆型的百分比。

  • group：用于分組的列標(biāo)題。

  • exportTable：返回用于形成圖形的數(shù)據(jù)框

exportTable=T，導(dǎo)出數(shù)據(jù)

> quantContig_output <- quantContig(combined, cloneCall="gene+nt",    scale = T, exportTable = T)
> quantContig_output
  contigs    values total   scaled
1    1552 group1_S1  7620 20.36745
2    2678 group1_S2  8336 32.12572
3    1093 group2_S3  3064 35.67232
4     941 group2_S4  8125 11.58154

4.2、克隆型豐度

通過豐度檢查克隆型的相對分布。abundanceContig()將根據(jù)樣本或運(yùn)行中的實例數(shù)生成具有克隆型總數(shù)的折線圖。

abundanceContig(combined, cloneCall = "gene", scale = F)
image.png

4.3、克隆型的長度

lengtheContig()我們可以通過調(diào)用函數(shù)查看 CDR3 序列的長度分布。重要的是，與其他基本可視化不同，cloneCall只能是“nt”或“aa”。由于如上所述調(diào)用克隆型的方法，長度應(yīng)顯示多峰曲線，這是對未返回的鏈序列和單個條形碼中的多個鏈?zhǔn)褂肗A的產(chǎn)物。

> lengthContig(combined, cloneCall="aa", chain = "both")#代碼來自官方，出現(xiàn)報錯
Error in seq.default(min(Con.df$length), max(Con.df$length), by = 5) :
  'from' must be a finite number
In addition: Warning messages:
1: In min(Con.df$length) : no non-missing arguments to min; returning Inf
2: In max(Con.df$length) : no non-missing arguments to max; returning -Inf
> ?lengthContig#發(fā)現(xiàn)chains出的參數(shù)改變
Usage:

     lengthContig(
       df,
       cloneCall = "aa",
       group = NULL,
       scale = FALSE,
       chains = "combined",
       exportTable = FALSE
     )

Arguments:

      df: The product of combineTCR(), combineBCR(), or
          expression2List()

cloneCall: How to call the clonotype - CDR3 nucleotide (nt), CDR3 amino
          acid (aa).

   group: The group header for which you would like to analyze the
          data.

   scale: Converts the graphs into denisty plots in order to show
          relative distributions.

  chains: Whether to keep clonotypes "combined" or visualize by chain.

exportTable: Returns the data frame used for forming the graph.
#通過修改參數(shù)得到正確圖形。
lengthContig(combined, cloneCall="aa", chain = "combined")#左圖
lengthContig(combined, cloneCall="aa", chain = "single") #右圖
新版本結(jié)果圖，左圖展示的是一起統(tǒng)計，右圖展示的是單鏈統(tǒng)計

4.4 比較克隆型

compareClonotypes()查看樣本之間的克隆型和動態(tài)變化。

compareClonotypes(combined, numbers = 15, samples = c("group1_S1", "group1_S2"),  cloneCall="aa", graph = "alluvial")

Arguments:
df: The product of combineTCR(), combineBCR(), or
          expression2List()
cloneCall: How to call the clonotype - CDR3 gene (gene), CDR3
          nucleotide (nt), CDR3 amino acid (aa), or CDR3
          gene+nucleotide (gene+nt).
samples: 可用于根據(jù)列表元素的名稱隔離特定樣本
clonotypes: 可用于分離特定的克隆型序列，確保調(diào)用與您想要可視化的序列匹配。
numbers: 要繪制的最高克隆型數(shù)，這將根據(jù)單個樣本的頻率計算。這也可以留空。
graph: The type of graph produced, either "alluvial" or "area".
exportTable: Returns the data frame used for forming the graph.

4.5 可視化基因使用

vizGenes()用來查看基因在單鏈中的使用情況，這個函數(shù)在新版本中沒有，‘1.5.2’這個版本中有。

vizGenes(combined[c(1,3,5)], #選擇樣本，
         gene = "V",  chain = "TRB",   y.axis = "J", # 表示TRB V 和 J 使用的差異
         plot = "heatmap", #繪制熱圖，還可以選“bar”，繪制柱狀圖
         scale = TRUE, 
         order = "gene")
image.png

4.6 克隆空間穩(wěn)態(tài)

clonalHomeostasis(combined, cloneCall = "gene",
                  cloneTypes = c(Rare = 1e-04,
                                 Small = 0.001,
                                 Medium = 0.01,
                                 Large = 0.1,
                                 Hyperexpanded = 1))
image.png

4.7 克隆比例

clonalProportion(combined, cloneCall = "nt")
image.png

4.8 Overlap Analysis

使用clonalOverlap()可以可視化樣本之間的相似性。目前可以執(zhí)行三種方法clonalOverlap()：

  • 1）overlap coefficient
  • 2）Morisita指數(shù)
  • 3）Jaccard指數(shù)
    前者正在研究根據(jù)較小樣本中獨(dú)特克隆型長度縮放的克隆型重疊。Morisita 指數(shù)更為復(fù)雜，它是個體在群體中分散程度的生態(tài)測量，包括人口規(guī)模。Jaccard Similarity Index 與重疊系數(shù)非常相似 - Jaccard Index 的分母不是使用較小樣本的長度，而是兩個比較的并集，導(dǎo)致數(shù)量要小得多。參數(shù)method = c("overlap", "morisita", "jaccard", "raw")
clonalOverlap(combined, cloneCall = "gene+nt",   method = "overlap")
自己的數(shù)據(jù)結(jié)果

4.9 多樣性分析

多樣性是使用五個指標(biāo)計算的：

    1. Shannon
    1. inverse Simpson
    1. Chao1
    1. Abundance-based Coverage Estimator (ACE)
    1. inverse Pielou’s measure of species evenness #此版本計算后并沒有這個指標(biāo)的結(jié)果。
      前兩者一般用于估計基線多樣性，Chao/ACE 指數(shù)用于估計樣本的豐富度。
p<-clonalDiversity(combined, cloneCall = "gene+nt",  n.boots = 100)
image.png

4.10 Scatter Compare

scatterClonotype(combined, cloneCall ="gene", 
                 x.axis = "PY_P", 
                 y.axis = "PY_T",
                 dot.size = "total",
                 graph = "proportion")
網(wǎng)站說明圖片，自己運(yùn)行報錯 could not find function "scatterClonotype"

參考：
http://bioconductor.riken.jp/packages/release/bioc/html/scRepertoire.html
https://ncborcherding.github.io/vignettes/vignette.html

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容