10X单细胞（10X空间转录组）转录组 & VDJ 联合分析（14）之CoNGA

hello，今天我们继续来分享分析代码部分，这一次，我们直接用单细胞Seurat 的分析结果联合VDJ进行分析，话不多说，直接开始。

import scanpy as sc
import numpy as np
import anndata2ri as ad
# Activate the anndata2ri conversion between SingleCellExperiment and AnnData
ad.activate()

#Loading the rpy2 extension enables cell magic to be used
#This runs R code in jupyter notebook cells
%load_ext rpy2.ipython

sc.settings.verbosity = 3
sc.logging.print_versions()
%%R
# IMPORTANT! Use Seurat v3.0.2 as subsequent versions are packaged with "Reticulate" which interferes with rpy2
library('Seurat') 
library('base')
library('patchwork')

#read your Seurat object in
V1 <- readRDS('../conga_datasets/processed/vdj_v1_hs_V1_sc_5gex.rds')

DimPlot(V1) | FeaturePlot( V1, "CD3E")
# Keep in mind, conga will keep only the T cells with paired TCR information

写出h5ad

%%R -o V1_sce
# R-magics line above tells anndata2ri to convert the dataset to AnnData

#  but first convert your object to SingleCellExperiment Format 
V1_sce <- as.SingleCellExperiment(V1)
V1_sce
#write to h5ad for now. Ideally be nice to get away from strictly using the path
V1_sce.write("../conga_datasets/processed/vdj_v1_hs_V1_sc_5gex.h5ad")

开始分析,数据大家换成自己的数据

%matplotlib inline

#begin conga workflow
import conga
import pandas as pd
import matplotlib.pyplot as plt
from conga.tcrdist.make_10x_clones_file import make_10x_clones_file


# call the Anndata file and conga on
gex_datafile = "../conga_datasets/processed/vdj_v1_hs_V1_sc_5gex.h5ad"
gex_datatype = 'h5ad' # other possibilities right now: ['10x_mtx', 'h5ad'] (h5ad from scanpy)
tcr_datafile = './testDrive/vdj_v1_hs_pbmc_t_filtered_contig_annotations.csv'
organism = 'human'
clones_file = 'tmp_hs_pbmc_clones.tsv'
outfile_prefix = 'tmp_hs_pbmc'
# this creates the TCRdist 'clones file'
make_10x_clones_file( tcr_datafile, organism, clones_file )
# this command will create another file with the kernel PCs for subsequent reading by conga
conga.preprocess.make_tcrdist_kernel_pcs_file_from_clones_file( clones_file, organism )

数据联合

adata = conga.preprocess.read_dataset(gex_datafile, gex_datatype, clones_file )

# store the organism info in adata
adata.uns['organism'] = organism

# top 50 TCRdist kPCS 
adata.obsm['X_pca_tcr'].shape

(2287, 50)

# CDR3-alpha regions:
adata.obs['cdr3a'].head(3)

AAACCTGAGTCTTGCA-1 CAQSDSNYQLIW
AAACCTGCAGATGGGT-1 CAVLTNDMRF
AAACCTGTCCTTGCCA-1 CATDFNTGANSKLTF
Name: cdr3a, dtype: object

前处理

adata = conga.preprocess.filter_and_scale( adata )

挑选单个克隆的细胞

adata = conga.preprocess.reduce_to_single_cell_per_clone( adata )
adata

compute pca to find rep cell for each clone (1630, 1324)
computing PCA
on highly variable genes
with n_comps=50
finished (0:00:00)
num_clones: 1539
Normalize and logging matrix...
unable to find feature_types varname
Index(['HES4', 'TNFRSF4', 'RP3-395M20.12', 'MEGF6', 'ACOT7', 'RP1-202O8.3',
'RP11-312B8.1', 'UTS2', 'TNFRSF9', 'SPSB1',
...
'PKNOX1', 'AP001055.6', 'ITGB2', 'ITGB2-AS1', 'COL18A1', 'COL6A1',
'COL6A2', 'C21orf58', 'S100B', 'MT-ND6'],
dtype='object', length=1324)
found 489 of 512 good_genes in var_names organism=human
reduce from 1630 cells to 1539 cells (one per clonotype)
normalize_and_log_the_raw_matrix:: matrix is already logged
AnnData object with n_obs × n_vars = 1539 × 1324
obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'RNA_snn_res.0.8', 'seurat_clusters', 'ident', 'va', 'ja', 'cdr3a', 'cdr3a_nucseq', 'vb', 'jb', 'cdr3b', 'cdr3b_nucseq', 'n_genes', 'percent_mito', 'n_counts', 'clone_sizes'
var: 'vst.mean', 'vst.variance', 'vst.variance.expected', 'vst.variance.standardized', 'vst.variable', 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'mean', 'std'
uns: 'organism', 'log1p', 'pca', 'raw_matrix_is_logged', 'X_igex_genes'
obsm: 'X_umap', 'X_pca_tcr', 'X_igex'
varm: 'PCs'
layers: 'logcounts'

