<jupyter_start><jupyter_text>CPU analysis of single-cell ATAC-seq data from 60K cellsCopyright (c) 2020, NVIDIA CORPORATION.Licensed under the Apache License, Version 2.0 (the "License") you may not use this file except in compliance with the License. You may obtain a copy of the License athttp://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.  Import requirements<jupyter_code>%%time
import numpy as np
import scanpy as sc
import anndata
import matplotlib.pyplot as plt
from collections import OrderedDict

import time
import os, wget

import utils

import warnings
warnings.filterwarnings('ignore', 'Expected ')
warnings.simplefilter('ignore')<jupyter_output>CPU times: user 1.95 s, sys: 3.21 s, total: 5.17 s
Wall time: 1.76 s<jupyter_text>Set parameters<jupyter_code># filtering peaks
n_top_peaks = 25000 # Number of peaks to retain

# PCA
n_components = 50 # Number of principal components to compute

# t-SNE
tsne_n_pcs = 20 # Number of principal components to use for t-SNE

# KNN
n_neighbors = 10 # Number of nearest neighbors for KNN graph
knn_n_pcs = 50 # Number of principal components to use for finding nearest neighbors

# UMAP
umap_min_dist = 0.4 
umap_spread = 1.0

# Differential peaks
n_diff_peaks = 50 # Number of differential peaks to compute for each cluster

# Marker genes
# Gene coordinates are in GRCh37, from grch37.ensembl.org
markers = ['MS4A1', 'LEF1', 'NKG7', 'TREM1', 'GATA1', 'PAX5']
marker_coords = OrderedDict()
marker_coords['MS4A1'] = ['chr11', 60223225, 60238233, 'fwd']
marker_coords['LEF1'] = ['chr4', 108968701, 109090112, 'rev']
marker_coords['NKG7'] = ['chr19', 51874860, 51875969, 'rev']
marker_coords['TREM1'] = ['chr6', 41235664, 41254457, 'rev']
marker_coords['GATA1'] = ['chrX', 48644962, 48652716, 'fwd']
marker_coords['PAX5'] = ['chr9', 36833272, 37034103, 'rev']

# Number of parallel jobs
sc.settings.n_jobs=32
start_time = time.time()<jupyter_output><empty_output><jupyter_text>Input data<jupyter_code>%%time
input_file = "../data/dsci_resting_nonzeropeaks.h5ad"

if not os.path.exists(input_file):
    print('Downloading import file...')
    os.makedirs('../data', exist_ok=True)
    wget.download('https://rapids-single-cell-examples.s3.us-east-2.amazonaws.com/dsci_resting_nonzeropeaks.h5ad',
                  input_file)

