1、编程基础 整个分析过程中需要掌握 R语言 和 Linux 的，因此，笔者已整理过一些大家常用的R语言与Linux的命令供大家学习。R语言基础学习手册 2、硬件基础 教程制作过程中发现，你至少需要准备一个 4 线程、 16 GB运行内存的电脑来运行以下测试数据，不然可能有电脑死机的风险：不推荐用自己的电脑做生信。实战的数据分析可能需要 上百GB 的内存，如果你的设备不能满足，可以考虑下列服务器的租赁。
最便宜的版本即可帮助你们完成90%以上的单细胞分析需求(土豪随意)：足够支持你完成硕博生涯的生信环境 如何拥有一台服务器，且不用自己搭建维护 为实验室准备一份生物信息学不动产 配置一个心仪的工作站(硬件+环境配置) 3、生信基础 scTCR (Single-cell T Cell Receptor sequencing)和 scBCR (Single-cell B Cell Receptor sequencing)基本都会和 scRNA-seq 的技术同屏出现，所以建议大家在学习这部分教程前先行学习我们此前的 scRNA-seq 系列教程：保姆级教程 《手把手教你做单细胞转录组测序数据分析》
二、免疫组库测序 在制作每个教程时，我们都建议大家了解分析背后的背景知识，这里我们先通过一篇发表在《 Nature Methods 》上的综述文章了解单细胞免疫组库测序的原理与应用: Single-cell immune repertoire analysis 1、TCR和BCR 适应性免疫系统依靠免疫组库的巨大多样性来识别和应对广泛的病原体和外来物质。
这是由大量表面结合的 T细胞抗原受体(TCR) 和 B细胞抗原受体(BCR) 介导的。后者也可在B细胞分化为浆细胞或浆母细胞后作为可溶性抗体分泌。这些受体广泛的抗原特异性也使T细胞和B细胞能够区分自我和非自我抗原，使免疫系统能够对威胁做出适当的反应，同时使宿主不受伤害。2、结构组成 每个TCR或BCR是由两条不同的链组成的二聚体。
TCR由αβ T细胞中的 α链(TRA) 与 β链(TRB) 配对组成，或γδ T细胞中的 γ链(TRG) 与 δ链(TRD) 配对组成。BCRs由 重链(IgH) 和 轻链 组成；轻链来自 κ(IgK) 或 λ(IgL) 基因座。这些链是由 RAG1 和 RAG2 蛋白协调的 可变(V)、多样性(D)和连接(J) 基因片段的基因重组的复杂过程的产物，分别发生在胸腺和骨髓中发育的T细胞和B细胞中。
在这些片段的交界处还可以引入随机的 无模板法(N) 和或 回文(P) 核苷酸插入，进一步增加了复杂性。这些重组事件发生在被称为 互补决定区(CDR)3 的交界处，而 CDR1 和CDR2 则完全在V基因区域内。由于 CDRs 是与同源抗原结合的区域，它们，特别是 CDR3 区域，已成为大多数下游分析的焦点。
此外，B细胞活化后，在活化诱导的胞苷脱氨酶(cytidine deaminase，AID)的介导下，BCRs可在整个受体中发生(随机)体细胞超突变(somatic hypermutations，SHMs)。AID也是BCR的类别转换重组(CSR)所必需的，这是一个生物学过程，它取代了编码同型类的BCR恒定基因，对B细胞成熟和整体体液免疫产生广泛的影响。
3、单细胞免疫组库测序 有多种”批量”高通量技术用于免疫组库测序，涉及从给定组织中分析混合的 BCR 或 TCR。Bulk’s免疫组库测序技术在很大程度上仅限于分析单链（例如,只有 TRB 或 IgH ），并不能捕获适应性免疫受体的二聚体性质。文库构建和测序策略也阻碍了使用”大容量”方法恢复真正的配对链测序。通过单细胞技术，使成对链的适应性免疫受体在规模上的分析成为可能。各种高通量测序 TCRs 和 BCRs 的技术。
针对 BCR 和 TCR 的VDJ序列的异质性，10X genomics公司推出了在cDNA的基础上，进一步进行扩增子测序，通过分析 BCR 和 TCR 的不同序列，分析免疫环境中的克隆型多样性，解析有机体抵制病原微生物入侵的免疫机制。10X genomics的 单细胞BCR/TCR 测序的原理，与单细胞转录组类似，通过微流控系统将带有条形码和引物的凝胶珠与单个细胞包裹在油滴中；
接下来在每个油滴内，凝胶珠溶解，细胞裂解释放 mRNA，通过逆转录产生带有 10X barcode 和 UMI 信息的 cDNA。破油（Breaking Emulsions）之后，cDNA 一分为二，后续同时进行基因表达和免疫组库的文库构建。其中 TCR 或者 BCR 的 V(D)J 序列通过设计在 C 区的巢式引物进行 PCR 富集；而 mRNA 的信息，与 10X Genomics 3’mRNA 文库不同，保留的是 5’端 的信息。
测序后即可一次性获得大量单细胞的基因表达和免疫组库数据。结合单细胞转录组测序数据与 BCR/TCR 多样性数据，解析免疫细胞中更高维度的异质性。
4、分析流程 处理 scTCR/BCR-seq 数据的典型工作流程如下：• 数据准备：获取10x单细胞V(D)J测序数据（ filtered_contig_annotations.csv / filtered_contig_annotations.json ）,同时获取对应的scRNA-seq表达矩阵（用于后续整合）,准备元数据（metadata，如样本分组信息） • 重新注释：标准化链类型名称（IgH/IgK/IgL）,接下来注释细胞克隆类型（如完整重链、轻链是否都存在），分别提取 V、J、C、CDR3 区域信息，准确识别TCR/BCR配对信息和克隆特征，确保分析可靠性。
• 克隆聚类与筛选：去除非配对链、低质量链，基于CDR3序列进行相似性聚类（相同CDR3表示同源克隆），统计克隆频数、构建免疫克隆谱图，揭示免疫细胞克隆扩增特征。• 单细胞免疫受体与表达谱整合分析：将克隆信息整合进Seurat对象或Scanpy对象中，联合UMAP图上可视化免疫克隆空间分布，分析免疫细胞类型与克隆扩增的关系。
• 高级整合分析：富集分析：是否特定克隆富集在某亚群 轨迹分析：追踪克隆扩增细胞的分化路径（使用 monocle3 / scvelo 拟时序分析克隆动态） 疾病相关性：对比健康 vs 疾病样本克隆谱（分析CDR3长度、多样性指数（Shannon, Simpson））
三、cellranger VDJ CellRanger vdj 可用于分析 Chromium Single Cell 5' V (D) J 库产生的测序数据。它从 V (D) J 库中获取 FASTQ 文件，并执行序列组装和配对克隆型调用。它使用 Chromium 细胞条形码和 UMI 来组装每个细胞的 V (D) J 转录本。
克隆型和 CDR3 序列以 .vloupe 文件的形式输出，可以加载到 Loupe V (D) J Browser 中。其核心功能为 cellranger vdj，主要输出 all_contig_annotations.json 、 filtered_contig_annotations.csv 、 clonotypes.csv 、 consensus.fasta 、 filtered_contig.fasta。
1、scBCR数据的上游定量 1.1 下载 10x genomics 的软件下载是要在线申请的，大家可以访问网站自行申请链接：https://support.10xgenomics.com/single-cell-atac/software/downloads/latest wget -O cellranger-9.0.1.tar.gz "https://cf.10xgenomics.com/releases/cell-exp/cellranger-9.0.1.tar.gz?Expires=1754582697&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA&Signature=***~M6dv50cJKytHQJBHS6SS0UdP7G3Jw-k7m9~m4zyMMvTXxynww3~n-G~de5CNnApnIqsloGjj--YS~~WAGButJs~8nKl~OUtnVOdHaar2WYY15yS2H93oSyDqMYfO-zSXtXiVXtyr3bZM4oZAlBWiRfu69pFRw82myjDwt0QBCSLyth1PoQABLHFI3A23nNlbkIuxSfBYAeQnr~mJdRhHk5eyrDSf9Y75ACtSwmv7vPWexatnXnOOWez36czc75usOa2nXjKeZBdbmV6QjGj0nxIyEEyzHccQ__" 这里的软件下载是动态加密地址，在资料目录中已经替大家下载好了名为 cellranger-9.0.1.tar.gz 的文件（这里推荐版本大于6.1）。
需要注意的是，cellranger 只能在 Linux 环境中运行(十小时学会Linux)，且需要至少8核心、64运行内存、1TB磁盘空间、64位 CentOS/RedHat 7.0 或 Ubuntu 14.04 及以上的环境才能够运行，没有计算环境的同学可参考：足够支持你完成硕博生涯的生信环境 1.2 参考基因组准备 除了软件外，大家还需要下载参考基因组数据，最常用的参考数据为人的 GRCh38 与小鼠的 mm10 参考基因组 wget https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-GRCh38-2020-A.tar.gz tar -xf refdata-gex-GRCh38-2020-A.tar.g wget https://cf.10xgenomics.com/supp/cell-atac/refdata-cellranger-arc-mm10-2020-A-2.0.0.tar.gz tar -xf refdata-cellranger-arc-mm10-2020-A-2.0.0.tar.gz wget https://cf.10xgenomics.com/supp/cell-vdj/refdata-cellranger-vdj-GRCh38-alts-ensembl-5.0.0.tar.gz tar -xf refdata-cellranger-vdj-GRCh38-alts-ensembl-5.0.0.tar.gz 1.3 测试数据下载 这里我们下载了10X 官网提供的人类 B 细胞 多组学 5’测序数据（包括表达和VDJ数据）以及VDJ数据库用作展示 mkdir dataset-multi-practice cd dataset-multi-practice wget https://cf.10xgenomics.com/samples/cell-vdj/6.0.0/sc5p_v2_hs_B_1k_multi_5gex_b_Multiplex/sc5p_v2_hs_B_1k_multi_5gex_b_Multiplex_fastqs.tar -xf sc5p_v2_hs_B_1k_multi_5gex_b_Multiplex_fastqs.tar 如果是非模式动物需要自定义reference_vdj 如果官方gtf文件中没有提供 feature_type 包含 V / D / J / C 区段，则需要使用 IMGT/GENE-DB 数据库自己构建对应的 Gtf 文件 cellranger mkvdjref --genome =< GENOME_NAME > \ --fasta =< path/to/sequences.fasta > \ --genes =< path/to/annotations.gtf > 它会生成一个可用于 --vdj-reference 参数的参考目录，目录结构类似于：refdata-cellranger-vdj-xxx-alts-ensembl-5.0.0/ ├── fasta/ │ └── regions.fa ├── reference.json └── reference_info.txt 还需要创建一个 multi_config.csv 文件去定义文件配置 nano multi_config.csv [gene-expression] # 表达数据相关配置 reference,/path/to/project expect-cells,1000 create-bam,true [vdj] # V(D)J 数据相关配置 reference,/path/to/project [libraries] # 所有原始 fastq 文件的路径和类型 fastq_id,fastqs,lanes,feature_types,subsample_rate sc5p_v2_hs_B_1k_5gex,/path/to/project | 2,gene expression, sc5p_v2_hs_B_1k_b,/path/to/project | 2,vdj, 1.4 软件运行 准备好必须的文件，就可以开始运行 mkdir runs/ cd runs/ cellranger multi --id = HumanB_Cell_multi --csv = ../multi_config.csv 生成这样日志就说明成功运行 Martian Runtime - v4.0.13 Serving UI at 相关链接见专题页 ? auth=*** Running preflight checks ( please wait ) ... 2025-08-07 06:40:20 [ runtime ] ( ready ) ID.HumanB_Cell_multi.SC_MULTI_CS.PARSE_MULTI_CONFIG 2025-08-07 06:40:20 [ runtime ] ( run:local ) ID.HumanB_Cell_multi.SC_MULTI_CS.PARSE_MULTI_CONFIG.fork0.chnk0.main 2025-08-07 06:40:33 [ runtime ] ( chunks_complete ) ID.HumanB_Cell_multi. …… 输出显示 “Pipestance completed successfully!” 时，任务就完成了 web_summary: /path/to/project metrics_summary: /path/to/project } Waiting 6 seconds for UI to do final refresh. Pipestance completed successfully! 检查一下生成的文件目录结构 ── runs └── HumanB_Cell_multi ├── _cmdline ├── extras ├── _filelist ├── _finalstate ├── HumanB_Cell_multi.mri.tgz ├── _invocation ├── _jobmode ├── _log ├── _mrosource ├── outs ├── _perf ├── _perf._truncated_ ├── SC_MULTI_CS ├── _sitecheck ├── _tags ├── _timestamp ├── _uuid ├── _vdrkill └── _versions 这样我们就成功生成了定量的结果文件，接下来可以使用 scRepertoire 对数据进行下游分析 2. scTCR数据的上游定量 除了 scBCR-seq 数据，cellranger同样支持 scTCR-seq 数据的定量分析，这里我们下载一个scTCR-seq测试数据进行类似的分析。
2.1 TCR测试数据下载 cd dataset-multi-practice wget https://cf.10xgenomics.com/samples/cell-vdj/5.0.0/sc5p_v2_hs_T_1k_multi_5gex_t/sc5p_v2_hs_T_1k_multi_5gex_t_fastqs.tar -xf sc5p_v2_hs_T_1k_multi_5gex_t_fastqs.tar 2.2 准备配置文件 类似的接下来准备 multi_config_T.csv nano multi_config_T.csv [gene-expression] # 表达数据相关配置 reference,/path/to/project expect-cells,1000 create-bam,true [vdj] # V(D)J 数据相关配置 reference,/path/to/project [libraries] # 所有原始 fastq 文件的路径和类型 fastq_id,fastqs,lanes,feature_types,subsample_rate sc5p_v2_hs_T_1k_5gex,/path/to/project | 2,gene expression, sc5p_v2_hs_T_1k_T,/path/to/project | 2,vdj, 2.3 软件运行 准备好必备的文件，就可以开始运行 mkdir -p runs/ cd runs/ cellranger multi --id = HumanT_Cell_multi --csv = ../multi_config_T.csv 检查一下生成的文件目录结构 tree -L 3 ── runs ├── multi_config.csv ├── multi_config_T.csv ├── runs │ ├── HumanB_Cell_multi │ │ ├── _cmdline │ │ ├── extras │ │ ├── _filelist │ │ ├── _finalstate │ │ ├── HumanB_Cell_multi.mri.tgz │ │ ├── _invocation │ │ ├── _jobmode │ │ ├── _log │ │ ├── _mrosource │ │ ├── outs │ │ ├── _perf │ │ ├── _perf._truncated_ │ │ ├── SC_MULTI_CS │ │ ├── _sitecheck │ │ ├── _tags │ │ ├── _timestamp │ │ ├── _uuid │ │ ├── _vdrkill │ │ └── _versions │ ├── HumanT_Cell_multi │ │ ├── _cmdline │ │ ├── extras │ │ ├── _filelist │ │ ├── _finalstate │ │ ├── HumanT_Cell_multi.mri.tgz │ │ ├── _invocation │ │ ├── _jobmode │ │ ├── _log │ │ ├── _mrosource │ │ ├── outs │ │ ├── _perf │ │ ├── _perf._truncated_ │ │ ├── SC_MULTI_CS │ │ ├── _sitecheck │ │ ├── _tags │ │ ├── _timestamp │ │ ├── _uuid │ │ ├── _vdrkill │ │ └── _versions │ └── output │ ├── clonal_quant.pdf │ ├── clusters_umap.pdf │ ├── dotplot_clusters.pdf │ └── scRep_example_full.rds ├── sc5p_v2_hs_B_1k_multi_5gex_b_fastqs │ ├── sc5p_v2_hs_B_1k_5gex_fastqs │ │ ├── sc5p_v2_hs_B_1k_5gex_S1_L001_I1_001.fastq.gz │ │ ├── sc5p_v2_hs_B_1k_5gex_S1_L001_I2_001.fastq.gz │ │ ├── sc5p_v2_hs_B_1k_5gex_S1_L001_R1_001.fastq.gz │ │ ├── sc5p_v2_hs_B_1k_5gex_S1_L001_R2_001.fastq.gz │ │ ├── sc5p_v2_hs_B_1k_5gex_S1_L002_I1_001.fastq.gz │ │ ├── sc5p_v2_hs_B_1k_5gex_S1_L002_I2_001.fastq.gz │ │ ├── sc5p_v2_hs_B_1k_5gex_S1_L002_R1_001.fastq.gz │ │ └── sc5p_v2_hs_B_1k_5gex_S1_L002_R2_001.fastq.gz │ └── sc5p_v2_hs_B_1k_b_fastqs │ ├── sc5p_v2_hs_B_1k_b_S1_L001_I1_001.fastq.gz │ ├── sc5p_v2_hs_B_1k_b_S1_L001_I2_001.fastq.gz │ ├── sc5p_v2_hs_B_1k_b_S1_L001_R1_001.fastq.gz │ ├── sc5p_v2_hs_B_1k_b_S1_L001_R2_001.fastq.gz │ ├── sc5p_v2_hs_B_1k_b_S1_L002_I1_001.fastq.gz │ ├── sc5p_v2_hs_B_1k_b_S1_L002_I2_001.fastq.gz │ ├── sc5p_v2_hs_B_1k_b_S1_L002_R1_001.fastq.gz │ └── sc5p_v2_hs_B_1k_b_S1_L002_R2_001.fastq.gz ├── sc5p_v2_hs_B_1k_multi_5gex_b_Multiplex_fastqs.tar ├── sc5p_v2_hs_T_1k_multi_5gex_t │ ├── sc5p_v2_hs_T_1k_5gex_fastqs │ │ ├── sc5p_v2_hs_T_1k_5gex_S1_L001_I1_001.fastq.gz │ │ ├── sc5p_v2_hs_T_1k_5gex_S1_L001_I2_001.fastq.gz │ │ ├── sc5p_v2_hs_T_1k_5gex_S1_L001_R1_001.fastq.gz │ │ ├── sc5p_v2_hs_T_1k_5gex_S1_L001_R2_001.fastq.gz │ │ ├── sc5p_v2_hs_T_1k_5gex_S1_L002_I1_001.fastq.gz │ │ ├── sc5p_v2_hs_T_1k_5gex_S1_L002_I2_001.fastq.gz │ │ ├── sc5p_v2_hs_T_1k_5gex_S1_L002_R1_001.fastq.gz │ │ └── sc5p_v2_hs_T_1k_5gex_S1_L002_R2_001.fastq.gz │ └── sc5p_v2_hs_T_1k_t_fastqs │ ├── sc5p_v2_hs_T_1k_t_S1_L001_I1_001.fastq.gz │ ├── sc5p_v2_hs_T_1k_t_S1_L001_I2_001.fastq.gz │ ├── sc5p_v2_hs_T_1k_t_S1_L001_R1_001.fastq.gz │ ├── sc5p_v2_hs_T_1k_t_S1_L001_R2_001.fastq.gz │ ├── sc5p_v2_hs_T_1k_t_S1_L002_I1_001.fastq.gz │ ├── sc5p_v2_hs_T_1k_t_S1_L002_I2_001.fastq.gz │ ├── sc5p_v2_hs_T_1k_t_S1_L002_R1_001.fastq.gz │ └── sc5p_v2_hs_T_1k_t_S1_L002_R2_001.fastq.gz └── sc5p_v2_hs_T_1k_multi_5gex_t_fastqs.tar 16 directories, 72 files
四、scRepertoire 框架使用户能够轻松地将 scRNA 和免疫分析结合起来，支持使用 10x、AIRR、BD、MiXCR、TRUST4 和 WAT3R 单细胞克隆格式，并与流行的基于 R 的单细胞数据 Seurat 进行交互。
这里我们同样和大家一起阅读一下随 scRepertoire 发布的文章，了解一下这个软件可以做的分析内容，感兴趣的同学可以阅读一下原文：Yang Q, Safina KR, Nguyen KDQ, Tuong ZK, Borcherding N (2025) scRepertoire 2: Enhanced and efficient toolkit for single-cell immune profiling. PLoS Comput Biol 21(6): e1012760. 1、 红斑性头痛/正常皮肤 演示数据 文章使用 scRNA-seq 与单细胞 TCR&BCR 测序（scAIRR）对6位患者的红斑性偏头痛（SKL）和邻近正常皮肤（SKN）进行联合分析。
经质控后，识别出13种免疫细胞类型 （Fig2A），其中浆细胞显示出最显著的克隆扩增现象。在样本192567SKL中，一个浆细胞克隆占据超过50%的比例 （Fig2B），CD4_TEM、CD8_CTL与MAIT细胞也有明显扩增。通过比较患者内不同组织间的TCR克隆共享情况，发现SKN与SKL存在显著重叠，提示组织驻留T细胞可能同时活跃于两种环境 （Fig2C）；
而在样本192566中，外周血（BL）与SKL共享TCR克隆，但与SKN互斥，暗示不同组织间存在特异的T细胞迁移或保留机制。BCR克隆未观察到跨组织共享，可能与采样量低或免疫区隔机制有关。使用Shannon与逆Simpson指数评估克隆多样性 （Fig2D），发现各组织间存在差异。
CDR3氨基酸位点分析表明样本192567中第6与11位以丝氨酸（S）为主 （Fig2E），其对应的Atchley因子（如AF3分子大小与AF5电荷）可能影响抗原结合能力。此外，还评估了氨基酸kmer使用频率、VJ配对与单基因使用等特征 （Fig2F），并利用Levenshtein距离聚类分析CDR3序列，识别潜在的抗原选择驱动机制 （Fig2G）。scRepertoire通过这些多维度分析，系统刻画了复杂组织中免疫克隆动态。
2、数据实操 我们使用前面的10X 参考数据集来进行代码实操 # 安装所需包 "BorchLab/scRepertoire")) # 或者通过BiocManager安装 # 加载所需包 library (scRepertoire) # Loading required package: ggplot2 library (Seurat) # Loading required package: SeuratObject # Loading required package: sp # 'SeuratObject' was built with package 'Matrix' 1.7.1 but the current # version is 1.7.2; it is recomended that you reinstall 'SeuratObject' as # the ABI for 'Matrix' may have changed # # Attaching package: 'SeuratObject' # The following objects are masked from 'package:base': # # intersect, t 2.1 基本克隆可视化 # 因为测试文件中样本数仅有一个，接下来可视化我们使用R包自带的数据文件 data ( "contig_list" ) data ( "scRep_example" ) # 将 Contigs 组合成克隆信息 # CombinedTCR是将原始 TCR 测序数据组织成结构化格式用于 scRepertoire 分析的关键第一步。
# 它允许灵活处理单链和配对链、条形码消歧和初始过滤，生成一个数据帧列表，其中每行代表一个单细胞及其相关的 TCR 克隆类型。
combined.TCR <- combineTCR (contig_list, samples = c ( "P17B" , "P17L" , "P18B" , "P18L" , "P19B" , "P19L" , "P20B" , "P20L" ), removeNA = FALSE , removeMulti = FALSE , filterMulti = FALSE ) # 另外如果想要加上其他信息，可以使用addVariable函数 combined.TCR <- addVariable (combined.TCR, variable.name = "Type" , variables = rep ( c ( "B" , "L" ), 4 )) head (combined.TCR[[ 1 ]]) # barcode sample TCR1 cdr3_aa1 # 1 P17B_AAACCTGAGTACGACG-1 P17B TRAV25.TRAJ20.TRAC CGCSNDYKLSF # 3 P17B_AAACCTGCAACACGCC-1 P17B TRAV38-2/DV8.TRAJ52.TRAC CAYRSAQAGGTSYGKLTF # 5 P17B_AAACCTGCAGGCGATA-1 P17B TRAV12-1.TRAJ9.TRAC CVVSDNTGGFKTIF # 7 P17B_AAACCTGCATGAGCGA-1 P17B TRAV12-1.TRAJ9.TRAC CVVSDNTGGFKTIF # 9 P17B_AAACGGGAGAGCCCAA-1 P17B TRAV20.TRAJ8.TRAC CAVRGEGFQKLVF # 10 P17B_AAACGGGAGCGTTTAC-1 P17B TRAV12-1.TRAJ9.TRAC CVVSDNTGGFKTIF # cdr3_nt1 # 1 TGTGGGTGTTCTAACGACTACAAGCTCAGCTTT # 3 TGTGCTTATAGGAGCGCGCAGGCTGGTGGTACTAGCTATGGAAAGCTGACATTT # 5 TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT # 7 TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT # 9 TGTGCTGTGCGAGGAGAAGGCTTTCAGAAACTTGTATTT # 10 TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT # TCR2 cdr3_aa2 # 1 TRBV5-1.None.TRBJ2-7.TRBC2 CASSLTDRTYEQYF # 3 TRBV10-3.None.TRBJ2-2.TRBC2 CAISEQGKGELFF # 5 TRBV9.None.TRBJ2-2.TRBC2 CASSVRRERANTGELFF # 7 TRBV9.None.TRBJ2-2.TRBC2 CASSVRRERANTGELFF # 9 <NA> <NA> # 10 TRBV9.None.TRBJ2-2.TRBC2 CASSVRRERANTGELFF # cdr3_nt2 # 1 TGCGCCAGCAGCTTGACCGACAGGACCTACGAGCAGTACTTC # 3 TGTGCCATCAGTGAACAGGGGAAAGGGGAGCTGTTTTTT # 5 TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT # 7 TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT # 9 <NA> # 10 TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT # CTgene # 1 TRAV25.TRAJ20.TRAC_TRBV5-1.None.TRBJ2-7.TRBC2 # 3 TRAV38-2/DV8.TRAJ52.TRAC_TRBV10-3.None.TRBJ2-2.TRBC2 # 5 TRAV12-1.TRAJ9.TRAC_TRBV9.None.TRBJ2-2.TRBC2 # 7 TRAV12-1.TRAJ9.TRAC_TRBV9.None.TRBJ2-2.TRBC2 # 9 TRAV20.TRAJ8.TRAC_NA # 10 TRAV12-1.TRAJ9.TRAC_TRBV9.None.TRBJ2-2.TRBC2 # CTnt # 1 TGTGGGTGTTCTAACGACTACAAGCTCAGCTTT_TGCGCCAGCAGCTTGACCGACAGGACCTACGAGCAGTACTTC # 3 TGTGCTTATAGGAGCGCGCAGGCTGGTGGTACTAGCTATGGAAAGCTGACATTT_TGTGCCATCAGTGAACAGGGGAAAGGGGAGCTGTTTTTT # 5 TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT_TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT # 7 TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT_TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT # 9 TGTGCTGTGCGAGGAGAAGGCTTTCAGAAACTTGTATTT_NA # 10 TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT_TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT # CTaa # 1 CGCSNDYKLSF_CASSLTDRTYEQYF # 3 CAYRSAQAGGTSYGKLTF_CAISEQGKGELFF # 5 CVVSDNTGGFKTIF_CASSVRRERANTGELFF # 7 CVVSDNTGGFKTIF_CASSVRRERANTGELFF # 9 CAVRGEGFQKLVF_NA # 10 CVVSDNTGGFKTIF_CASSVRRERANTGELFF # CTstrict # 1 TRAV25.TRAJ20.TRAC;TGTGGGTGTTCTAACGACTACAAGCTCAGCTTT_TRBV5-1.None.TRBJ2-7.TRBC2;TGCGCCAGCAGCTTGACCGACAGGACCTACGAGCAGTACTTC # 3 TRAV38-2/DV8.TRAJ52.TRAC;TGTGCTTATAGGAGCGCGCAGGCTGGTGGTACTAGCTATGGAAAGCTGACATTT_TRBV10-3.None.TRBJ2-2.TRBC2;TGTGCCATCAGTGAACAGGGGAAAGGGGAGCTGTTTTTT # 5 TRAV12-1.TRAJ9.TRAC;TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT_TRBV9.None.TRBJ2-2.TRBC2;TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT # 7 TRAV12-1.TRAJ9.TRAC;TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT_TRBV9.None.TRBJ2-2.TRBC2;TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT # 9 TRAV20.TRAJ8.TRAC;TGTGCTGTGCGAGGAGAAGGCTTTCAGAAACTTGTATTT_NA;NA # 10 TRAV12-1.TRAJ9.TRAC;TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT_TRBV9.None.TRBJ2-2.TRBC2;TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT # Type # 1 B # 3 B # 5 B # 7 B # 9 B # 10 B 构建好VDJ对象后，我们可以进行基础的分析可视化 # 计算和量化细胞克隆（clone）多样性和丰度，帮助你了解免疫细胞中TCR/BCR克隆的分布情况 clonalQuant (combined.TCR, cloneCall= "strict" , chain = "both" , group.by = "sample" , scale = TRUE ) # 计算并展示免疫细胞克隆（TCR/BCR）在样本中绝对或相对的丰度分布情况 clonalAbundance (combined.TCR, cloneCall = "gene" , scale = FALSE ) clonalAbundance (combined.TCR, cloneCall = "gene" , scale = TRUE ) # 可视化 CDR3 序列的长度分布 clonalLength (combined.TCR, cloneCall= "aa" , acid sequence）定义克隆，chain = "both" ) # 以密度图的形式进行缩放 clonalLength (combined.TCR, cloneCall= "aa" , acid sequence）定义克隆，chain = "TRA" , scale = TRUE ) # 比较不同样本或条件之间的免疫克隆（TCR/BCR）多样性和重叠情况 clonalCompare (combined.TCR, top.clones = 10 , samples = c ( "P17B" , "P17L" ), cloneCall= "aa" , graph = "alluvial" ) # 我们还可以选择突出显示特定的克隆，此外，我们可以简化图形，将克隆标记为 “Clone: 1”、“Clone: 2” 等 clonalCompare (combined.TCR, top.clones = 10 , highlight.clones = c ( "CVVSDNTGGFKTIF_CASSVRRERANTGELFF" , "NA_CASSVRRERANTGELFF" ), relabel.clones = TRUE , samples = c ( "P17B" , "P17L" ), cloneCall= "aa" , graph = "alluvial" ) # 如果只想看特定的克隆标记 clonalCompare (combined.TCR, clones = c ( "CVVSDNTGGFKTIF_CASSVRRERANTGELFF" , "NA_CASSVRRERANTGELFF" ), relabel.clones = TRUE , samples = c ( "P17B" , "P17L" ), cloneCall= "aa" , graph = "alluvial" ) # 如果想要比较两个样本（或条件）之间克隆的丰度关系 clonalScatter (combined.TCR, cloneCall = "gene" , # 用V和J基因组合定义克隆 x.axis = "P18B" , # X轴样本名（样本1） y.axis = "P18L" , # Y轴样本名（样本2） dot.size = "total" , # 点的大小表示该克隆在两个样本中总的细胞数（丰度总和） graph = "proportion" ) # 绘制的是克隆丰度比例的散点图 2.2 可视化克隆动力学 # 接下来我们可以查看不同克隆在整个免疫受体库中所占比例的分布情况 # 默认的分类是：# Rare: 0 to 0.0001 # Small: 0.0001 to 0.001 # Medium: 0.001 to 0.01 # Large: 0.01 to 0.1 # Hyperexpanded: 0.1 to 1 clonalHomeostasis (combined.TCR, cloneCall = "gene" ) # 可以看到免疫库是被少量超级克隆主导，还是很多稀有克隆构成 # 也可以通过改变截点来改变分布比例 clonalHomeostasis (combined.TCR, cloneCall = "gene" , cloneSize = c ( Rare = 0.001 , Small = 0.01 , Medium = 0.1 , Large = 0.3 , Hyperexpanded = 1 )) # 当然也可以观察除了sample外其他分类中不同克隆类型的分布比例 clonalHomeostasis (combined.TCR, group.by = "Type" , cloneCall = "gene" ) # 看所有样本中不同排名区间的克隆类型在不同样本中所占比例 clonalProportion (combined.TCR, cloneCall = "gene" ) # 看所有样本中不同排名区间的克隆类型在不同样本中所占比例 clonalProportion (combined.TCR, cloneCall = "nt" , clonalSplit = c ( 1 , 5 , 10 , 100 , 1000 , 10000 )) 2.3 整体序列分布概要 # CDR3 氨基酸位置组成百分比的热图 # 在 CDR3 序列的第 1、2、3、…、20 位上，分别有哪些氨基酸出现，以及它们出现的比例是多少 percentAA (combined.TCR, chain = "TRB" , aa.length = 20 ) # 我们还可以量化 CDR3 序列中氨基酸残基的熵 / 多样性水平。
# positionalEntropy 将百分比 AA 的残基量化与多样性计算相结合。
没 # 有方差的位置将报告一个值为 0, 用于比较同一位置不同样本的氨基酸多样性 positionalEntropy (combined.TCR, chain = "TRB" , method = "norm.entropy" , aa.length = 20 ) # method: # shannon - Shannon Index # inv.simpson - Inverse Simpson Index # gini.simpson - Gini-Simpson Index # norm.entropy - Normalized Entropy # pielou - Pielou’s Evenness # hill1, hill2, hill3 - Hill Numbers # 这个包也提供了查看 CDR3 氨基酸序列在每个位置的理化性质变化，与 positionalEntropy 类似，positionalProperty 的重要区别在于，它会删除 NA 值，因为这些值会使平均计算失效，并显示一个带有 95% 置信区间的带状结构，该区间围绕着所选特性的平均值 positionalProperty (combined.TCR[ c ( 1 , 2 )], chain = "TRB" , aa.length = 20 , method = "Atchley" ) # method——这里主要用于评估氨基酸的性质,提供不同的方法类型"Atchley", "Kidera", "stScales", "tScales", "VHSE" 方法名 核心内容 主要特征/维度 主要用途/来源 Atchley Factors 主成分分析提取的 5 个正交因子 (1) 极性/疏水性, (2) 二级结构倾向, (3) 分子体积, (4) 电荷性质, (5) 化学组成 Atchley et al., 2005，用于蛋白质序列的进化与结构分析 Kidera Factors 从 188 个性质中主成分分析得出的 10 个综合因子 综合了 10 种无相关性的物化特性 Kidera et al., 1985，用于蛋白质折叠与稳定性研究 MSWHIM 基于 WHIM（Weighted Holistic Invariant Molecular）描述符 从三维结构推导的分子形状、极性和质量分布 Todeschini et al., 化学信息学、分子建模 ProtFP Protein Fingerprint 编码氨基酸的疏水性、电子效应、立体结构参数 Van Westen et al., 蛋白质功能预测 stScales ST scales（Schneider & Wrede scales） 极性、疏水性、结构倾向 Schneider & Wrede, 1994 tScales T scale descriptors 化学拓扑指数，用图论描述氨基酸性质 Todeschini et al., 2000 VHSE Vectors of Hydrophobic, Steric, and Electronic properties 8 个综合因子：疏水性、立体效应、电子效应 Mei et al., 2005，用于 QSAR/蛋白质分析 # vizGenes是percentGeneUsage 的更灵活版本，主要作用是帮你看不同样本或条件下，免疫受体基因（V 基因、D 基因、J 基因等）的使用比例 vizGenes (combined.TCR, x.axis = "TRBV" , y.axis = NULL , # No specific y-axis variable, will group all samples plot = "barplot" ) # PercentGenes 函数是 percentGeneUsage 的专用别名，旨在量化特定免疫受体链中单个 V、D 或 J 基因座的使用。
默认情况下，它返回一个热图可视化 TRBV 基因片段的百分比，从而便于在样本间直观地比较基因使用谱 percentGenes (combined.TCR, chain = "TRB" , gene = "Vgene" ) # PercentGenes 返回的原始数据 (当 exportTable = TRUE 时) 可以作为进一步下游分析的强大输入，例如主成分分析等降维技术。
这允许你总结复杂的基因使用模式，并识别具有相似或不同库的样本 df.genes <- percentGenes (combined.TCR, chain = "TRB" , gene = "Vgene" , exportTable = TRUE ) # 生成基因矩阵进行pca降维 pc <- prcomp ( t (df.genes)) # 提取用于画图的数据 df_plot <- as.data.frame ( cbind (pc $ rotation[, 1 : 2 ], rownames (df.genes))) colnames (df_plot) <- c ( "PC1" , "PC2" , "Samples" ) df_plot $ PC1 <- as.numeric (df_plot $ PC1) df_plot $ PC2 <- as.numeric (df_plot $ PC2) ggplot (df_plot, aes ( x = PC1, y = PC2)) + geom_point ( aes ( fill = Samples), shape = 21 , size = 5 ) + guides ( fill= guide_legend ( title= "Samples" )) + scale_fill_manual ( values = hcl.colors ( nrow (df_plot), "inferno" )) + theme_classic + labs ( title = "PCA of TRBV Gene Usage" ) 可以看出P17B在V区域的特征较为独特 # 我们还可以可视化VJ基因的独特组合，查看不同组合出现的频率高低 percentVJ (combined.TCR[ 1 : 2 ], chain = "TRB" ) # 同样也可以通过降维查看不同样本的pca分布 df.vj <- percentVJ (combined.TCR, chain = "TRB" , exportTable = TRUE ) # Export proportions for PCA # 对V-J配对矩阵进行PCA降维 pc.vj <- prcomp ( t (df.vj)) # 提取画图数据 df_plot_vj <- as.data.frame ( cbind (pc.vj $ rotation[, 1 : 2 ], rownames (df.vj))) colnames (df_plot_vj) <- c ( "PC1" , "PC2" , "Sample" ) df_plot_vj $ PC1 <- as.numeric (df_plot_vj $ PC1) df_plot_vj $ PC2 <- as.numeric (df_plot_vj $ PC2) # 绘制PCA降维图 ggplot (df_plot_vj, aes ( x = PC1, y = PC2)) + geom_point ( aes ( fill = Sample), shape = 21 , size = 5 ) + guides ( fill= guide_legend ( title= "Samples" )) + scale_fill_manual ( values = hcl.colors ( nrow (df_plot_vj), "inferno" )) + theme_classic + labs ( title = "PCA of TRBV-TRBJ Gene Pairings" ) 可以看出P17B的VJ组合特征较为独特 基于单细胞TCR/BCR数据，我们也可以基于序列相似性对克隆类型做类似scRNA-seq的聚类分析 # 基于TRA链的氨基酸序列对前两个样本进行聚类 sub_combined <- clonalCluster (combined.TCR[ c ( 1 , 2 )], chain = "TRA" , sequence = "aa" , threshold = 0.85 ) <1 表示归一化相似度，而值>=1 表示原始编辑距离 # 聚类过程遵循以下关键步骤：# 序列选择：该函数选择特定链的核苷酸 (序列 =“nt”) 或氨基酸 (序列 =“aa”) CDR3 序列。
# 距离计算：它计算所有序列对之间的编辑距离。默认情况下，它还要求序列共享相同的 V 基因 (use.v = TRUE)。# 网络构建：在满足相似度阈值的任意两个序列之间创建一条边，形成网络图。# 聚类：一种基于图的聚类算法用于识别网络中的连通组件或社区。默认情况下，它将所有直接或间接连通的序列识别为单个聚类 (cluster.method =“components”)。
# 输出：该函数可以将结果簇 ID 添加到输入对象中，返回一个用于网络分析的图表对象，或导出一个稀疏邻接矩阵。
head (sub_combined[[ 1 ]][, c ( "barcode" , "TCR1" , "TRA_cluster" )]) # barcode TCR1 TRA_cluster # 1 P17B_AAACCTGAGTACGACG-1 TRAV25.TRAJ20.TRAC <NA> # 2 P17B_AAACCTGCAACACGCC-1 TRAV38-2/DV8.TRAJ52.TRAC <NA> # 3 P17B_AAACCTGCAGGCGATA-1 TRAV12-1.TRAJ9.TRAC TRA:Cluster.7 # 4 P17B_AAACCTGCATGAGCGA-1 TRAV12-1.TRAJ9.TRAC TRA:Cluster.7 # 5 P17B_AAACGGGAGAGCCCAA-1 TRAV20.TRAJ8.TRAC <NA> # 6 P17B_AAACGGGAGCGTTTAC-1 TRAV12-1.TRAJ9.TRAC TRA:Cluster.7 2.4 比较样本克隆多样性、差异与相似性 clonalRarefaction 函数基于克隆丰度数据，利用 Hill 数估计稀释曲线（rarefaction），即推断克隆群的物种丰富度和多样性。
Hill 数指标包括：- 0：物种丰富度（克隆数量） - 1：香农多样性指数 - 2：辛普森多样性指数 函数支持三种绘图类型：1. 基于样本量的稀释/外推曲线 2. 样本完整性曲线（采样覆盖度） 3. 基于覆盖度的稀释/外推曲线 主要参数有：- plot.type：选择绘图类型 - hill.numbers：指定Hill数指标（0、1或2） - n.boots：bootstrap重采样次数，用于置信区间估计，默认20次 示例代码（使用物种丰富度指标，绘制基于样本量的稀释曲线）：clonalRarefaction (combined.TCR, plot.type = 1 , # 绘图类型 hill.numbers = 0 , # Hill指数 n.boots = 2 ) # 置信区间 clonalRarefaction (combined.TCR, plot.type = 2 , hill.numbers = 0 , n.boots = 2 ) clonalRarefaction 提供了一种复杂的免疫库组成建模方法，可以区分稀有克隆和扩增克隆的分布。
通过应用拼接统计模型，它提供了库的基本克隆结构的更准确表示，使得可靠的比较和对免疫系统动力学的更深入的理解成为可能。软件还提供了 clonalSizeDistribution 是用来分析免疫克隆大小分布的一种高级方法。它基于Koch等人提出的统计模型，把大克隆和小克隆分开建模：大克隆符合幂律分布（也就是少数大克隆占优势），小克隆则符合泊松分布（数量多但个体小）。通过拟合这个模型，可以用数学距离（欧氏距离）比较不同样本的克隆分布差异。
主要参数是聚类的方法，比如常用的“ward.D2”，用来把克隆分组，方便比较。举个例子，下面代码用氨基酸克隆信息，采用“ward.D2”方法来做分析：clonalSizeDistribution (combined.TCR, cloneCall = "aa" , method = "ward.D2" ) 同样软件也提供比较样本间相似性的方法，clonalOverlap 用于衡量和可视化不同样本之间克隆序列的相似度或重叠程度。
通过不同的方法参数，可以选择不同的相似度指标：overlap：重叠系数 morisita：Morisita重叠指数 jaccard：Jaccard指数 cosine：余弦相似度 raw：重叠克隆的精确数量 clonalOverlap (combined.TCR, cloneCall = "strict" , method = "morisita" ) CloneOverlap 是一种评估不同免疫库之间相似性和共享克隆类型的工具。
通过提供各种定量方法和清晰的热图可视化，它使研究人员能够识别重叠程度，为共享免疫反应、交叉反应性或不同处理对克隆组成的影响提供见解。2.5 合并克隆和单细胞对象 接下来我们可以将scTCR数据整合到scRNA-seq数据中, 在单细胞 RNA 测序工作流程中，通常首先识别高度可变的特征来进行降维。然后将这些特征直接用于 UMAP/tSNE 投影或作为主成分分析的输入。同样的方法通常也应用于聚类。
然而，在免疫聚焦的数据集中，来自 TCR 和 BCR 的 VDJ 基因通常是最可变的基因之一。这种变异性是由于淋巴细胞内的克隆扩张和多样性而自然产生的。因此，UMAP 预测和聚类结果可能受到克隆信息的影响，而不是细胞类型之间更广泛的转录差异。为了缓解这个问题，一个常见的策略是在进行聚类和降维之前，从高度可变的特征集中排除 VDJ 基因。
我们引入了一组功能，通过从 Seurat Object 或基因名称向量中移除 VDJ 相关基因来促进这个过程 (对基于 SCE 的工作流程有用) # 导入全部的示例数据 scRep_example <- readRDS ( "/path/to/project" ) # 计算高变基因 scRep_example <- FindVariableFeatures (scRep_example, selection.method = "vst" , nfeatures = 2000 ) # 可以使用如下函数去除VDJ基因 # quietVDJgenes – Removes both TCR and BCR VDJ genes. # quietTCRgenes – Removes only TCR VDJ genes. # quietBCRgenes – Removes only BCR VDJ genes, but retains BCR VDJ pseudogenes in the variable features. # 先检查scRNA-seq的前10个高变基因 VariableFeatures (scRep_example)[ 1 : 10 ] # [1] "S100A8" "S100A9" "LYZ" "S100A12" "CXCL8" "IL22" # [7] "CCL4L2" "IGLV3-19" "CTSL" "IGKV1-5" 接下来移除 TCR VDJ 基因 library (Trex) # 去除 TCR VDJ 基因 scRep_example <- quietTCRgenes (scRep_example) # 检查去除TCR基因后的前10个高变基因 VariableFeatures (scRep_example)[ 1 : 10 ] # [1] "S100A8" "S100A9" "LYZ" "S100A12" "CXCL8" "IL22" # [7] "CCL4L2" "IGLV3-19" "CTSL" "IGKV1-5" 使用包含 30,000 个细胞的完整单细胞对象。
我们将使用 Seurat 和单细胞实验 (SCE) 与散点图进行串联可视化 # 将seurat对象转换为sce对象 sce <- Seurat :: as.SingleCellExperiment (scRep_example) 在通过 combineBCR 或 combineTCR 将 contig 数据处理成克隆后，我们可以使用 combineExpression 将克隆信息添加到单单元格对象中 重要的是，附加的主要要求是匹配 Seurat 或 Single-Cell Experiment 对象元数据的行名称中的连续单元条形码和条形码。
如果这些不匹配，附加将失败。我们建议更改 Single-Cell Experiment 对象的条形码，以便于使用。
# 我们可以使用刚创建的 Single-Cell Experiment 对象来查看默认的 cloneSize 分组，# 其 group.by 设置为 combineTCR 中使用的 sample 变量 scRep_example <- combineExpression (combined.TCR, scRep_example, cloneCall= "gene" , group.by = "sample" , proportion = FALSE , cloneSize= c ( Single= 1 , Small= 5 , Medium= 20 , Large= 100 , Hyperexpanded= 500 )) # 构建绘制颜色 colorblind_vector <- hcl.colors ( n= 7 , palette = "inferno" , fixup = TRUE ) Seurat :: DimPlot (scRep_example, group.by = "cloneSize" ) + scale_color_manual ( values= rev (colorblind_vector[ c ( 1 , 3 , 4 , 5 , 7 )])) 如果你同时拥有 TCR 和 BCR 富集数据，或者希望同时包含 γ-delta 和 α-beta T 单元的信息，可以创建一个包含 combineTCR 和 combineBCR 输出的单个列表，然后使用 combineExpression TCR <- combineTCR (...) BCR <- combineBCR (...) list.receptors <- c (TCR, BCR) seurat <- combineExpression (list.receptors, seurat, cloneCall= "gene" , proportion = TRUE ) 2.6 整合后数据可视化 整合scTCR数据和scRNA-seq数据的可视化,以降维图为参考，clonalOverlay 生成克隆扩展细胞位置的叠加。
它突出显示 UMAP 或 tSNE 图中高克隆频率或比例的区域 clonalOverlay (scRep_example, reduction = "umap" , cutpoint = 1 , # 以1个作为划分克隆和未克隆的区别 bins = 10 , # 把克隆扩增数分成多少个区间 facet.by = "orig.ident" ) + guides ( color = "none" ) 类似于 clonalOverlay，clonalNetwork 用来展示不同细胞群之间克隆的网络互动 library (ggraph) # # Attaching package: 'ggraph' # The following object is masked from 'package:sp': # # geometry # 克隆关系网络图，可视化不同细胞群之间共享 TCR/BCR 克隆的程度 clonalNetwork (scRep_example, reduction = "umap" , group.by = "seurat_clusters" , filter.clones = NULL , filter.identity = NULL , cloneCall = "aa" ) filter.identity 参数来查看相对于单个簇的克隆关系 Identity filter clonalNetwork (scRep_example, reduction = "umap" , group.by = "seurat_clusters" , filter.identity = 3 , cloneCall = "aa" ) 还可以高亮特定的克隆类型在UMAP图中显示 scRep_example <- highlightClones (scRep_example, cloneCall= "aa" , sequence = c ( "CAERGSGGSYIPTF_CASSDPSGRQGPRWDTQYF" , "CARKVRDSSYKLIF_CASSDSGYNEQFF" )) Seurat :: DimPlot (scRep_example, group.by = "highlight" ) + guides ( color= guide_legend ( nrow= 3 , byrow= TRUE )) + ggplot2 :: theme ( plot.title = element_blank , legend.position = "bottom" ) CloneOccupy 按簇可视化单元格的数量，并将其分类为特定的克隆频率范围。
clonalOccupy (scRep_example, x.axis = "seurat_clusters" ) 这个包还提供了alluvialClones来帮助查看跨多个分类变量的克隆，从而检查这些变量之间的交换 alluvialClones (scRep_example, cloneCall = "aa" , y.axes = c ( "orig.ident" , "seurat_clusters" ), color = c ( "CVVSDNTGGFKTIF_CASSVRRERANTGELFF" , "NA_CASSVRRERANTGELFF" )) + scale_fill_manual ( values = c ( "grey" , colorblind_vector[ 3 ])) # Scale for fill is already present. # Adding another scale for fill, which will replace the existing scale. 像冲积图一样，我们也可以使用 circlize R 包中的和弦图来可视化簇之间的相互连接 library (circlize) # ======================================== # circlize version 0.4.16 # CRAN page: https://cran.r-project.org/package=circlize # Github page: https://github.com/jokergoo/circlize # Documentation: https://jokergoo.github.io/circlize_book/book/ # # If you use it in published research, please cite: # Gu, Z. circlize implements and enhances circular visualization # in R. Bioinformatics 2014. # # This message can be suppressed by: # suppressPackageStartupMessages(library(circlize)) # ======================================== library (scales) # 提取克隆共享矩阵信息 circles <- getCirclize (scRep_example, group.by = "seurat_clusters" ) # 给每个cluster赋予颜色 grid.cols <- hue_pal ( length ( unique (scRep_example $ seurat_clusters))) names (grid.cols) <- unique (scRep_example $ seurat_clusters) # 绘制弦图 chordDiagram (circles, self.link = 1 , grid.col = grid.cols) 量化克隆偏差 指数 全称及含义 高值代表的生物学含义 低值代表的生物学含义 expa Clonal Expansion（克隆扩增） 衡量一个细胞群中是否存在少数克隆占据多数细胞的现象。
计算方式为 \(1 - \frac{H}{H_{max}}\)，其中 \(H\) 为 Shannon 熵。熵低 → 集中度高 → 扩增指数高。少数高频克隆主导群体，提示可能存在 抗原驱动的免疫应答 或 肿瘤特异性反应。克隆分布均匀，没有明显扩增，提示群体多样性高，缺乏明显抗原选择。
migr Cross-tissue Migration（跨组织迁移） 衡量同一克隆型（clonotype）在不同组织（由 type 参数定义）中的分布均匀程度。基于熵计算。克隆在多个组织均有分布，提示可能存在 免疫细胞在组织间的迁移 （如外周血 ↔︎ 肿瘤）。克隆集中于单一组织，提示克隆 局部驻留 或缺乏迁移。
tran State Transition（状态转变） 衡量同一克隆型在不同功能状态（由 cluster 定义，如效应型、记忆型、耗竭型）中的分布程度。基于熵计算。克隆在多个功能状态间分布，提示可能发生 功能状态转换 （如记忆 T 细胞激活为效应 T 细胞）。克隆主要存在于单一功能状态，提示功能状态稳定，转变较少。
参考文献 Zhang et al., 2018, Nature 2.7 量化克隆偏差 # 计算和绘制 STARTRAC 的三种参数 StartracDiversity (scRep_example, type = "seurat_clusters" , group.by = "orig.ident" ) 最后 clonalBias 也可以用来量化单个克隆在不同细胞群或功能区间中的偏向性的新指标，衡量某个克隆是否倾向于集中在一个特定的细胞群（cluster）或功能区间（compartment） # 这里实例数据中存在有的样本在某些clusters中数量为0，导致这步报错，因此我们挑选两个不为0的样本进行测试 library (dplyr) # # Attaching package: 'dplyr' # The following objects are masked from 'package:stats': # # filter, lag # The following objects are masked from 'package:base': # # intersect, setdiff, setequal, union # 先筛选元数据，只保留选中组 Idents (scRep_example) <- "orig.ident" selected_ids <- c ( "P17L" , "P18L" ) # 子集数据，只保留对应细胞 sc <- subset (scRep_example, idents= selected_ids) # 再做clonalBias (sc, cloneCall = "aa" , split.by = "orig.ident" , group.by = "seurat_clusters" , n.boots = 10 , min.expand = 5 ) # Smoothing formula not specified. Using: y ~ qss(x, lambda = 3) # Warning in rq.fit.sfn(x, y, tau = tau, rhs = rhs, control = control, ...): tiny diagonals replaced with Inf when calling blkfct 2.8 基于scBCR-seq数据的分析 基本分析思路与scTCR-seq数据一致，这里快速展现一下代码和结果 载入 scRNA 数据（Seurat） # 读取表达矩阵 seurat_obj <- Read10X_h5 ( "./dataset-multi-practice/runs/HumanB_Cell_multi/outs/per_sample_outs/HumanB_Cell_multi/count/sample_filtered_feature_bc_matrix.h5" ) seurat_obj <- CreateSeuratObject ( counts = seurat_obj, project = "scRNA_VDJ" ) # 后续常规预处理 seurat_obj <- NormalizeData (seurat_obj) seurat_obj <- FindVariableFeatures (seurat_obj) seurat_obj <- ScaleData (seurat_obj) seurat_obj <- RunPCA (seurat_obj) pct <- seurat_obj[[ "pca" ]] @ stdev / sum (seurat_obj[[ "pca" ]] @ stdev) * 100 cumu <- cumsum (pct) co1 <- which (cumu > 90 & pct < 5 )[ 1 ] co2 <- sort ( which ((pct[ 1 : length (pct) - 1 ] - pct[ 2 : length (pct)]) > 0.1 ), decreasing = TRUE )[ 1 ] + 1 pcs <- min (co1, co2) bestpc <- 1 : pcs seurat_obj <- RunUMAP (seurat_obj, dims = bestpc) seurat_obj <- FindNeighbors (seurat_obj, reduction = "pca" , dims = bestpc) seurat_obj <- FindClusters (seurat_obj, resolution = 0.3 ) # Modularity Optimizer version 1.3.0 by Ludo Waltmand Nees Jan van Eck # # Number of nodes: 1149 # Number of edges: 34932 # # Running Louvain algorithm... # Maximum modularity in 10 random starts: 0.8555 # Number of communities: 4 # Elapsed time: 0 seconds 载入 V(D)J 数据（scRepertoire） # 读取VDJ数据 contig <- read.csv ( "./dataset-multi-practice/runs/HumanB_Cell_multi/outs/per_sample_outs/HumanB_Cell_multi/vdj_b/filtered_contig_annotations.csv" ) # 整合 Seurat 对象与克隆信息 combined <- combineBCR (contig, samples = "sample1" , ID = "BCR" ) # 或 "T-AB" 用于TCR # 添加到 Seurat 对象中（使用cell barcode匹配） # 因为combined对象的barcode列是由"samples_ID_barcode"构成 colnames (seurat_obj) <- paste0 ( "sample1_BCR_" , colnames (seurat_obj)) rownames (seurat_obj @ meta.data) <- colnames (seurat_obj) seurat_obj $ barcode <- colnames (seurat_obj) seurat_obj <- combineExpression (combined, seurat_obj, cloneCall = "strict" ) # BCR(TCR)信息被加入到seurat对象meta信息中 head (seurat_obj @ meta.data) # orig.ident nCount_RNA nFeature_RNA # sample1_BCR_AAACCTGAGGGCTCTC-1 scRNA_VDJ 4424 1848 # sample1_BCR_AAACCTGCAAATTGCC-1 scRNA_VDJ 1540 949 # sample1_BCR_AAACCTGGTAAGGATT-1 scRNA_VDJ 5404 2236 # sample1_BCR_AAACCTGGTAATAGCA-1 scRNA_VDJ 5188 2147 # sample1_BCR_AAACCTGGTACGCACC-1 scRNA_VDJ 5168 1983 # sample1_BCR_AAACCTGTCCAACCAA-1 scRNA_VDJ 10327 3424 # RNA_snn_res.0.3 seurat_clusters # sample1_BCR_AAACCTGAGGGCTCTC-1 0 0 # sample1_BCR_AAACCTGCAAATTGCC-1 3 3 # sample1_BCR_AAACCTGGTAAGGATT-1 2 2 # sample1_BCR_AAACCTGGTAATAGCA-1 0 0 # sample1_BCR_AAACCTGGTACGCACC-1 0 0 # sample1_BCR_AAACCTGTCCAACCAA-1 0 0 # barcode # sample1_BCR_AAACCTGAGGGCTCTC-1 sample1_BCR_AAACCTGAGGGCTCTC-1 # sample1_BCR_AAACCTGCAAATTGCC-1 sample1_BCR_AAACCTGCAAATTGCC-1 # sample1_BCR_AAACCTGGTAAGGATT-1 sample1_BCR_AAACCTGGTAAGGATT-1 # sample1_BCR_AAACCTGGTAATAGCA-1 sample1_BCR_AAACCTGGTAATAGCA-1 # sample1_BCR_AAACCTGGTACGCACC-1 sample1_BCR_AAACCTGGTACGCACC-1 # sample1_BCR_AAACCTGTCCAACCAA-1 sample1_BCR_AAACCTGTCCAACCAA-1 # CTgene # sample1_BCR_AAACCTGAGGGCTCTC-1 IGHV4-61.NA.IGHJ5.IGHM_IGLV1-40.IGLJ2.IGLC3 # sample1_BCR_AAACCTGCAAATTGCC-1 <NA> # sample1_BCR_AAACCTGGTAAGGATT-1 IGHV1-2.NA.IGHJ5.IGHM_IGKV4-1.IGKJ2.IGKC # sample1_BCR_AAACCTGGTAATAGCA-1 IGHV1-18.IGHD6-13.IGHJ5.IGHM_IGLV1-44.IGLJ3.IGLC2 # sample1_BCR_AAACCTGGTACGCACC-1 IGHV7-4-1.NA.IGHJ4.IGHM_IGLV7-46.IGLJ2.IGLC2 # sample1_BCR_AAACCTGTCCAACCAA-1 IGHV3-23.NA.IGHJ4.IGHM_IGKV4-1.IGKJ1.IGKC # CTnt # sample1_BCR_AAACCTGAGGGCTCTC-1 TGTGCGAGAGGGGATAGCAGTGGCTGGCGAGGGGGCAACTGGTTCGACCCCTGG_TGCCAGTCCTATGACAGCAGCCTGAGTGACGTGTTC # sample1_BCR_AAACCTGCAAATTGCC-1 <NA> # sample1_BCR_AAACCTGGTAAGGATT-1 TGTGCGATGGGATATTGTATTAATAATAACTGTTACGAGGGGTGGTTCGACCCCTGG_TGTCAGCAATACTATGATACTCCCCGTACTTTT # sample1_BCR_AAACCTGGTAATAGCA-1 TGTGCGAGAGCCAAACGTTGGGGGTATAGCAGCAGCTGGTGCGACTACTGG_TGTGCAGCATGGGATGACAGCCTGAATGGTGGGGTGTTC # sample1_BCR_AAACCTGGTACGCACC-1 TGTGCGAGAGCCCTGGGGGCTATTGAACTCTTTGATTACTGG_TGCTTGCTCTCCTATAGTGGTGCCCATGTGGTATTC # sample1_BCR_AAACCTGTCCAACCAA-1 TGTGCGAAAGAAATGGTTCGGGGGTGGGTTGACTACTGG_TGTCAGCAATATTATAGTACTCGGCGGACGTTC # CTaa # sample1_BCR_AAACCTGAGGGCTCTC-1 CARGDSSGWRGGNWFDPW_CQSYDSSLSDVF # sample1_BCR_AAACCTGCAAATTGCC-1 <NA> # sample1_BCR_AAACCTGGTAAGGATT-1 CAMGYCINNNCYEGWFDPW_CQQYYDTPRTF # sample1_BCR_AAACCTGGTAATAGCA-1 CARAKRWGYSSSWCDYW_CAAWDDSLNGGVF # sample1_BCR_AAACCTGGTACGCACC-1 CARALGAIELFDYW_CLLSYSGAHVVF # sample1_BCR_AAACCTGTCCAACCAA-1 CAKEMVRGWVDYW_CQQYYSTRRTF # CTstrict # sample1_BCR_AAACCTGAGGGCTCTC-1 IGH.1.IGHV4-61_IGLC:Cluster.21.IGLV1-40 # sample1_BCR_AAACCTGCAAATTGCC-1 <NA> # sample1_BCR_AAACCTGGTAAGGATT-1 IGH.2.IGHV1-2_IGLC.1.IGKV4-1 # sample1_BCR_AAACCTGGTAATAGCA-1 IGH.3.IGHV1-18_IGLC:Cluster.22.IGLV1-44 # sample1_BCR_AAACCTGGTACGCACC-1 IGH.4.IGHV7-4-1_IGLC:Cluster.40.IGLV7-46 # sample1_BCR_AAACCTGTCCAACCAA-1 IGH.5.IGHV3-23_IGLC:Cluster.19.IGKV4-1 # clonalProportion clonalFrequency # sample1_BCR_AAACCTGAGGGCTCTC-1 0.0009025271 1 # sample1_BCR_AAACCTGCAAATTGCC-1 NA # sample1_BCR_AAACCTGGTAAGGATT-1 0.0009025271 1 # sample1_BCR_AAACCTGGTAATAGCA-1 0.0009025271 1 # sample1_BCR_AAACCTGGTACGCACC-1 0.0009025271 1 # sample1_BCR_AAACCTGTCCAACCAA-1 0.0009025271 1 # cloneSize # sample1_BCR_AAACCTGAGGGCTCTC-1 Small (1e-04 < X <= 0.001) # sample1_BCR_AAACCTGCAAATTGCC-1 <NA> # sample1_BCR_AAACCTGGTAAGGATT-1 Small (1e-04 < X <= 0.001) # sample1_BCR_AAACCTGGTAATAGCA-1 Small (1e-04 < X <= 0.001) # sample1_BCR_AAACCTGGTACGCACC-1 Small (1e-04 < X <= 0.001) # sample1_BCR_AAACCTGTCCAACCAA-1 Small (1e-04 < X <= 0.001) # 因为测试数据只有一个样本，但是我们可以模拟出两个样本来后续画图 n <- ncol (seurat_obj) # 细胞总数 seurat_obj $ sample <- NA # 新建一列 seurat_obj $ sample[ 1 : floor (n / 2 )] <- "s1" seurat_obj $ sample[( floor (n / 2 ) + 1 ) : n] <- "s2" head (seurat_obj @ meta.data) # orig.ident nCount_RNA nFeature_RNA # sample1_BCR_AAACCTGAGGGCTCTC-1 scRNA_VDJ 4424 1848 # sample1_BCR_AAACCTGCAAATTGCC-1 scRNA_VDJ 1540 949 # sample1_BCR_AAACCTGGTAAGGATT-1 scRNA_VDJ 5404 2236 # sample1_BCR_AAACCTGGTAATAGCA-1 scRNA_VDJ 5188 2147 # sample1_BCR_AAACCTGGTACGCACC-1 scRNA_VDJ 5168 1983 # sample1_BCR_AAACCTGTCCAACCAA-1 scRNA_VDJ 10327 3424 # RNA_snn_res.0.3 seurat_clusters # sample1_BCR_AAACCTGAGGGCTCTC-1 0 0 # sample1_BCR_AAACCTGCAAATTGCC-1 3 3 # sample1_BCR_AAACCTGGTAAGGATT-1 2 2 # sample1_BCR_AAACCTGGTAATAGCA-1 0 0 # sample1_BCR_AAACCTGGTACGCACC-1 0 0 # sample1_BCR_AAACCTGTCCAACCAA-1 0 0 # barcode # sample1_BCR_AAACCTGAGGGCTCTC-1 sample1_BCR_AAACCTGAGGGCTCTC-1 # sample1_BCR_AAACCTGCAAATTGCC-1 sample1_BCR_AAACCTGCAAATTGCC-1 # sample1_BCR_AAACCTGGTAAGGATT-1 sample1_BCR_AAACCTGGTAAGGATT-1 # sample1_BCR_AAACCTGGTAATAGCA-1 sample1_BCR_AAACCTGGTAATAGCA-1 # sample1_BCR_AAACCTGGTACGCACC-1 sample1_BCR_AAACCTGGTACGCACC-1 # sample1_BCR_AAACCTGTCCAACCAA-1 sample1_BCR_AAACCTGTCCAACCAA-1 # CTgene # sample1_BCR_AAACCTGAGGGCTCTC-1 IGHV4-61.NA.IGHJ5.IGHM_IGLV1-40.IGLJ2.IGLC3 # sample1_BCR_AAACCTGCAAATTGCC-1 <NA> # sample1_BCR_AAACCTGGTAAGGATT-1 IGHV1-2.NA.IGHJ5.IGHM_IGKV4-1.IGKJ2.IGKC # sample1_BCR_AAACCTGGTAATAGCA-1 IGHV1-18.IGHD6-13.IGHJ5.IGHM_IGLV1-44.IGLJ3.IGLC2 # sample1_BCR_AAACCTGGTACGCACC-1 IGHV7-4-1.NA.IGHJ4.IGHM_IGLV7-46.IGLJ2.IGLC2 # sample1_BCR_AAACCTGTCCAACCAA-1 IGHV3-23.NA.IGHJ4.IGHM_IGKV4-1.IGKJ1.IGKC # CTnt # sample1_BCR_AAACCTGAGGGCTCTC-1 TGTGCGAGAGGGGATAGCAGTGGCTGGCGAGGGGGCAACTGGTTCGACCCCTGG_TGCCAGTCCTATGACAGCAGCCTGAGTGACGTGTTC # sample1_BCR_AAACCTGCAAATTGCC-1 <NA> # sample1_BCR_AAACCTGGTAAGGATT-1 TGTGCGATGGGATATTGTATTAATAATAACTGTTACGAGGGGTGGTTCGACCCCTGG_TGTCAGCAATACTATGATACTCCCCGTACTTTT # sample1_BCR_AAACCTGGTAATAGCA-1 TGTGCGAGAGCCAAACGTTGGGGGTATAGCAGCAGCTGGTGCGACTACTGG_TGTGCAGCATGGGATGACAGCCTGAATGGTGGGGTGTTC # sample1_BCR_AAACCTGGTACGCACC-1 TGTGCGAGAGCCCTGGGGGCTATTGAACTCTTTGATTACTGG_TGCTTGCTCTCCTATAGTGGTGCCCATGTGGTATTC # sample1_BCR_AAACCTGTCCAACCAA-1 TGTGCGAAAGAAATGGTTCGGGGGTGGGTTGACTACTGG_TGTCAGCAATATTATAGTACTCGGCGGACGTTC # CTaa # sample1_BCR_AAACCTGAGGGCTCTC-1 CARGDSSGWRGGNWFDPW_CQSYDSSLSDVF # sample1_BCR_AAACCTGCAAATTGCC-1 <NA> # sample1_BCR_AAACCTGGTAAGGATT-1 CAMGYCINNNCYEGWFDPW_CQQYYDTPRTF # sample1_BCR_AAACCTGGTAATAGCA-1 CARAKRWGYSSSWCDYW_CAAWDDSLNGGVF # sample1_BCR_AAACCTGGTACGCACC-1 CARALGAIELFDYW_CLLSYSGAHVVF # sample1_BCR_AAACCTGTCCAACCAA-1 CAKEMVRGWVDYW_CQQYYSTRRTF # CTstrict # sample1_BCR_AAACCTGAGGGCTCTC-1 IGH.1.IGHV4-61_IGLC:Cluster.21.IGLV1-40 # sample1_BCR_AAACCTGCAAATTGCC-1 <NA> # sample1_BCR_AAACCTGGTAAGGATT-1 IGH.2.IGHV1-2_IGLC.1.IGKV4-1 # sample1_BCR_AAACCTGGTAATAGCA-1 IGH.3.IGHV1-18_IGLC:Cluster.22.IGLV1-44 # sample1_BCR_AAACCTGGTACGCACC-1 IGH.4.IGHV7-4-1_IGLC:Cluster.40.IGLV7-46 # sample1_BCR_AAACCTGTCCAACCAA-1 IGH.5.IGHV3-23_IGLC:Cluster.19.IGKV4-1 # clonalProportion clonalFrequency # sample1_BCR_AAACCTGAGGGCTCTC-1 0.0009025271 1 # sample1_BCR_AAACCTGCAAATTGCC-1 NA # sample1_BCR_AAACCTGGTAAGGATT-1 0.0009025271 1 # sample1_BCR_AAACCTGGTAATAGCA-1 0.0009025271 1 # sample1_BCR_AAACCTGGTACGCACC-1 0.0009025271 1 # sample1_BCR_AAACCTGTCCAACCAA-1 0.0009025271 1 # cloneSize sample # sample1_BCR_AAACCTGAGGGCTCTC-1 Small (1e-04 < X <= 0.001) s1 # sample1_BCR_AAACCTGCAAATTGCC-1 <NA> s1 # sample1_BCR_AAACCTGGTAAGGATT-1 Small (1e-04 < X <= 0.001) s1 # sample1_BCR_AAACCTGGTAATAGCA-1 Small (1e-04 < X <= 0.001) s1 # sample1_BCR_AAACCTGGTACGCACC-1 Small (1e-04 < X <= 0.001) s1 # sample1_BCR_AAACCTGTCCAACCAA-1 Small (1e-04 < X <= 0.001) s1 基本克隆可视化 # 计算和量化细胞克隆（clone）多样性和丰度，帮助你了解免疫细胞中BCR克隆的分布情况 clonalQuant (combined, cloneCall= "strict" , chain = "both" , group.by = "sample" , scale = TRUE ) # 计算并展示免疫细胞克隆（BCR）在样本中绝对或相对的丰度分布情况 clonalAbundance (combined, cloneCall = "gene" , scale = FALSE ) clonalAbundance (combined, cloneCall = "gene" , scale = TRUE ) # 可视化 CDR3 序列的长度分布 clonalLength (combined, cloneCall= "aa" , acid sequence）定义克隆，chain = "both" ) clonalCompare (combined, top.clones = 10 , samples = c ( "s1" , "s2" ), cloneCall= "aa" , graph = "alluvial" ) # 以密度图的形式进行缩放 clonalLength (combined, cloneCall= "aa" , acid sequence）定义克隆，chain = "IGH" , scale = TRUE ) # 接下来我们可以查看不同克隆在整个免疫受体库中所占比例的分布情况 clonalHomeostasis (combined, cloneCall = "gene" ) # 也可以通过改变截点来改变分布比例 clonalHomeostasis (combined, cloneCall = "gene" , cloneSize = c ( Rare = 0.001 , Small = 0.01 , Medium = 0.1 , Large = 0.3 , Hyperexpanded = 1 )) # 看所有样本中不同排名区间的克隆类型在不同样本中所占比例 clonalProportion (combined, cloneCall = "gene" ) # 看所有样本中不同排名区间的克隆类型在不同样本中所占比例 clonalProportion (combined, cloneCall = "nt" , clonalSplit = c ( 1 , 5 , 10 , 100 , 1000 , 10000 )) # CDR3 氨基酸位置组成百分比的热图 # 在 CDR3 序列的第 1、2、3、…、20 位上，分别有哪些氨基酸出现，以及它们出现的比例是多少 percentAA (combined, chain = "both" , aa.length = 20 ) # 我们还可以量化 CDR3 序列中氨基酸残基的熵 / 多样性水平。
positionalEntropy 将百分比 AA 的残基量化与多样性计算相结合。
0, 用于比较同一位置不同样本的氨基酸多样性 positionalEntropy (combined, chain = "IGH" , method = "norm.entropy" , aa.length = 20 ) # positionalProperty类似positionalEntropy positionalProperty (combined[ 1 ], chain = "IGH" , aa.length = 20 , method = "Atchley" ) # 看不同样本或条件下，免疫受体基因（V 基因、D 基因、J 基因等）的使用比例 vizGenes (combined, x.axis = "IGLC" , y.axis = NULL , # No specific y-axis variable, will group all samples plot = "barplot" ) # 量化特定免疫受体链中单个 V、D 或 J 基因座的使用 percentGenes (combined, chain = "IGH" , gene = "Vgene" ) # 可视化VJ基因的独特组合，查看不同组合出现的频率高低 percentVJ (combined[ 1 ], chain = "IGH" ) # 基于IGH链的氨基酸序列进行聚类 sub_combined <- clonalCluster (combined, chain = "IGL" , sequence = "aa" , threshold = 0.85 ) <1 表示归一化相似度，而值>=1 表示原始编辑距离 head (sub_combined[[ 1 ]][, c ( "barcode" , "IGLC" , "IGL_cluster" )]) # barcode IGLC IGL_cluster # 1 sample1_BCR_TGCCCATGTGTGACGA-1 IGLV2-14.IGLJ2.IGLC2 <NA> # 2 sample1_BCR_CAGTCCTGTCCGAGTC-1 IGLV2-8.IGLJ3.IGLC2 <NA> # 3 sample1_BCR_TGTATTCAGAACAATC-1 IGLV2-8.IGLJ1.IGLC1 IGL:Cluster.7 # 4 sample1_BCR_AGCTCTCCACATAACC-1 IGLV2-14.IGLJ2.IGLC2 <NA> # 5 sample1_BCR_CGTCTACCAAGAGGCT-1 IGLV2-14.IGLJ1.IGLC1 <NA> # 6 sample1_BCR_TGACTTTAGGATCGCA-1 IGLV2-14.IGLJ1.IGLC1 IGL:Cluster.9 # 查看不同克隆类型在scRNA-seq的细胞umap中分布 seurat_obj # An object of class Seurat # 36601 features across 1149 samples within 1 assay # Active assay: RNA (36601 features, 2000 variable features) # 3 layers present: counts, data, scale.data # 2 dimensional reductions calculated: pca, umap # 构建绘制颜色 colorblind_vector <- hcl.colors ( n= 7 , palette = "inferno" , fixup = TRUE ) Seurat :: DimPlot (seurat_obj, group.by = "cloneSize" ) + scale_color_manual ( values= rev (colorblind_vector[ c ( 1 , 3 , 4 , 5 , 7 )])) # clonalOverlay 生成克隆扩展细胞位置的叠加。
它突出显示 UMAP 或 tSNE 图中高克隆频率或比例的区域 clonalOverlay (seurat_obj, reduction = "umap" , cutpoint = 1 , # 以1个作为划分克隆和未克隆的区别 bins = 10 , # 把克隆扩增数分成多少个区间 facet.by = "sample" ) + guides ( color = "none" ) # clonalNetwork用来展示不同细胞群之间克隆的网络互动 clonalNetwork (seurat_obj, reduction = "umap" , group.by = "seurat_clusters" , filter.clones = NULL , filter.identity = NULL , cloneCall = "aa" ) # 可以使用 filter.identity 参数来查看相对于单个簇的克隆关系 clonalNetwork (seurat_obj, reduction = "umap" , group.by = "seurat_clusters" , filter.identity = 3 , cloneCall = "aa" ) # 高亮特定的克隆类型在UMAP图中显示 seurat_obj <- highlightClones (seurat_obj, cloneCall= "aa" , sequence = c ( "CAMGYCINNNCYEGWFDPW_CQQYYDTPRTF" , "CARALGAIELFDYW_CLLSYSGAHVVF" )) Seurat :: DimPlot (seurat_obj, group.by = "highlight" ) + guides ( color= guide_legend ( nrow= 3 , byrow= TRUE )) + ggplot2 :: theme ( plot.title = element_blank , legend.position = "bottom" ) # 按簇可视化单元格的数量，并将其分类为特定的克隆频率范围 clonalOccupy (seurat_obj, x.axis = "seurat_clusters" ) # 帮助查看跨多个分类变量的克隆，从而检查这些变量之间的交换 alluvialClones (seurat_obj, cloneCall = "aa" , y.axes = c ( "sample" , "seurat_clusters" ), color = c ( "CAMGYCINNNCYEGWFDPW_CQQYYDTPRTF" , "CARALGAIELFDYW_CLLSYSGAHVVF" )) + scale_fill_manual ( values = c ( "grey" , colorblind_vector[ 3 ])) # Warning in colnames(meta) == color: longer object length is not a multiple of # shorter object length # Scale for fill is already present. # Adding another scale for fill, which will replace the existing scale. # 使用 circlize R 包中的和弦图来可视化簇之间的相互连接 circles <- getCirclize (seurat_obj, group.by = "seurat_clusters" ) # 给每个cluster赋予颜色 grid.cols <- hue_pal ( length ( unique (seurat_obj $ seurat_clusters))) names (grid.cols) <- unique (seurat_obj $ seurat_clusters) # 绘制弦图 chordDiagram (circles, self.link = 1 , grid.col = grid.cols) # 计算和绘制 STARTRAC 的三种参数 StartracDiversity (seurat_obj, type = "seurat_clusters" , group.by = "orig.ident" ) # clonalBias也可以用来量化单个克隆在不同细胞群或功能区间中的偏向性的新指标 # 衡量某个克隆是否倾向于集中在一个特定的细胞群（cluster）或功能区间（compartment） Idents (seurat_obj) <- "sample" clonalBias (seurat_obj, cloneCall = "aa" , split.by = "sample" , group.by = "seurat_clusters" , n.boots = 10 , min.expand = 1 ) # Smoothing formula not specified. Using: y ~ qss(x, lambda = 3) # Warning: Computation failed in `stat_quantile`. # Caused by error in `validObject`: # ! invalid class "matrix.csr" object: ra and ja don't have equal lengths 版本信息 大家如果遇到报错，请先检查一下各包的版本是否与我一致：sessionInfo # R version 4.4.2 (2024-10-31) # Platform: x86_64-pc-linux-gnu # Running under: Ubuntu 20.04.4 LTS # # Matrix products: default # BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 # LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3; LAPACK version 3.9.0 # # locale: # [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C # [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 # [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 # [7] LC_PAPER=en_US.UTF-8 LC_NAME=C # [9] LC_ADDRESS=C LC_TELEPHONE=C # [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C # # time zone: Etc/UTC # tzcode source: system (glibc) # # attached base packages: # [1] stats graphics grDevices utils datasets methods base # # other attached packages: # [1] dplyr_1.1.4 scales_1.3.0 circlize_0.4.16 ggraph_2.2.1 # [5] Trex_0.99.12 Seurat_5.2.1 SeuratObject_5.0.2 sp_2.2-0 # [9] scRepertoire_2.2.1 ggplot2_3.5.1 # # loaded via a namespace (and not attached): # [1] cubature_2.1.4 RcppAnnoy_0.0.22 # [3] splines_4.4.2 later_1.4.1 # [5] tibble_3.2.1 polyclip_1.10-7 # [7] fastDummies_1.7.5 lifecycle_1.0.4 # [9] hdf5r_1.3.12 globals_0.16.3 # [11] lattice_0.22-6 MASS_7.3-64 # [13] magrittr_2.0.3 plotly_4.10.4 # [15] sass_0.4.9 rmarkdown_2.29 # [17] jquerylib_0.1.4 yaml_2.3.10 # [19] httpuv_1.6.15 sctransform_0.4.1 # [21] spam_2.11-1 spatstat.sparse_3.1-0 # [23] reticulate_1.40.0 cowplot_1.1.3 # [25] pbapply_1.7-2 RColorBrewer_1.1-3 # [27] keras_2.15.0 abind_1.4-8 # [29] zlibbioc_1.52.0 Rtsne_0.17 # [31] GenomicRanges_1.58.0 purrr_1.0.2 # [33] BiocGenerics_0.52.0 hash_2.2.6.3 # [35] tweenr_2.0.3 evmix_2.12 # [37] GenomeInfoDbData_1.2.13 IRanges_2.40.1 # [39] S4Vectors_0.44.0 ggrepel_0.9.6 # [41] irlba_2.3.5.1 listenv_0.9.1 # [43] spatstat.utils_3.1-4 openintro_2.5.0 # [45] iNEXT_3.0.2 goftest_1.2-3 # [47] MatrixModels_0.5-3 airports_0.1.0 # [49] RSpectra_0.16-2 spatstat.random_3.4-1 # [51] fitdistrplus_1.2-2 parallelly_1.42.0 # [53] codetools_0.2-20 DelayedArray_0.32.0 # [55] ggforce_0.4.2 shape_1.4.6.1 # [57] tidyselect_1.2.1 UCSC.utils_1.2.0 # [59] farver_2.1.2 viridis_0.6.5 # [61] base64enc_0.1-3 matrixStats_1.5.0 # [63] stats4_4.4.2 spatstat.explore_3.4-3 # [65] jsonlite_1.8.9 tidygraph_1.3.1 # [67] progressr_0.15.1 ggridges_0.5.6 # [69] ggalluvial_0.12.5 survival_3.8-3 # [71] tools_4.4.2 stringdist_0.9.15 # [73] ica_1.0-3 Rcpp_1.0.14 # [75] glue_1.8.0 tfruns_1.5.3 # [77] gridExtra_2.3 SparseArray_1.6.0 # [79] xfun_0.50 MatrixGenerics_1.18.1 # [81] GenomeInfoDb_1.42.1 withr_3.0.2 # [83] fastmap_1.2.0 SparseM_1.84-2 # [85] digest_0.6.37 R6_2.5.1 # [87] mime_0.12 colorspace_2.1-1 # [89] scattermore_1.2 tensor_1.5 # [91] spatstat.data_3.1-6 tidyr_1.3.1 # [93] generics_0.1.3 data.table_1.16.4 # [95] usdata_0.3.1 graphlayouts_1.2.2 # [97] httr_1.4.7 htmlwidgets_1.6.4 # [99] S4Arrays_1.6.0 whisker_0.4.1 # [101] uwot_0.2.2 pkgconfig_2.0.3 # [103] gtable_0.3.6 tensorflow_2.16.0 # [105] lmtest_0.9-40 SingleCellExperiment_1.28.1 # [107] XVector_0.46.0 htmltools_0.5.8.1 # [109] dotCall64_1.2 Biobase_2.66.0 # [111] png_0.1-8 spatstat.univar_3.1-3 # [113] ggdendro_0.2.0 knitr_1.49 # [115] rstudioapi_0.17.1 tzdb_0.4.0 # [117] reshape2_1.4.4 rjson_0.2.23 # [119] nlme_3.1-168 GlobalOptions_0.1.2 # [121] cachem_1.1.0 zoo_1.8-12 # [123] stringr_1.5.1 KernSmooth_2.23-26 # [125] parallel_4.4.2 miniUI_0.1.1.1 # [127] pillar_1.10.1 grid_4.4.2 # [129] vctrs_0.6.5 RANN_2.6.2 # [131] VGAM_1.1-12 promises_1.3.2 # [133] xtable_1.8-4 cluster_2.1.8 # [135] evaluate_1.0.3 isoband_0.2.7 # [137] zeallot_0.1.0 readr_2.1.5 # [139] truncdist_1.0-2 cli_3.6.3 # [141] compiler_4.4.2 rlang_1.1.5 # [143] crayon_1.5.3 future.apply_1.11.3 # [145] labeling_0.4.3 plyr_1.8.9 # [147] stringi_1.8.4 deldir_2.0-4 # [149] viridisLite_0.4.2 assertthat_0.2.1 # [151] munsell_0.5.1 gsl_2.1-8 # [153] lazyeval_0.2.2 spatstat.geom_3.4-1 # [155] quantreg_5.99.1 Matrix_1.7-2 # [157] RcppHNSW_0.6.0 hms_1.1.3 # [159] patchwork_1.3.0 bit64_4.6.0-1 # [161] future_1.34.0 shiny_1.10.0 # [163] SummarizedExperiment_1.36.0 evd_2.3-7.1 # [165] ROCR_1.0-11 igraph_2.1.4 # [167] memoise_2.0.1 bslib_0.9.0 # [169] bit_4.5.0.1 cherryblossom_0.1.0
五、参考 A toolkit for single-cell immune receptor analysis scRepertoire 2: Enhanced and efficient toolkit for single-cell immune profiling scRepertoire: An R-based toolkit for single-cell immune receptor analysis