好了砸逊，我們今天要繼續(xù)這個專題的分享赤屋，內(nèi)容很多壁袄，好好深入學(xué)習(xí)的話，估計(jì)需要一輩子??嗜逻。

第一點(diǎn)，TCR（BCR的聚類方式）Clustering by V-gene, J-gene and junction length

也就是說逆日，聚類的選擇有三種萄凤，V,J,還有一起用來聚類，（為什么不能單獨(dú)D也聚類呢靡努？這個大家可以思考一下）坪圾。

模式呢，有四種

第一種惑朦，Amino acid model兽泄，這種模式就是之前分享TCRdist采用的模式，但是不是簡簡單單的看看是不是相同的氨基酸殘基就可以了漾月，而是要對20種氨基酸進(jìn)行劃類病梢，同是堿性氨基酸的替換，那么判定兩者之間的distance很小梁肿，這個需要很強(qiáng)的背景蜓陌，要考慮氨基酸的理化性質(zhì)、電荷性質(zhì)等等栈雳，比較難护奈，感興趣大家可以研究一下。當(dāng)然哥纫，距離的度量還是Hamming distance霉旗。

第二種，Hamming distance model蛀骇，這種模式與上面的模式不同厌秒，這個模式采用的是nucleotide sequences，也是比較簡單的比較方法擅憔，但是不如第一種來的更加準(zhǔn)確鸵闪，畢竟行駛生物學(xué)功能的是蛋白質(zhì)啊~~

第三種，人和小鼠 1-mer 模型暑诸。hh_s1f 和 mk_rs5nf 距離模型是從 YVanderHeidenU+13 中的人類 5-mer 靶向模型和 CDiNiroVanderHeiden+16 中的小鼠 5-mer 靶向模型的平均和對稱化得出的單核苷酸距離矩陣蚌讼。 它們大致類似于轉(zhuǎn)換/顛換模型辟灰。 (這種模型在單細(xì)胞角度分析非常少見)。

圖片.png

第四種 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) 的可能性與中心堿基突變?yōu)榻o定核苷酸 sub(j) 的可能性的乘積 \rightarrow i)。 這個概率矩陣通過以下步驟轉(zhuǎn)換為距離矩陣 D：

圖片.png

由于距離矩陣 D 不是對稱的朗和，因此可以指定 --sym 參數(shù)來計(jì)算 D(j\rightarrow i) 和 D(i\rightarrow j) 的平均值 (avg) 或最小值 (min)眶拉。 由 D 定義的每個核苷酸差異的距離對連接處的所有 5 聚體進(jìn)行求和，以產(chǎn)生兩個連接序列之間的距離谒臼。

這四種距離模式要根據(jù)情況來進(jìn)行采用拾氓，通常我們采用氨基酸模式就可以咙鞍。

第二點(diǎn)续滋，Reassigning heavy chain IG V gene alleles

B 細(xì)胞免疫球蛋白受體的高通量測序?yàn)檫m應(yīng)性免疫提供了前所未有的見解疲酌。 分析這些數(shù)據(jù)的關(guān)鍵步驟涉及通過將包含每個免疫球蛋白序列的種系 V湿颅、D 和 J 基因片段等位基因與已知 V(D)J 等位基因的數(shù)據(jù)庫進(jìn)行匹配來分配肖爵。 然而，對于利用以前未檢測到的等位基因的序列揉稚，該過程將失敗，其在群體中的頻率尚不清楚驻呐。 （未知基因的判斷目前在10X種也沒有解決）猜拾。

TIgGER 是一種計(jì)算方法挎袜，通過首先從 V(D)J 重排序列確定個體攜帶的完整基因片段集（包括新的等位基因）盯仪，顯著改善了 V(D)J 等位基因分配全景。 TIgGER 然后可以從這些序列中推斷出受試者的基因型，并使用該基因型來糾正最初的 V(D)J 等位基因分配奔浅。

TIgGER 的應(yīng)用繼續(xù)在人類中發(fā)現(xiàn)高頻率新等位基因鲁驶，突出了對這種方法的迫切需要钥弯。 （然而脆霎，TIgGER 可以并且已經(jīng)用于其他物種的數(shù)據(jù)睛蛛。）

核心能力

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

這部分我們了解即可荸频，非常難旭从。

下面我們用示例來解決一下BCR序列的聚類問題

