10X单细胞（10X空间转录组）BCR（TCR）数据分析之(11)dandelion

好了，我们今天要继续这个专题的分享，内容很多，好好深入学习的话，估计需要一辈子😄。

第一点，TCR（BCR的聚类方式）Clustering by V-gene, J-gene and junction length

也就是说，聚类的选择有三种，V,J,还有一起用来聚类，（为什么不能单独D也聚类呢？这个大家可以思考一下）。

模式呢，有四种

第一种，Amino acid model，这种模式就是之前分享TCRdist采用的模式，但是不是简简单单的看看是不是相同的氨基酸残基就可以了，而是要对20种氨基酸进行划类，同是碱性氨基酸的替换，那么判定两者之间的distance很小，这个需要很强的背景，要考虑氨基酸的理化性质、电荷性质等等，比较难，感兴趣大家可以研究一下。当然，距离的度量还是Hamming distance。

第二种，Hamming distance model，这种模式与上面的模式不同，这个模式采用的是nucleotide sequences，也是比较简单的比较方法，但是不如第一种来的更加准确，毕竟行驶生物学功能的是蛋白质啊~~

第三种，人和小鼠 1-mer 模型。hh_s1f 和 mk_rs5nf 距离模型是从 YVanderHeidenU+13 中的人类 5-mer 靶向模型和 CDiNiroVanderHeiden+16 中的小鼠 5-mer 靶向模型的平均和对称化得出的单核苷酸距离矩阵。 它们大致类似于转换/颠换模型。 (这种模型在单细胞角度分析非常少见)。

第四种 Human and mouse 5-mer models，hh_s5f 和 mk_rs5nf 距离模型分别基于 YVanderHeidenU+13 中的人类 5-mer 靶向模型和 CDiNiroVanderHeiden+16 中的小鼠 5-mer 靶向模型。 靶向矩阵 T 具有跨列的 5 聚体和作为行的 5 聚体中心碱基突变的核苷酸。 给定核苷酸 5-mer 对 T[i,j] 的值是该 5-mer 突变 mut(j) 的可能性与中心碱基突变为给定核苷酸 sub(j) 的可能性的乘积 \rightarrow i)。 这个概率矩阵通过以下步骤转换为距离矩阵 D：

由于距离矩阵 D 不是对称的，因此可以指定 --sym 参数来计算 D(j\rightarrow i) 和 D(i\rightarrow j) 的平均值 (avg) 或最小值 (min)。 由 D 定义的每个核苷酸差异的距离对连接处的所有 5 聚体进行求和，以产生两个连接序列之间的距离。

这四种距离模式要根据情况来进行采用，通常我们采用氨基酸模式就可以。

第二点，Reassigning heavy chain IG V gene alleles

B 细胞免疫球蛋白受体的高通量测序为适应性免疫提供了前所未有的见解。 分析这些数据的关键步骤涉及通过将包含每个免疫球蛋白序列的种系 V、D 和 J 基因片段等位基因与已知 V(D)J 等位基因的数据库进行匹配来分配。 然而，对于利用以前未检测到的等位基因的序列，该过程将失败，其在群体中的频率尚不清楚。 （未知基因的判断目前在10X种也没有解决）。

TIgGER 是一种计算方法，通过首先从 V(D)J 重排序列确定个体携带的完整基因片段集（包括新的等位基因），显着改善了 V(D)J 等位基因分配。 TIgGER 然后可以从这些序列中推断出受试者的基因型，并使用该基因型来纠正最初的 V(D)J 等位基因分配。

TIgGER 的应用继续在人类中发现高频率新等位基因，突出了对这种方法的迫切需要。 （然而，TIgGER 可以并且已经用于其他物种的数据。）

核心能力

• Detecting novel alleles
• Inferring a subject’s genotype
• Correcting preliminary allele calls

这部分我们了解即可，非常难。

下面我们用示例来解决一下BCR序列的聚类问题

在寻找克隆/克隆型的topic上，有很多方法可用于对 BCR 进行聚类，几乎都涉及一些基于序列相似性的度量。 BCR community还维护了许多非常完善的指导方针和标准。 例如，immcantation 使用许多基于模型的方法 Gupta2015 根据长度归一化的连接汉明距离的分布对克隆进行分组，而其他方法则使用整个 BCR V(D)J 序列来定义克隆。

1、加载模块

import os
import pandas as pd
import pandas as pd
import numpy as np
import scanpy as sc
import warnings
import functools
import seaborn as sns
import scipy.stats
import anndata
import dandelion as ddl
ddl.logging.print_header()
# change directory to somewhere more workable
os.chdir(os.path.expanduser('/Users/kt16/Downloads/dandelion_tutorial/'))
# I'm importing scanpy here to make use of its logging module.
import scanpy as sc
sc.settings.verbosity = 3
import warnings
warnings.filterwarnings('ignore')
sc.logging.print_header()

