免疫组库基础分析15：利用 NAIR 包探索疾病表型或治疗相关的 TCR/BCR 簇

三兔测序学社

发布于 2026-04-17 16:49:50

1220

编者语：

这一篇介绍了NAIR包最重要的功能，就是寻找找到与疾病状态、治疗响应或临床结局显著相关的一组相似免疫受体序列。这也是编者自己主要用的NAIR包核心代码，就是为了这一点醋，才写了前面三篇NAIR包的介绍博文。当然，这样的编写是NAIR包说明书本身的顺序，但实际生物信息学学习过程，其实可以直接从这一篇开始应用即可，再回头去学习代码基础参数与可选参数功能等内容，工具用过完成一个任务即可，不必了解如何制作过程。当然，由于每一个人的分析目的或者默认分析结果不理想，则需要个性化的研究，还是得老老实实阅读英文说明书。

在免疫组库测序（AIRR-Seq）研究中，我们常常不只关心单个 TCR/BCR 序列，更希望找到与疾病状态、治疗响应或临床结局显著相关的一组相似免疫受体序列，也就是免疫受体簇。NAIR 包提供了一套完整流程，能从批量（bulk)免疫组库样本数据中，系统性挖掘这类具有生物学意义的关联簇。

本文将完整介绍这一分析流程，从关联序列筛选、相似克隆识别，到全局网络构建与聚类分析。

一、引言

当我们拥有多份批量适应性免疫受体库测序（bulk -AIRR-Seq）样本时，可以借助 NAIR 包，研究与二分类变量相关的 TCR/BCR 簇。这里的二分类变量可以是疾病 / 健康、治疗组 / 对照组、预后良好 / 不良等。

整个分析的核心思路是：先找到组间差异显著的序列，再扩展到与之相似的序列克隆，最后构建网络并划分簇，得到最终的关联簇。

二、分析流程概览

1.识别关联序列

根据二分类变量将研究对象分为两组，使用费希尔精确检验，筛选出在两组间频率存在统计学显著差异的 TCR/BCR 序列。

2.识别与关联序列相似的克隆

对每个关联序列，设定序列相似度阈值（例如氨基酸差异不超过 1 个），定义其 “邻域序列”，并从所有样本中找出属于这些邻域的克隆。

3.构建关联克隆的全局网络

将上一步筛选出的所有克隆整合为一个全局网络，通过聚类分析将网络划分为不同簇，这些簇即为与目标表型相关的关联簇。

4.后续分析与可视化

在网络图中标记聚类编号、单独分析感兴趣的特定簇、绘制不同分组的着色方案等。

三、安装并加载 NAIR 包

首先安装并加载分析所需的核心包，这是所有步骤的基础。

# 安装NAIR包（首次使用执行）
install.packages("NAIR")
# 加载包
library(NAIR)

四、模拟演示数据

为了清晰展示整个流程，我们先构建一套模拟数据，完全贴合真实研究的数据结构。

set.seed(42)
library(NAIR)
data_dir <- getwd()
dir_input_samples <- file.path(data_dir, "input_samples")
dir.create(dir_input_samples, showWarnings = FALSE)
# Number of samples by control/treatment group
n_control <- n_treatment <- 15 
n_samples <- n_control + n_treatment
sample_size <- 30 # (seqs per sample)          
# sequences (first five are chosen to be associated with treatment)
base_seqs <- 
  c("CASSGAYEQYF", "CSVDLGKGNNEQFF", 
    "CASSIEGQLSTDTQYF", # 3 of the associated
    "CASSEEGQLSTDTQYF", # sequences differ by
    "CASSPEGQLSTDTQYF", # only one amino acid
    "RASSLAGNTEAFF", "CASSHRGTDTQYF", "CASDAGVFQPQHF"
  ) 
# relative generation probabilities by control/treatment group
pgen_c <- matrix(rep(c(rep(1, 5), rep(30, 3)), times = n_control), 
                 nrow = n_control, byrow = TRUE)
pgen_t <- matrix(rep(c(1, 1, rep(1/3, 3), rep(2, 3)), times = n_treatment), 
                 nrow = n_treatment, byrow = TRUE)
pgen <- rbind(pgen_c, pgen_t)
simulateToyData(    
  samples = n_samples, 
  sample_size = sample_size,
  prefix_length = 1, 
  prefix_chars = c("", ""),
  prefix_probs = cbind(rep(1, n_samples), rep(0, n_samples)),
  affixes = base_seqs, 
  affix_probs = pgen, 
  num_edits = 0,
  output_dir = dir_input_samples, 
  no_return = TRUE
)
# 输入文件与分组标签
group_labels <- c(rep("control", n_control), rep("treatment", n_treatment))
input_files <- file.path(dir_input_samples,paste0("Sample", 1:n_samples, ".rds"))
length(input_files)

五、步骤 1：寻找关联序列

第一步核心：筛选两组间分布差异显著的受体序列，使用findAssociatedSeqs()函数完成。根据二分类变量将研究对象分为两组，使用费希尔精确检验，筛选出在两组间频率存在统计学显著差异的 TCR/BCR 序列。

associated_seqs <- findAssociatedSeqs(file_list = input_files, 
#每一个样本的路径信息
                                      input_type = "rds", 
                                      group_ids = group_labels, #分组信息
                                      seq_col = "CloneSeq", 
                                      min_seq_length = NULL, 
                                      drop_matches = NULL, 
                                      min_sample_membership = NULL, 
                                      pval_cutoff = 0.1  #p值小于0.1的克隆进行下游分析
)

