一、介绍

scRepertoire 旨在从 10x Genomics Cell Ranger 管道中获取过滤器 contig 输出，处理该数据以根据两个 TCR 或 Ig 链分配克隆型，并分析克隆型动态。后者可以分为 1) 仅克隆型分析功能，例如独特的克隆型或克隆空间量化和 2) 使用 Seurat、SingleCellExperiment 或 Monocle 3 包与 mRNA 表达数据能够交互。

1、安装

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

2、加载库

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

二、数据处理

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 链的量化，下一步是通过细胞条形码创建一个包含 TCR 基因（由 VDJC 基因组成）和 CDR3 序列的单个列表对象。这是使用 执行的combineTCR()，其中输入是剥离的 contig_list。还可以通过样本和 ID 信息重新标记条形码，以防止重复。

#若有分组信息，则填入samples中，ID可填入样本名称，cells：T-AB 表示T 细胞的α-β；
T-GD 表示T细胞的γ-δ ，一般此处选择T-AB。
combined <- combineTCR(contig_list, 
                samples = c("group1", "group1", "group2", "group2"), 
                ID = c("S1", "S2", "S3", "S4"), cells ="T-AB")

处理BCR时，对应函数combineBCR()，但有两个主要注意事项：1）每个条形码最多只能有 2 个序列，如果存在更多，则选择具有最高读数的 2 个。2) 克隆型的严格定义 (CTstrict) 基于 V 基因和核苷酸序列的 >85% 标准化 Levenshtein 距离。Levenshtein 距离是在所有恢复的 BCR 序列中计算的，无论运行如何。

3、添加附加变量

如果要添加的变量不仅仅是样本和 ID 怎么办？我们可以使用addVariable()函数添加它们。我们需要的只是您要添加的变量的名称以及特定的字符或数值（变量）。例如，在这里我们添加了处理和测序样品的批次。

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()函数后删除特定的列表元素。为了进行子集化，我们需要确定要用于子集化的向量（名称）以及子集的变量值（变量）。

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

4、可视化及高级分析

4.1、Quantify Clonotypes

探索克隆型的第一个函数是quantContig()返回唯一克隆型的总数或相对数量。先来个图我们再说参数：

quantContig(combined, cloneCall="gene+nt", scale = T)+theme(axis.text.x = element_text(angle = 30,vjust = 0.85,hjust = 0.75))#X坐标加点斜体，出图好看
quantContig(combined, cloneCall="gene+nt", scale = F)+theme(axis.text.x = element_text(angle = 30,vjust = 0.85,hjust = 0.75))

image.png

• cloneCall :这里有四个选项
  (1) “gene” ：使用包含 TCR/Ig 的 VDJC 基因
  (2) “nt”：使用 CDR3 区域的核苷酸序列
  (3) “aa” ：使用 CDR3 区域的氨基酸序列
  (4) “gene+nt” 使用包含 TCR/Ig + CDR3 区域的核苷酸序列的 VDJC 基因。这是克隆型的正确定义
  需要注意的重要一点是，克隆型基本上是使用两个基因座的基因或 nt/aa CDR3 序列的组合来调用的。在 scRepertoire 的这种实施中，克隆型调用并未在 CDR3 序列中包含小的变化。因此，基因方法将是最敏感的，而nt或aa的使用则较为敏感，而对克隆型最特异的是gene+nt。此外，克隆型调用试图合并两个基因座，即TCRA和TCRB链，并且如果单个细胞条形码具有识别的多个序列（即, 2 个 TCRA 链在一个细胞中表达）。使用 10x 方法，有一个条形码子集只返回一个免疫受体链，未返回的链被分配一个NA值。

• scale：将图表转换为独特克隆型的百分比。

• group：用于分组的列标题。

• exportTable：返回用于形成图形的数据框

exportTable=T，导出数据

> 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()将根据样本或运行中的实例数生成具有克隆型总数的折线图。

abundanceContig(combined, cloneCall = "gene", scale = F)

image.png

4.3、克隆型的长度

lengtheContig()我们可以通过调用函数查看 CDR3 序列的长度分布。重要的是，与其他基本可视化不同，cloneCall只能是“nt”或“aa”。由于如上所述调用克隆型的方法，长度应显示多峰曲线，这是对未返回的链序列和单个条形码中的多个链使用NA的产物。

> lengthContig(combined, cloneCall="aa", chain = "both")#代码来自官方，出现报错
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#发现chains出的参数改变
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.
#通过修改参数得到正确图形。
lengthContig(combined, cloneCall="aa", chain = "combined")#左图
lengthContig(combined, cloneCall="aa", chain = "single") #右图

新版本结果图，左图展示的是一起统计，右图展示的是单链统计

4.4 比较克隆型

compareClonotypes()查看样本之间的克隆型和动态变化。

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: 可用于根据列表元素的名称隔离特定样本
clonotypes: 可用于分离特定的克隆型序列，确保调用与您想要可视化的序列匹配。
numbers: 要绘制的最高克隆型数，这将根据单个样本的频率计算。这也可以留空。
graph: The type of graph produced, either "alluvial" or "area".
exportTable: Returns the data frame used for forming the graph.

4.5 可视化基因使用

vizGenes()用来查看基因在单链中的使用情况，这个函数在新版本中没有，‘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 克隆空间稳态

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()可以可视化样本之间的相似性。目前可以执行三种方法clonalOverlap()：

• 1）overlap coefficient
• 2）Morisita指数
• 3）Jaccard指数
  前者正在研究根据较小样本中独特克隆型长度缩放的克隆型重叠。Morisita 指数更为复杂，它是个体在群体中分散程度的生态测量，包括人口规模。Jaccard Similarity Index 与重叠系数非常相似 - Jaccard Index 的分母不是使用较小样本的长度，而是两个比较的并集，导致数量要小得多。参数method = c("overlap", "morisita", "jaccard", "raw")

clonalOverlap(combined, cloneCall = "gene+nt",   method = "overlap")

自己的数据结果

4.9 多样性分析

多样性是使用五个指标计算的：

• 1. Shannon
• 2. inverse Simpson
• 3. Chao1
• 4. Abundance-based Coverage Estimator (ACE)
• 5. inverse Pielou’s measure of species evenness #此版本计算后并没有这个指标的结果。
     前两者一般用于估计基线多样性，Chao/ACE 指数用于估计样本的丰富度。

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")

网站说明图片，自己运行报错 could not find function "scatterClonotype"

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