2、读取文件(10X的标准文件)

folder_location = 'sc5p_v2_hs_PBMC_10k'
# or file_location = 'sc5p_v2_hs_PBMC_10k/'
vdj = ddl.read_10x_vdj(folder_location, filename_prefix='sc5p_v2_hs_PBMC_10k_b_filtered', verbose = True)
vdj
Dandelion class object with n_obs = 994 and n_contigs = 2601
    data: 'cell_id', 'is_cell_10x', 'sequence_id', 'high_confidence_10x', 'sequence_length_10x', 'locus', 'v_call', 'd_call', 'j_call', 'c_call', 'complete_vdj', 'productive', 'junction_aa', 'junction', 'consensus_count', 'duplicate_count', 'clone_id', 'raw_consensus_id_10x', 'sequence'
    metadata: 'clone_id', 'clone_id_by_size', 'locus_VJ', 'locus_VDJ', 'productive_VJ', 'productive_VDJ', 'v_call_VJ', 'v_call_VDJ', 'j_call_VJ', 'j_call_VDJ', 'c_call_VJ', 'c_call_VDJ', 'duplicate_count_VJ', 'duplicate_count_VDJ', 'duplicate_count_VJ_1', 'duplicate_count_VJ_2', 'duplicate_count_VDJ_1', 'duplicate_count_VDJ_2', 'duplicate_count_VJ_3', 'duplicate_count_VDJ_3', 'duplicate_count_VDJ_4', 'duplicate_count_VDJ_5', 'duplicate_count_VJ_4', 'duplicate_count_VJ_5', 'junction_aa_VJ', 'junction_aa_VDJ', 'status', 'status_summary', 'productive', 'productive_summary', 'isotype', 'isotype_summary', 'vdj_status', 'vdj_status_summary', 'constant_status_summary'
    distance: None
    edges: None
    layout: None
    graph: None
# read in the airr_rearrangement.tsv file
file_location = 'sc5p_v2_hs_PBMC_10k/sc5p_v2_hs_PBMC_10k_b_airr_rearrangement.tsv'###这个文件需要下载
vdj = ddl.read_10x_airr(file_location)

3、读取转录组数据(scnapy的读取)

adata = sc.read_10x_h5('sc5p_v2_hs_PBMC_10k/filtered_feature_bc_matrix.h5', gex_only=True)
adata.obs['sample_id'] = 'sc5p_v2_hs_PBMC_10k'
adata.var_names_make_unique()
adata
AnnData object with n_obs × n_vars = 10553 × 36601
    obs: 'sample_id'
    var: 'gene_ids', 'feature_types', 'genome'

Run QC on the transcriptome data.

ddl.pp.recipe_scanpy_qc(adata)
adata
AnnData object with n_obs × n_vars = 10553 × 36601
    obs: 'sample_id', 'scrublet_score', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'is_doublet', 'filter_rna'
    var: 'gene_ids', 'feature_types', 'genome'

运行 bcr 数据的过滤。 请注意，我使用 Dandelion 对象作为输入而不是 Pandas 数据框（两种类型的输入都可以使用。事实上，.tsv 的文件路径也可以）。

# The function will return both objects.
vdj, adata = ddl.pp.filter_contigs(vdj, adata)

Check the output V(D)J table

The vdj table is returned as a Dandelion class object in the .data slot; if a file was provided for filter_bcr above, a new file will be created in the same folder with the filtered prefix. Note that this V(D)J table is indexed based on contigs (sequence_id).

vdj
Dandelion class object with n_obs = 344 and n_contigs = 687
    data: 'cell_id', 'sequence_id', 'sequence', 'sequence_aa', 'productive', 'rev_comp', 'v_call', 'v_cigar', 'd_call', 'd_cigar', 'j_call', 'j_cigar', 'c_call', 'c_cigar', 'sequence_alignment', 'germline_alignment', 'junction', 'junction_aa', 'junction_length', 'junction_aa_length', 'v_sequence_start', 'v_sequence_end', 'd_sequence_start', 'd_sequence_end', 'j_sequence_start', 'j_sequence_end', 'c_sequence_start', 'c_sequence_end', 'consensus_count', 'duplicate_count', 'is_cell', 'locus', 'umi_count'
    metadata: 'locus_VDJ', 'locus_VJ', 'productive_VDJ', 'productive_VJ', 'v_call_VDJ', 'v_call_VJ', 'j_call_VDJ', 'j_call_VJ', 'c_call_VDJ', 'c_call_VJ', 'duplicate_count_VDJ', 'duplicate_count_VJ', 'duplicate_count_VDJ_1', 'duplicate_count_VJ_1', 'junction_aa_VDJ', 'junction_aa_VJ', 'status', 'status_summary', 'productive', 'productive_summary', 'isotype', 'isotype_summary', 'vdj_status', 'vdj_status_summary', 'constant_status_summary'
    distance: None
    edges: None
    layout: None
    graph: None