adata = anndata.read_h5ad(input_file)
adata_raw = adata.copy()<jupyter_output>CPU times: user 302 ms, sys: 365 ms, total: 668 ms
Wall time: 666 ms<jupyter_text>Preprocessing<jupyter_code>preprocessing_start = time.time()<jupyter_output><empty_output><jupyter_text>tf-idf normalizationWe use log1p(pseudocount / tf)\*idf, motivated by discussion in [this blog post](http://andrewjohnhill.com/blog/2019/05/06/dimensionality-reduction-for-scatac-data/) as well as our results on various scATAC datasets.<jupyter_code>%%time
lognormalized = utils.logtf_idf(adata.X, pseudocount=10**3)
adata.X = lognormalized<jupyter_output>CPU times: user 2.67 s, sys: 684 ms, total: 3.35 s
Wall time: 3.35 s<jupyter_text>Frequency-based peak selection<jupyter_code>%%time
adata = utils.filter_peaks(adata, n_top_peaks)
print(adata.X.shape)
print("Preprocessing time: %.2fsec" % (time.time() - preprocessing_start))<jupyter_output>Preprocessing time: 4.64sec<jupyter_text>Cluster & Visualize  PCA<jupyter_code>%%time
sc.pp.pca(adata, n_comps = n_components)<jupyter_output>CPU times: user 13min 35s, sys: 1h 14s, total: 1h 13min 49s
Wall time: 3min 10s<jupyter_text>UMAP visualization<jupyter_code>%%time
sc.pp.neighbors(adata, n_neighbors=n_neighbors, n_pcs=knn_n_pcs)
%%time
sc.tl.umap(adata, min_dist=umap_min_dist, spread=umap_spread)<jupyter_output>CPU times: user 55.1 s, sys: 33 s, total: 1min 28s
Wall time: 33 s<jupyter_text>Graph clustering<jupyter_code>%%time
sc.tl.louvain(adata)
%%time
sc.tl.leiden(adata)
%%time

fig, axs = plt.subplots(nrows = 1, ncols = 3, figsize = (12, 3), dpi=100)

sc.pl.umap(adata, color=['louvain'], ax=axs[0], show=False)
sc.pl.umap(adata, color=['leiden'], ax=axs[1], show=False)
sc.pl.umap(adata, color=['cell_type'], ax=axs[2], show=False)

plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>tSNE visualization<jupyter_code>%%time
sc.tl.tsne(adata, n_pcs=tsne_n_pcs)
%%time

fig, axs = plt.subplots(nrows = 1, ncols = 3, figsize = (12, 3), dpi=100)

sc.pl.tsne(adata, color=['louvain'], ax=axs[0], show=False)
sc.pl.tsne(adata, color=['leiden'], ax=axs[1], show=False)
sc.pl.tsne(adata, color=['cell_type'], ax=axs[2], show=False)

plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Visualize marker gene activities We use the raw adata object to compute gene activity scores. This allows us to consider all peaks, not just the subset we used for clustering and UMAP.<jupyter_code>%%time
# find which peaks overlap with each marker gene (+ 5kb upstream)
peaks_per_gene = {gene : [] for gene in markers}

for (i, peak) in enumerate(adata_raw.var_names):
    chrom, coords = peak.strip().split(':')
    start, end = [int(c) for c in coords.split('-')]
    for gene in markers:
        if utils.overlap(marker_coords[gene], [chrom, start, end], 5000, 0):
            peaks_per_gene[gene].append((i, peak))

for gene in peaks_per_gene:
    print(f"{gene} has {len(peaks_per_gene[gene])} nearby peaks")
%%time
# Calculate per-cell gene activity scores for each marker gene
cell_gene_scores = np.zeros((adata.X.shape[0], len(marker_coords.keys())))
for (i, gene) in enumerate(peaks_per_gene.keys()):
    indices = [x[0] for x in peaks_per_gene[gene]]
    peak_scores = adata_raw.X[:, indices].toarray()
    cell_gene_scores[:, i] = np.sum(peak_scores, axis=1) + 1

# Store per-cell gene activity scores with adata
for (i, gene) in enumerate(markers):
    adata.obs[gene] = np.log2(cell_gene_scores[:, i].tolist())
    adata.obs.loc[np.isneginf(adata.obs[gene]), gene] = 0
%%time

fig, axs = plt.subplots(nrows=2, ncols=3, figsize = (8, 4), dpi=100)
axs = axs.flatten()
for (i, gene) in enumerate(markers):
    sc.pl.umap(adata, color_map="Blues", color=gene, ax = axs[i], show=False)
plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Find differential peaks We use a logistic regression model to find peaks that show differential activity across clusters. Here, we needed to increase the maximum number of iterations to 1200 for the model to converge.<jupyter_code>%%time
sc.tl.rank_genes_groups(adata, 
                        groupby="louvain", 
                        n_genes=n_diff_peaks, 
                        groups='all', 
                        reference='rest',
                        method='logreg',
                        max_iter=1200
                       )
%%time

fig, axs = plt.subplots(len(adata.uns['rank_genes_groups']['names'][0])//3 + 1, 3, figsize=(9,9), dpi=100)
axs = axs.flatten()

for (i, peak) in enumerate(adata.uns['rank_genes_groups']['names'][0]):
    sc.pl.umap(adata, color=peak, 
           ax=axs[i], show=False,
           vmax=min(25, np.max(adata[:,peak].X)),
           cmap='viridis'
          )
plt.tight_layout()
plt.show()
print("Full time: %s" % (time.time() - start_time))<jupyter_output>Full time: 1561.1989142894745
{'fraction_non_alphanumeric': 0.13243606998654106, 'fraction_numerical': 0.04226110363391655, 'mean_word_length': 4.883610451306414, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 42, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '430760', 'license': 'Apache-2.0', 'n_tokens_mistral': 2972, 'n_tokens_neox': 2675, 'n_words': 673, 'path': 'notebooks/dsci_bmmc_60k_cpu.ipynb', 'repo_name': 'trivialfis/rapids-single-cell-examples', 'size': '7430'}