cellranger VDJ 数据过滤

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

df = pd.read_csv("/mnt/g/20220309-scBCR/HY01-1F11_ALL.csv",sep=",",low_memory=False)
df

	n	datasets	origins	donors	entropy_cell	near	far	dref	dref_aa	ext	...	group_id	group_ncells	clonotype_id	clonotype_ncells	nchains	exact_subclonotype_id	barcodes	HY01-1F11_barcodes	barcode	HY01-1F11_barcode
0	59	HY01-1F11	s1	d1	0.0	NaN	NaN	14	6	NaN	...	1	1477	1	1477	3	1	AAAGCAAAGTACGCGA-1,AACACGTGTTGGTTTG-1,ACACCGGG...	AAAGCAAAGTACGCGA-1,AACACGTGTTGGTTTG-1,ACACCGGG...	AAAGCAAAGTACGCGA-1	AAAGCAAAGTACGCGA-1
1	59	HY01-1F11	s1	d1	0.0	NaN	NaN	14	6	NaN	...	1	1477	1	1477	3	1	AAAGCAAAGTACGCGA-1,AACACGTGTTGGTTTG-1,ACACCGGG...	AAAGCAAAGTACGCGA-1,AACACGTGTTGGTTTG-1,ACACCGGG...	AACACGTGTTGGTTTG-1	AACACGTGTTGGTTTG-1
2	59	HY01-1F11	s1	d1	0.0	NaN	NaN	14	6	NaN	...	1	1477	1	1477	3	1	AAAGCAAAGTACGCGA-1,AACACGTGTTGGTTTG-1,ACACCGGG...	AAAGCAAAGTACGCGA-1,AACACGTGTTGGTTTG-1,ACACCGGG...	ACACCGGGTACAGTGG-1	ACACCGGGTACAGTGG-1
3	59	HY01-1F11	s1	d1	0.0	NaN	NaN	14	6	NaN	...	1	1477	1	1477	3	1	AAAGCAAAGTACGCGA-1,AACACGTGTTGGTTTG-1,ACACCGGG...	AAAGCAAAGTACGCGA-1,AACACGTGTTGGTTTG-1,ACACCGGG...	ACAGCCGAGAATCTCC-1	ACAGCCGAGAATCTCC-1
4	59	HY01-1F11	s1	d1	0.0	NaN	NaN	14	6	NaN	...	1	1477	1	1477	3	1	AAAGCAAAGTACGCGA-1,AACACGTGTTGGTTTG-1,ACACCGGG...	AAAGCAAAGTACGCGA-1,AACACGTGTTGGTTTG-1,ACACCGGG...	ACCCACTTCGCCGTGA-1	ACCCACTTCGCCGTGA-1
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
12016	1	HY01-1F11	s1	d1	0.0	NaN	NaN	0	0	NaN	...	3300	1	1	1	2	1	ACGCAGCTCCAAGCCG-1	ACGCAGCTCCAAGCCG-1	ACGCAGCTCCAAGCCG-1	ACGCAGCTCCAAGCCG-1
12017	1	HY01-1F11	s1	d1	0.0	NaN	NaN	2	0	NaN	...	3301	1	1	1	2	1	GCGCGATTCTGCGGCA-1	GCGCGATTCTGCGGCA-1	GCGCGATTCTGCGGCA-1	GCGCGATTCTGCGGCA-1
12018	1	HY01-1F11	s1	d1	0.0	NaN	NaN	7	6	NaN	...	3302	1	1	1	2	1	AGAATAGTCGCTGATA-1	AGAATAGTCGCTGATA-1	AGAATAGTCGCTGATA-1	AGAATAGTCGCTGATA-1
12019	1	HY01-1F11	s1	d1	0.0	NaN	NaN	4	4	NaN	...	3303	1	1	1	3	1	AGTGTCATCACTTATC-1	AGTGTCATCACTTATC-1	AGTGTCATCACTTATC-1	AGTGTCATCACTTATC-1
12020	1	HY01-1F11	s1	d1	0.0	NaN	NaN	8	6	NaN	...	3304	1	1	1	2	1	GTGTTAGTCTAACTGG-1	GTGTTAGTCTAACTGG-1	GTGTTAGTCTAACTGG-1	GTGTTAGTCTAACTGG-1