Check the AnnData object as well

adata
AnnData object with n_obs × n_vars = 10553 × 36601
    obs: 'sample_id', 'scrublet_score', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'is_doublet', 'filter_rna', 'has_contig', 'filter_contig_quality', 'filter_contig_VDJ', 'filter_contig_VJ', 'contig_QC_pass', 'filter_contig'
    var: 'gene_ids', 'feature_types', 'genome'

看一看基础信息

pd.crosstab(adata.obs['has_contig'], adata.obs['filter_contig'])

现在实际过滤 AnnData 对象并运行标准工作流程，从过滤基因和规范化数据开始

因为“过滤后”的 AnnData 对象作为过滤后但未处理的对象返回，我们仍然需要规范化并在此处运行通常的过程。 以下只是一个标准的处理过程，就是简单的单细胞数据处理。

# filter genes
sc.pp.filter_genes(adata, min_cells=3)
# Normalize the counts
sc.pp.normalize_total(adata, target_sum=1e4)
# Logarithmize the data
sc.pp.log1p(adata)
# Stash the normalised counts
adata.raw = adata

识别高变基因并回归

sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
sc.pl.highly_variable_genes(adata)
adata = adata[:, adata.var.highly_variable]
sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])
sc.pp.scale(adata, max_value=10)

接下来一些常规的分析套路

sc.tl.pca(adata, svd_solver='arpack')
sc.pl.pca_variance_ratio(adata, log=True, n_pcs = 50)

# Computing the neighborhood graph
sc.pp.neighbors(adata)
# Embedding the neighborhood graph
sc.tl.umap(adata)
# Clustering the neighborhood graph
sc.tl.leiden(adata)

sc.pl.umap(adata, color=['IGHM', 'JCHAIN'])

接下来的重点就是我们的BCR,Finding clones

The following is dandelion’s implementation of a rather conventional method to define clones, tl.find_clones.

Clone definition is based on the following criterias:

(1) Identical IGH V-J gene usage.
(2) Identical CDR3 junctional sequence length.
(3) CDR3 Junctional sequences attains a minimum of % sequence similarity, based on hamming distance. The similarity cut-off is tunable (default is 85%).
(4) Light chain usage. If cells within clones use different light chains, the clone will be splitted following the same conditions for heavy chains in (1-3) as above.

The ‘clone_id’ name follows a {A}{B}{C}_{D} format and largely reflects the conditions above where:

{A} indicates if the contigs use the same IGH V/J genes.
{B} indicates if IGH junctional sequences are equal in length.
{C} indicates if clones are splitted based on junctional hamming distance threshold
{D} indicates light chain pairing.

如果在克隆中仅检测到一组轻链使用，则不会注释最后一个位置。

Running tl.find_clones

该函数将采用文件路径、pandas DataFrame（例如，如果已经使用 pandas 读取过滤后的文件）或Dandelion类对象。 计算连接汉明距离的默认模式是使用 CDR3 连接氨基酸序列，通过 key 选项指定（None defaults to junction_aa）。 可以将其切换为使用 CDR3 连接核苷酸序列（key = 'junction'，甚至完整的 V(D)J 氨基酸序列（key = 'sequence_alignment_aa"），只要列名存在于 .data slot中。

如果要使用等位基因来定义 V-J 基因用法，请指定：

by_alleles = True

尽管可能需要调整一些参数，但使用相同的设置也可能对 TCR 进行聚类。

ddl.tl.find_clones(vdj)
vdj

as per convention，这将返回一个名为“clone_id”的新列。 如果提供文件路径作为输入，它还会自动将文件保存到文件名的基本目录中。 否则，将返回Dandelion对象。

vdj.metadata

clone_id	clone_id_by_size	sample_id	locus_VDJ	locus_VJ	productive_VDJ	productive_VJ	v_call_genotyped_VDJ	v_call_genotyped_VJ	j_call_VDJ	...	junction_aa_VJ	status	status_summary	productive	productive_summary	isotype	isotype_summary	vdj_status	vdj_status_summary	constant_status_summary
sc5p_v2_hs_PBMC_10k_AAACCTGTCCGTTGTC	148_2_2_202	1072	sc5p_v2_hs_PBMC_10k	IGH	IGK	T	T	IGHV1-69	IGKV1-8	IGHJ3	...	CQQYYSYPRTF	IGH + IGK	IGH + IGK	T + T	T + T	IgM	IgM	Single + Single	Single	Single

（可选）Running tl.define_clones（这个在之前分享过，文章在10X单细胞（10X空间转录组）TCR数据分析之TCRdist3(5)）