降维聚类

adata = conga.preprocess.cluster_and_tsne_and_umap( adata )

adata

computing PCA
on highly variable genes
with n_comps=50
finished (0:00:00)
computing neighbors
using 'X_pca' with n_pcs = 40
finished: added to .uns['neighbors']
.obsp['distances'], distances for each pair of neighbors
.obsp['connectivities'], weighted adjacency matrix (0:00:01)
computing UMAP
finished: added
'X_umap', UMAP coordinates (adata.obsm) (0:00:02)
running Louvain clustering
using the "louvain" package of Traag (2017)
finished: found 8 clusters and added
'louvain_gex', the cluster labels (adata.obs, categorical) (0:00:00)
ran louvain clustering: louvain_gex
computing neighbors
using 'X_pca' with n_pcs = 50
finished: added to .uns['neighbors']
.obsp['distances'], distances for each pair of neighbors
.obsp['connectivities'], weighted adjacency matrix (0:00:00)
computing UMAP
finished: added
'X_umap', UMAP coordinates (adata.obsm) (0:00:02)
running Louvain clustering
using the "louvain" package of Traag (2017)
finished: found 10 clusters and added
'louvain_tcr', the cluster labels (adata.obs, categorical) (0:00:00)
ran louvain clustering: louvain_tcr
AnnData object with n_obs × n_vars = 1539 × 1324
obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'RNA_snn_res.0.8', 'seurat_clusters', 'ident', 'va', 'ja', 'cdr3a', 'cdr3a_nucseq', 'vb', 'jb', 'cdr3b', 'cdr3b_nucseq', 'n_genes', 'percent_mito', 'n_counts', 'clone_sizes', 'louvain_gex', 'clusters_gex', 'louvain_tcr', 'clusters_tcr'
var: 'vst.mean', 'vst.variance', 'vst.variance.expected', 'vst.variance.standardized', 'vst.variable', 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'mean', 'std'
uns: 'organism', 'log1p', 'pca', 'raw_matrix_is_logged', 'X_igex_genes', 'neighbors', 'umap', 'louvain'
obsm: 'X_pca_tcr', 'X_igex', 'X_pca_gex', 'X_umap_gex', 'X_gex_2d', 'X_umap_tcr', 'X_tcr_2d'
varm: 'PCs'
layers: 'logcounts'
obsp: 'distances', 'connectivities'

联合分析

plt.figure(figsize=(12,6))
plt.subplot(121)
xy = adata.obsm['X_gex_2d']
clusters = np.array(adata.obs['clusters_gex'])
cmap = plt.get_cmap('tab20')
colors = [ cmap.colors[x] for x in clusters]
plt.scatter( xy[:,0], xy[:,1], c=colors)
plt.title('GEX UMAP colored by GEX clusters')

plt.subplot(122)
xy = adata.obsm['X_tcr_2d']
clusters = np.array(adata.obs['clusters_tcr'])
cmap = plt.get_cmap('tab20')
colors = [ cmap.colors[x] for x in clusters]
plt.scatter( xy[:,0], xy[:,1], c=colors)
plt.title('TCR UMAP colored by TCR clusters');
# these are the nbrhood sizes, as a fraction of the entire dataset:
nbr_fracs = [0.01, 0.1]

# we use this nbrhood size for computing the nndists
nbr_frac_for_nndists = 0.01

all_nbrs, nndists_gex, nndists_tcr = conga.preprocess.calc_nbrs(
    adata, nbr_fracs, also_calc_nndists=True, nbr_frac_for_nndists=nbr_frac_for_nndists)

# stash these in obs array, they are used in a few places...                                                                                                                
adata.obs['nndists_gex'] = nndists_gex
adata.obs['nndists_tcr'] = nndists_tcr

conga.preprocess.setup_tcr_cluster_names(adata) #stores in adata.uns

graph_vs_graph

results = conga.correlations.run_graph_vs_graph(adata, all_nbrs)
conga_scores = adata.obs['conga_scores']
good_mask = ( conga_scores <= 1.0 )
adata.obs['good_score_mask'] = good_mask
print(f'found {np.sum(good_mask)} conga hits')
results.sort_values('conga_score', inplace=True)
results.head(3)

conga_score	num_neighbors_gex	num_neighbors_tcr	overlap	overlap_corrected	mait_fraction	clone_index	nbr_frac	overlap_type	num_neighbors	cluster_size
5.722190e-49	NaN	NaN	72	59	0.888889	103	0.1	cluster_nbr	153.0	110.0
1.738370e-47	NaN	NaN	71	58	0.887324	77	0.1	cluster_nbr	153.0	110.0
1.738370e-47	NaN	NaN	71	58	0.887324	52	0.1	cluster_nbr	153.0	110.0

# write the results to a tsv file
clusters_gex = np.array(adata.obs['clusters_gex'])
clusters_tcr = np.array(adata.obs['clusters_tcr'])