12021 rows × 380 columns

df1 = df.drop_duplicates(subset=["group_id","exact_subclonotype_id"],keep="first")

df2 = df1[(df1.nchains >= 2) & (df1.group_ncells >= 2)]
df3 = df2.loc[:,['group_id', 'group_ncells', 'exact_subclonotype_id', 'n', 'barcodes', 'nchains', 'u1', 'r1', 'v_name1', 'd_name1', 'j_name1','const1', 'cdr3_aa1', 'cdr3_dna1', 'fwr1_aa1', 'fwr1_dna1', 'cdr1_aa1', 'cdr1_dna1', 'fwr2_aa1', 'fwr2_dna1', 'cdr2_aa1', 'cdr2_dna1', 'fwr3_aa1', 'fwr3_dna1', 'fwr4_aa1', 'fwr4_dna1', 'vj_aa_nl1', 'vj_seq_nl1', 'u2', 'r2', 'v_name2', 'd_name2', 'j_name2','const2', 'cdr3_aa2', 'cdr3_dna2', 'fwr1_aa2', 'fwr1_dna2', 'cdr1_aa2', 'cdr1_dna2', 'fwr2_aa2', 'fwr2_dna2', 'cdr2_aa2', 'cdr2_dna2', 'fwr3_aa2', 'fwr3_dna2', 'fwr4_aa2', 'fwr4_dna2', 'vj_aa_nl2', 'vj_seq_nl2', 'u3', 'r3', 'v_name3', 'd_name3', 'j_name3','const3', 'cdr3_aa3', 'cdr3_dna3', 'fwr1_aa3', 'fwr1_dna3', 'cdr1_aa3', 'cdr1_dna3', 'fwr2_aa3', 'fwr2_dna3', 'cdr2_aa3', 'cdr2_dna3', 'fwr3_aa3', 'fwr3_dna3', 'fwr4_aa3', 'fwr4_dna3', 'vj_aa_nl3', 'vj_seq_nl3', 'u4', 'r4', 'v_name4', 'd_name4', 'j_name4','const4', 'cdr3_aa4', 'cdr3_dna4', 'fwr1_aa4', 'fwr1_dna4', 'cdr1_aa4', 'cdr1_dna4', 'fwr2_aa4', 'fwr2_dna4', 'cdr2_aa4', 'cdr2_dna4', 'fwr3_aa4', 'fwr3_dna4', 'fwr4_aa4', 'fwr4_dna4', 'vj_aa_nl4', 'vj_seq_nl4']]

df3.to_csv("/mnt/g/20220309-scBCR/HY01-1F11_filter.csv",index=False)

cell count

print("过滤后共有 %d 个细胞被保留下来进行进一步分析"%np.sum(df3.n))

过滤后共有 9256 个细胞被保留下来进行进一步分析

clonotype count

print("过滤后共有 %d 个clonotye被保留下来进行进一步分析"%len(set(df3.group_id)))

过滤后共有 539 个clonotye被保留下来进行进一步分析

plt.style.use('ggplot')
fig, ax = plt.subplots()

ax.bar(df3[df3.group_id <=10].group_id, df3[df3.group_id <= 10].group_ncells)
ax.set_ylabel('num_cells')
ax.set_title('Top 10 Clonotype')
#ax.set_xticks(ind, labels=['G1', 'G2', 'G3', 'G4', 'G5'])
#ax.legend()
plt.show()

CDR3 lengths distribute