函数 pp.calculate_threshold 将运行 shazam 的 distToNearest 函数并返回显示长度归一化汉明距离分布和自动阈值的图。

同样， pp.calculate_threshold 将采用文件路径、pandas DataFrame 或 Dandelion 对象作为输入。 如果提供了Dandelion 对象，阈值将被插入到 .threshold slot中。 更精细的控制请直接使用shazam的distToNearest和changeo的DefineClones.py函数。

ddl.pp.calculate_threshold(vdj)

当然，阈值也可以人工选择，之前已经分享过了。类似于下面的代码

ddl.pp.calculate_threshold(vdj, manual_threshold = 0.1)

Generation of V(D)J network

dandelion 生成一个网络以促进结果的可视化，其灵感来自 Bashford-Rogers13。 这使用完整的 V(D)J contig 序列而不是仅仅连接序列来绘制每个克隆的树状网络。 后面会通过scanpy实现实际的可视化。

tl.generate_network

首先我们需要生成网络。 tl.generate_network 将采用已定义克隆的 V(D)J 表，特别是在 'clone_id' 列下。 默认模式是使用氨基酸序列来构建 Levenshtein 距离矩阵，但可以使用 key 选项进行切换。

如果从 immcantation 的方法或任何其他方法解析的预处理表，只要它是 AIRR 格式，也可以使用该表。

可以指定 clone_key 选项来为选择的克隆 ID 定义生成网络，只要它作为 .data slot中的列存在。

ddl.tl.generate_network(vdj)

This step works reasonably fast here but will take quite a while when a lot of contigs are provided.

这篇的最后我们来看看可视化

tl.transfer

ddl.tl.transfer(adata, vdj) # this will include singletons. To show only expanded clones, specify expanded_only=True
adata

import matplotlib as mpl
mpl.rcParams.update(mpl.rcParamsDefault)
ddl.pl.barplot(vdj,
               color = 'v_call_genotyped_VDJ',
               figsize = (12, 4))

ddl.pl.barplot(vdj,
               color = 'v_call_genotyped_VDJ',
               figsize = (12, 4),
               sort_descending = None,
               palette = 'tab20')

ddl.pl.barplot(vdj,
               color = 'v_call_genotyped_VDJ',
               normalize = False,
               figsize = (12, 4),
               sort_descending = None,
               palette = 'tab20')

import matplotlib.pyplot as plt
ddl.pl.stackedbarplot(vdj,
                      color = 'isotype',
                      groupby = 'status',
                      xtick_rotation =0,
                      figsize = (4,4))
plt.legend(bbox_to_anchor = (1,1),
           loc='upper left',
           frameon=False)

ddl.pl.stackedbarplot(vdj,
                      color = 'v_call_genotyped_VDJ',
                      groupby = 'isotype')
plt.legend(bbox_to_anchor = (1,1),
           loc='upper left',
           frameon=False)

ddl.pl.stackedbarplot(vdj,
                      color = 'v_call_genotyped_VDJ',
                      groupby = 'isotype',
                      normalize = True)
plt.legend(bbox_to_anchor = (1,1),
           loc='upper left',
           frameon=False)

图片.png

ddl.pl.stackedbarplot(vdj,
                      color = 'v_call_genotyped_VDJ',
                      groupby = 'vdj_status')
plt.legend(bbox_to_anchor = (1,1),
           loc='upper left',
           frameon=False)

ddl.pl.stackedbarplot(vdj,
                      color = 'v_call_genotyped_VDJ',
                      groupby = 'sample_id')
plt.legend(bbox_to_anchor = (1, 0.5),
           loc='center left',
           frameon=False)

ddl.pl.spectratype(vdj,
                   color = 'junction_length',
                   groupby = 'c_call',
                   locus='IGH',
                   width = 2.3)
plt.legend(bbox_to_anchor = (1,1),
           loc='upper left',
           frameon=False)

ddl.pl.spectratype(vdj,
                   color = 'junction_aa_length',
                   groupby = 'c_call',
                   locus='IGH')
plt.legend(bbox_to_anchor = (1,1),
           loc='upper left',
           frameon=False)

图片.png

ddl.pl.spectratype(vdj,
                   color = 'junction_aa_length',
                   groupby = 'c_call',
                   locus=['IGK','IGL'])
plt.legend(bbox_to_anchor = (1,1),
           loc='upper left',
           frameon=False)

图片.png

pl.clone_overlap(重点来了)

ddl.tl.clone_overlap(adata,
                     groupby = 'leiden',
                     colorby = 'leiden')
ddl.pl.clone_overlap(adata,
                     groupby = 'leiden',
                     colorby = 'leiden',
                     return_graph=True,
                     group_label_offset=.5)

图片.png

今天我们就分享到这里吧，问题非常多，但是欲速则不达

生活很好，有你更好