在尋找克隆/克隆型的topic上和悦，有很多方法可用于對 BCR 進(jìn)行聚類鸽素，幾乎都涉及一些基于序列相似性的度量。 BCR community還維護(hù)了許多非常完善的指導(dǎo)方針和標(biāo)準(zhǔn)蚜迅。 例如坐梯，immcantation 使用許多基于模型的方法 Gupta2015 根據(jù)長度歸一化的連接漢明距離的分布對克隆進(jìn)行分組吵血，而其他方法則使用整個 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的標(biāo)準(zhǔn)文件)

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支竹、讀取轉(zhuǎn)錄組數(shù)據(jù)(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'

運(yùn)行 bcr 數(shù)據(jù)的過濾柳洋。 請注意，我使用 Dandelion 對象作為輸入而不是 Pandas 數(shù)據(jù)框（兩種類型的輸入都可以使用募书。事實(shí)上莹捡，.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'

看一看基礎(chǔ)信息

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

圖片.png

現(xiàn)在實(shí)際過濾 AnnData 對象并運(yùn)行標(biāo)準(zhǔn)工作流程，從過濾基因和規(guī)范化數(shù)據(jù)開始

因?yàn)椤斑^濾后”的 AnnData 對象作為過濾后但未處理的對象返回寥茫，我們?nèi)匀恍枰?guī)范化并在此處運(yùn)行通常的過程纱耻。 以下只是一個標(biāo)準(zhǔn)的處理過程弄喘，就是簡單的單細(xì)胞數(shù)據(jù)處理蘑志。

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

接下來一些常規(guī)的分析套路

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

圖片.png

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

圖片.png

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

圖片.png

接下來的重點(diǎn)就是我們的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.

如果在克隆中僅檢測到一組輕鏈?zhǔn)褂醚蚴迹瑒t不會注釋最后一個位置突委。

Running tl.find_clones

該函數(shù)將采用文件路徑匀油、pandas DataFrame（例如桥滨，如果已經(jīng)使用 pandas 讀取過濾后的文件）或Dandelion類對象弛车。 計(jì)算連接漢明距離的默認(rèn)模式是使用 CDR3 連接氨基酸序列喻括，通過 key 選項(xiàng)指定（None defaults to junction_aa）唬血。 可以將其切換為使用 CDR3 連接核苷酸序列（key = 'junction'拷恨，甚至完整的 V(D)J 氨基酸序列（key = 'sequence_alignment_aa"）挑随，只要列名存在于 .data slot中兜挨。

如果要使用等位基因來定義 V-J 基因用法拌汇，請指定：

by_alleles = True

盡管可能需要調(diào)整一些參數(shù)噪舀，但使用相同的設(shè)置也可能對 TCR 進(jìn)行聚類。

ddl.tl.find_clones(vdj)
vdj

as per convention界逛，這將返回一個名為“clone_id”的新列。 如果提供文件路徑作為輸入净响，它還會自動將文件保存到文件名的基本目錄中。 否則赞别，將返回Dandelion對象仿滔。

vdj.metadata

（可選）Running tl.define_clones（這個在之前分享過，文章在10X單細(xì)胞（10X空間轉(zhuǎn)錄組）TCR數(shù)據(jù)分析之TCRdist3(5)）

函數(shù) pp.calculate_threshold 將運(yùn)行 shazam 的 distToNearest 函數(shù)并返回顯示長度歸一化漢明距離分布和自動閾值的圖。

同樣盐固， pp.calculate_threshold 將采用文件路徑刁卜、pandas DataFrame 或 Dandelion 對象作為輸入。 如果提供了Dandelion 對象，閾值將被插入到 .threshold slot中。 更精細(xì)的控制請直接使用shazam的distToNearest和changeo的DefineClones.py函數(shù)羔挡。

ddl.pp.calculate_threshold(vdj)

圖片.png

當(dāng)然低矮，閾值也可以人工選擇，之前已經(jīng)分享過了肠虽。類似于下面的代碼

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

圖片.png

Generation of V(D)J network

dandelion 生成一個網(wǎng)絡(luò)以促進(jìn)結(jié)果的可視化，其靈感來自 Bashford-Rogers13找颓。 這使用完整的 V(D)J contig 序列而不是僅僅連接序列來繪制每個克隆的樹狀網(wǎng)絡(luò)益老。 后面會通過scanpy實(shí)現(xiàn)實(shí)際的可視化。

tl.generate_network

首先我們需要生成網(wǎng)絡(luò)。 tl.generate_network 將采用已定義克隆的 V(D)J 表呛牲，特別是在 'clone_id' 列下着茸。 默認(rèn)模式是使用氨基酸序列來構(gòu)建 Levenshtein 距離矩陣敬特，但可以使用 key 選項(xiàng)進(jìn)行切換掰伸。

如果從 immcantation 的方法或任何其他方法解析的預(yù)處理表合搅，只要它是 AIRR 格式赌髓，也可以使用該表春弥。

可以指定 clone_key 選項(xiàng)來為選擇的克隆 ID 定義生成網(wǎng)絡(luò)，只要它作為 .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))

圖片.png

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

圖片.png

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

圖片.png

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)

圖片.png

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

圖片.png

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)

圖片.png

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)

圖片.png

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)

圖片.png

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(重點(diǎn)來了)

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

今天我們就分享到這里吧忍燥，問題非常多梅垄，但是欲速則不達(dá)

生活很好，有你更好