结果如下：在后续的处理步骤中，只需要关注并使用 ReceptorSeq 这一列的数据

associated_seqs[, 1:5]
#>        ReceptorSeq fisher_pvalue shared_by_n_samples samples_g0 samples_g1
#> 8   CSVDLGKGNNEQFF  1.052106e-05                  18          3         15
#> 7      CASSGAYEQYF  1.157316e-04                  17          3         14
#> 4 CASSEEGQLSTDTQYF  5.197401e-03                  10          1          9
#> 5 CASSIEGQLSTDTQYF  6.559548e-02                  16          5         11

六、步骤 2：查找关联克隆

基于关联序列，找到所有相似的克隆，构建序列邻域数据集

# 建立一个结果文件夹
dir_nbds <- file.path(data_dir, "assoc_seq_nbds")
# Identify clones in a neighborhood around each associated sequence
findAssociatedClones(file_list = input_files, 
                     input_type = "rds", 
                     group_ids = group_labels, 
                     seq_col = "CloneSeq", 
                     assoc_seqs = associated_seqs$ReceptorSeq,#第一步中组间有差异的TCR、BCR序列
                     min_seq_length = NULL, 
                     drop_matches = NULL,
                     output_dir = dir_nbds  #输入文件夹
)

结果如下

# first neighborhood data file
one_nbd <- readRDS(list.files(dir_nbds, full.names = TRUE)[[1]])
# first few rows
head(one_nbd)
#>                   CloneSeq CloneFrequency CloneCount SampleID    GroupID
#> Sample5.14  CSVDLGKGNNEQFF     0.02668662       2924  Sample5  reference
#> Sample7.17  CSVDLGKGNNEQFF     0.01957279       2113  Sample7  reference
#> Sample10.17 CSVDLGKGNNEQFF     0.03594648       4097 Sample10  reference
#> Sample16.10 CSVDLGKGNNEQFF     0.03404553       3802 Sample16 comparison
#> Sample16.28 CSVDLGKGNNEQFF     0.02736537       3056 Sample16 comparison
#> Sample17.3  CSVDLGKGNNEQFF     0.04114818       4603 Sample17 comparison
#>                   AssocSeq
#> Sample5.14  CSVDLGKGNNEQFF
#> Sample7.17  CSVDLGKGNNEQFF
#> Sample10.17 CSVDLGKGNNEQFF
#> Sample16.10 CSVDLGKGNNEQFF
#> Sample16.28 CSVDLGKGNNEQFF
#> Sample17.3  CSVDLGKGNNEQFF

七、步骤 3：构建关联簇的全局网络

整合所有相似克隆，构建全局网络并划分关联簇，这是核心分析步骤。

#创建步骤2输出文件路径向量
nbd_files <- list.files(dir_nbds, full.names = TRUE)
# Combine neighborhoods and perform network analysis
all_clusters <- 
  buildAssociatedClusterNetwork(
    file_list = nbd_files, 
#输入文件路径
    seq_col = "CloneSeq",
    size_nodes_by = 1.5,
    print_plots = TRUE
  )
# 查看网络核心结果
names(all_clusters )
head(all_clusters$node_data) # 节点信息
head(all_clusters$cluster_data) # 聚类信息

八、步骤 4：附加分析与可视化

对网络结果进行精细化分析，提升结果的实用性和可读性。

# 全局网络图，节点按关联序列着色
all_clusters <- addPlots(all_clusters,
                         color_nodes_by = "AssocSeq",
                         color_title = "Neighborhood Sequence",
                         size_nodes_by = 1.5,
                         print_plots = TRUE
)
#标记cluster簇编号
all_clusters <- labelClusters(all_clusters, size = 10)
all_clusters$plots[[1]]
# 关注特定的克隆簇
buildNet(
data = all_clusters$node_data[all_clusters$node_data$cluster_id == 2, ],
#选择簇2
seq_col = "CloneSeq", 
color_nodes_by = c("CloneSeq", "SampleID"), 
#选择序列或者样本ID进行分组可视化
color_scheme = c("plasma", "turbo"),
size_nodes_by = 3, 
plot_title = "Cluster 2",
print_plots = TRUE
)

九、结果保存

将所有分析结果保存，方便后续论文写作与数据复用。

# 保存网络结果
saveRDS(all_clusters, "final_cluster_network.rds")
# 保存聚类表格
write.csv(all_clusters$cluster_data, "cluster_statistics.csv", row.names = FALSE)
write.csv(all_clusters$node_data, "node_information.csv", row.names = FALSE)

注意事项：

步骤	参数名称	作用	建议关注点
Step 1	group_ids	定义样本分组	必须准确，错一位结果全错。
Step 1	pval_cutoff	设定显著性门槛	根据数据量调整，多重检验时可能需要更严格。
Step 2	nbd_radius	定义“相似”范围	设为 0 只找完全相同的序列；设为 1 找突变体。
Step 3	cluster_fun	选择聚类算法	默认算法通常可用，若聚类效果不好可尝试更换。
Step 3	color_nodes_by	可视化着色	调试时可尝试不同值（如 SampleID）以排查数据问题。