new_cols = {x: y for x, y in zip(df3.loc[:,["const2","cdr3_aa2","v_name2","j_name2","n"]].columns,df3.loc[:,["const1","cdr3_aa1","v_name1","j_name1","n"]].columns)}
df_out = pd.concat([df3.loc[:,["const1","cdr3_aa1","v_name1","j_name1","n"]],df3.loc[:,["const2","cdr3_aa2","v_name2","j_name2","n"]].rename(columns=new_cols)],ignore_index=True)

new_cols = {x: y for x, y in zip(df3.loc[:,["const3","cdr3_aa3","v_name3","j_name3","n"]].columns,df3.loc[:,["const1","cdr3_aa1","v_name1","j_name1","n"]].columns)}
df_out = pd.concat([df_out,df3.loc[:,["const3","cdr3_aa3","v_name3","j_name3","n"]].rename(columns=new_cols)],ignore_index=True)

new_cols = {x: y for x, y in zip(df3.loc[:,["const4","cdr3_aa4","v_name4","j_name4","n"]].columns,df3.loc[:,["const1","cdr3_aa1","v_name1","j_name1","n"]].columns)}
df_out = pd.concat([df_out,df3.loc[:,["const4","cdr3_aa4","v_name4","j_name4","n"]].rename(columns=new_cols)],ignore_index=True)

df_out = df_out.dropna(how="any")
df_out

	const1	cdr3_aa1	v_name1	j_name1	n
0	IGHG1	CAPIHYDYGTWFAYW	IGHV14-3	IGHJ3	59
1	IGHG1	CAPISYDYGTWFAYW	IGHV14-3	IGHJ3	201
2	IGHG1	CAPIHYDYGTWFAYW	IGHV14-3	IGHJ3	174
3	IGHG1	CAPIYYDYGTWFAYW	IGHV14-3	IGHJ3	173
4	IGHG1	CAPIHYDYGTWFAYW	IGHV14-3	IGHJ3	92
...	...	...	...	...	...
5999	IGKC	CQQYWSTPYTF	IGKV13-85	IGKJ2	1
6001	IGKC	CQQYNSYPLTF	IGKV6-15	IGKJ5	1
6002	IGKC	CQQYNSYPLTF	IGKV6-15	IGKJ5	1
6006	IGKC	CQQYNSYPFTF	IGKV6-15	IGKJ4	2
6107	IGLC1	CALWYSTIWVF	IGLV1	IGLJ1	1

3254 rows × 5 columns

a = np.histogram([len(x) for x in df_out.cdr3_aa1],bins=np.arange(25))

plt.style.use('ggplot')
fig, ax = plt.subplots()

ax.bar(a[1][0:24],a[0])
ax.set_ylabel('num')
ax.set_xlabel('lengths')
#ax.set_title('CDR3 lengths distribute')
plt.show()

IGH CDR3 lengths distribute

df_out_H = df_out[[x.startswith('IGH') for x in df_out.const1]]

a = np.histogram([len(x) for x in df_out_H.cdr3_aa1],bins=np.arange(25))

plt.style.use('ggplot')
fig, ax = plt.subplots()

ax.bar(a[1][0:24],a[0])
ax.set_ylabel('num')
ax.set_xlabel('lengths')
#ax.set_title('CDR3 lengths distribute')
plt.show()

IGK & IGL CDR3 lengths distribute

df_out_L = df_out[[x.startswith('IGK') or x.startswith('IGL')  for x in df_out.const1]]

a = np.histogram([len(x) for x in df_out_L.cdr3_aa1],bins=np.arange(25))

plt.style.use('ggplot')
fig, ax = plt.subplots()

ax.bar(a[1][0:24],a[0])
ax.set_ylabel('num')
ax.set_xlabel('lengths')
#ax.set_title('CDR3 lengths distribute')
plt.show()

CDR3 abundance

df_out_2 = df_out.loc[:,["cdr3_aa1","n"]]
abu = df_out_2.groupby(df_out_2.cdr3_aa1).sum()

abu = abu.sort_values(by=["n"],ascending=False)

plt.style.use('ggplot')
fig, ax = plt.subplots()