indices = results['clone_index']
results['gex_cluster'] = clusters_gex[indices]
results['tcr_cluster'] = clusters_tcr[indices]
for tag in 'va ja cdr3a vb jb cdr3b'.split():
    results[tag] = list(adata.obs[tag][indices])
tsvfile = outfile_prefix+'_graph_vs_graph_hits.tsv'
print('saving graph-vs-graph results to file:',tsvfile)

results.to_csv(tsvfile, sep='\t', index=False)
#put the conga hits on top
colors = np.sqrt(np.maximum(-1*np.log10(conga_scores),0.0))
reorder = np.argsort(colors)

plt.figure(figsize=(12,6))
plt.subplot(121)
xy = adata.obsm['X_gex_2d']
plt.scatter( xy[reorder,0], xy[reorder,1], c=colors[reorder], vmin=0, vmax=np.sqrt(5))
plt.title('GEX UMAP colored by conga score')

plt.subplot(122)
xy = adata.obsm['X_tcr_2d']
plt.scatter( xy[reorder,0], xy[reorder,1], c=colors[reorder], vmin=0, vmax=np.sqrt(5))
plt.title('TCR UMAP colored by conga score');

最后的结果

nbrs_gex, nbrs_tcr = all_nbrs[0.1]

min_cluster_size = 5

# calc tcr sequence features of good cluster pairs
good_bicluster_tcr_scores = conga.correlations.calc_good_cluster_tcr_features(
    adata, good_mask, clusters_gex, clusters_tcr, conga.tcr_scoring.all_tcr_scorenames, verbose=False,
    min_count=min_cluster_size)

# run rank_genes on most common biclusters
rank_genes_uns_tag = 'rank_genes_good_biclusters'
conga.correlations.run_rank_genes_on_good_biclusters(
    adata, good_mask, clusters_gex, clusters_tcr, min_count=min_cluster_size, key_added= rank_genes_uns_tag)

gex_header_tcr_score_names = ['mhci2', 'cdr3len', 'cd8', 'nndists_tcr']

logo_pngfile = outfile_prefix+'_bicluster_logos.png'

conga.plotting.make_logo_plots(
    adata, nbrs_gex, nbrs_tcr, min_cluster_size, logo_pngfile,
    good_bicluster_tcr_scores=good_bicluster_tcr_scores,
    rank_genes_uns_tag = rank_genes_uns_tag,
    gex_header_tcr_score_names = gex_header_tcr_score_names )

pval_threshold = 1.
results = []
for nbr_frac in nbr_fracs:
    nbrs_gex, nbrs_tcr = all_nbrs[nbr_frac]
    print('finding biased GEX features for nbrhoods with size', nbr_frac, nbrs_gex.shape)
    results.append( conga.correlations.tcr_nbrhood_rank_genes_fast( adata, nbrs_tcr, pval_threshold, verbose=False))
    results[-1]['nbr_frac'] = nbr_frac

# save the results to a file
tsvfile = outfile_prefix+'_tcr_nbr_graph_vs_gex_features.tsv'
print('making:', tsvfile)
results_df = pd.concat(results, ignore_index=True)
results_df.to_csv(tsvfile, index=False, sep='\t')

pngfile = outfile_prefix+'_tcr_nbr_graph_vs_gex_features.png'
print('making:', pngfile)
conga.plotting.plot_ranked_strings_on_cells(
    adata, results_df, 'X_tcr_2d', 'clone_index', 'mwu_pvalue_adj', 1.0, 'feature', pngfile)

pngfile = outfile_prefix+'_tcr_nbr_graph_vs_gex_features_panels.png'
print('making:', pngfile)
conga.plotting.make_feature_panel_plots(adata, 'tcr', all_nbrs, results_df, pngfile, 'gex')

pval_threshold = 1.
results = []
tcr_score_names = conga.tcr_scoring.all_tcr_scorenames # the TCR features to use
for nbr_frac in nbr_fracs:
    nbrs_gex, nbrs_tcr = all_nbrs[nbr_frac]
    results.append( conga.correlations.gex_nbrhood_rank_tcr_scores(
        adata, nbrs_gex, tcr_score_names, pval_threshold, verbose=False ))
    results[-1]['nbr_frac'] = nbr_frac
results_df = pd.concat(results, ignore_index=True)

tsvfile = outfile_prefix+'_gex_nbr_graph_vs_tcr_features.tsv'
print('making:', tsvfile)
results_df.to_csv(tsvfile, index=False, sep='\t')

pngfile = outfile_prefix+'_gex_nbr_graph_vs_tcr_features.png'
print('making:', pngfile)

conga.plotting.plot_ranked_strings_on_cells(
    adata, results_df, 'X_gex_2d', 'clone_index', 'mwu_pvalue_adj', 1.0, 'feature', pngfile,
    direction_column='ttest_stat')

pngfile = outfile_prefix+'_gex_nbr_graph_vs_tcr_features_panels.png'
print('making:', pngfile)
conga.plotting.make_feature_panel_plots(adata, 'gex', all_nbrs, results_df, pngfile, 'tcr')

代码奉上，生活很好，有你更好