ax.bar(np.arange(20),abu.n[0:20])
ax.set_xticklabels([])
plt.show()

V gene usage

重链中V gene使用频率

df_out_3 = df_out_H.loc[:,["v_name1","n"]]
abu = df_out_3.groupby(df_out_3.v_name1).sum()

abu = abu.sort_values(by=["n"],ascending=False)
abu

	n
v_name1
IGHV1-19	1582
IGHV14-3	1514
IGHV1-69	1080
IGHV1-56	980
IGHV2-5	917
...	...
IGHV1-18	2
IGHV9-4	2
IGHV5-9-3	1
IGHV1-43	1
IGHV1-42	1

81 rows × 1 columns

plt.style.use('ggplot')
fig, ax = plt.subplots()

ax.bar(np.arange(20),abu.n[0:20])
ax.set_xticklabels([])
plt.show()

轻链中V gene使用频率

df_out_3 = df_out_L.loc[:,["v_name1","n"]]
abu = df_out_3.groupby(df_out_3.v_name1).sum()

abu = abu.sort_values(by=["n"],ascending=False)
abu

	n
v_name1
IGKV8-30	2057
IGKV10-96	1550
IGKV12-46	1167
IGKV6-17	955
IGKV8-27	864
...	...
IGKV4-53	1
IGKV6-29	1
IGKV3-3	1
IGKV1-132	1
IGLV3	1

75 rows × 1 columns

plt.style.use('ggplot')
fig, ax = plt.subplots()

ax.bar(np.arange(20),abu.n[0:20])
ax.set_xticklabels([])
plt.show()

J gene usage

重链中J gene使用频率

df_out_4 = df_out_H.loc[:,["j_name1","n"]]
abu = df_out_4.groupby(df_out_4.j_name1).sum()
abu = abu.sort_values(by=["n"],ascending=False)
abu

	n
j_name1
IGHJ3	4409
IGHJ4	2169
IGHJ2	1739
IGHJ1	1233

plt.style.use('ggplot')
fig, ax = plt.subplots()

ax.bar(np.arange(len(abu.n)),abu.n[0:])
ax.set_xticklabels([])
plt.show()

轻链中J gene使用频率

df_out_4 = df_out_L.loc[:,["j_name1","n"]]
abu = df_out_4.groupby(df_out_4.j_name1).sum()
abu = abu.sort_values(by=["n"],ascending=False)
abu

	n
j_name1
IGKJ2	5702
IGKJ1	2143
IGKJ5	1355
IGKJ4	624
IGLJ1	24
IGLJ2	17
IGLJ3	3

plt.style.use('ggplot')
fig, ax = plt.subplots()

ax.bar(np.arange(len(abu.n)),abu.n[0:])
ax.set_xticklabels([])
plt.show()

V_J pairs

主要是看V gene和J gene联用频率

df_pair = df_out.groupby(["v_name1","j_name1"]).sum()
df_pair

		n
v_name1	j_name1
IGHV1-12	IGHJ3	2
IGHV1-14	IGHJ1	4
	IGHJ2	27
	IGHJ3	4
	IGHJ4	33
...	...	...
IGKV9-129	IGKJ4	3
IGLV1	IGLJ1	24
	IGLJ3	3
IGLV2	IGLJ2	16
IGLV3	IGLJ2	1

361 rows × 1 columns

df_pair.sort_values(by="n",ascending=False)

		n
v_name1	j_name1
IGKV8-30	IGKJ2	2041
IGHV1-19	IGHJ3	1575
IGKV10-96	IGKJ2	1520
IGHV14-3	IGHJ3	1504
IGHV1-69	IGHJ1	1015
...	...	...
IGHV5-6-3	IGHJ1	1
IGHV5-6-5	IGHJ3	1
IGKV1-132	IGKJ1	1
IGHV5-9-3	IGHJ2	1
IGLV3	IGLJ2	1

361 rows × 1 columns