<jupyter_start><jupyter_text>Prelim Dataset downloaded from : https://support.10xgenomics.com/single-cell-multiome-atac-gex/datasets/1.0.0/pbmc_unsorted_10k Data is available at `s3://fh-pi-setty-m-eco-public/single-cell-primers/multiome/`Review the notebook `PBMC-RNA-standalone.ipynb` for setup instructions.  Load data Download all the multiome data using the following command<jupyter_code>!aws s3 sync s3://fh-pi-setty-m-eco-public/single-cell-primers/multiome/ data/multiome/<jupyter_output>download: s3://fh-pi-setty-m-eco-public/single-cell-primers/multiome/ArchR/pbmc_multiome_atac/export/peak_counts/cells.csv to data/multiome/ArchR/pbmc_multiome_atac/export/peak_counts/cells.csv
download: s3://fh-pi-setty-m-eco-public/single-cell-primers/multiome/ArchR/pbmc_multiome_atac/export/cell_metadata.csv to data/multiome/ArchR/pbmc_multiome_atac/export/cell_metadata.csv
download: s3://fh-pi-setty-m-eco-public/single-cell-primers/multiome/ArchR/pbmc_multiome_atac/export/peak_counts/peaks.csv to data/multiome/ArchR/pbmc_multiome_atac/export/peak_counts/peaks.csv
download: s3://fh-pi-setty-m-eco-public/single-cell-primers/multiome/ArchR/pbmc_multiome_atac/export/svd.csv to data/multiome/ArchR/pbmc_multiome_atac/export/svd.csv
download: s3://fh-pi-setty-m-eco-public/single-cell-primers/multiome/ArchR/pbmc_multiome_atac/export/peak_counts/counts.mtx to data/multiome/ArchR/pbmc_multiome_atac/export/peak_counts/counts.mtx
download: s3://fh-pi-setty-m-eco-public/single-cell-primers/mult[...]<jupyter_text>Imports<jupyter_code>import os
import pandas as pd
import numpy as np
import scanpy as sc
import pyranges as pr
import warnings
import palantir 
import phenograph
import harmony
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

sns.set_style('ticks')
matplotlib.rcParams['figure.figsize'] = [4, 4]
matplotlib.rcParams['figure.dpi'] = 100
matplotlib.rcParams['image.cmap'] = 'Spectral_r'
warnings.filterwarnings(action="ignore", module="matplotlib", message="findfont")<jupyter_output><empty_output><jupyter_text>Utility functions<jupyter_code>def log_transform(ad, ps=0.1):
    ad.X.data = np.log2(ad.X.data + ps) - np.log2(ps)
from joblib import Parallel, delayed
from scipy.sparse import hstack, csr_matrix, issparse
import gc

def _dot_func(x, y):
    return x.dot(y)


def impute_data(dm_res, ad, n_steps=3):
    T_steps = dm_res['T'] ** n_steps
    T_steps = T_steps.astype(np.float32)

    # RUn in parallel
    seq = np.append(np.arange(0, ad.X.shape[1], 100), [ad.X.shape[1]])
    res = Parallel(n_jobs=-1)(delayed(_dot_func)(T_steps, ad.X[:, seq[i - 1]:seq[i]]) for i in range(1, len(seq)))
    imputed_data = hstack(res)
    imputed_data = imputed_data.todense()
    imputed_data[imputed_data < 1e-2] = 0
    gc.collect()

    return imputed_data<jupyter_output><empty_output><jupyter_text>Load data<jupyter_code>data_dir = 'data/multiome/'
# Load 10x filtered matrix containing both RNA and ATAC features 
comb_ad = sc.read_10x_h5(data_dir + 'pbmc_unsorted_10k_filtered_feature_bc_matrix.h5', gex_only=False)
comb_ad.obs_names = 'pbmc_multiome#' + comb_ad.obs_names
comb_ad
# Barcode metrics
bc_metrics = pd.read_csv(data_dir + '/pbmc_unsorted_10k_per_barcode_metrics.csv', 
                index_col=0)
bc_metrics.index = 'pbmc_multiome' + '#' +bc_metrics.index 

# FRIP
bc_metrics['FRIP'] = bc_metrics['atac_peak_region_fragments']/bc_metrics['atac_fragments']
# Separate RNA and ATAC components
rna_ad = comb_ad[:, comb_ad.var['feature_types'] == 'Gene Expression']
rna_ad.var_names_make_unique()

atac_ad = comb_ad[:, comb_ad.var['feature_types'] == 'Peaks']
atac_ad.var_names_make_unique()
rna_ad
atac_ad
sample_cells = rna_ad.obs_names<jupyter_output><empty_output><jupyter_text>Preprocess This section follows the preprocessing steps of the single-cell RNA-seq analysis<jupyter_code>ad = rna_ad<jupyter_output><empty_output><jupyter_text>QC<jupyter_code>ad.var['mt'] = ad.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(ad, qc_vars=['mt'], inplace=True, percent_top=None, log1p=False, )<jupyter_output><empty_output><jupyter_text>Filter cell based on RNA<jupyter_code>plt.figure()
ax = plt.gca()
ax.hist(np.log10(ad.obs['total_counts']), 50)
ax.set_xlabel('log10(molecules per cell)')
ax.set_ylabel('Frequency')
sns.despine()
# Manually adjust the threshold percentage depending on the dataset
# Doublet detection (Eg: scrublet) will work better for higher molecule count cells
THRESHOLDS = np.log10(np.percentile((ad.obs['total_counts']), [3.5, 99.5]))

plt.figure()
ax = plt.gca()
ax.hist(np.log10(ad.obs['total_counts']), 50)
ax.vlines(THRESHOLDS, ax.get_ylim()[0], ax.get_ylim()[1], color='black', linestyle='--')
ax.set_xlabel('log10(molecules per cell)')
ax.set_ylabel('Frequency')
sns.despine()
# Filter cells 
sample_cells = sample_cells[(np.log10(ad.obs['total_counts'][sample_cells]) > THRESHOLDS[0]) & \
    (np.log10(ad.obs['total_counts'][sample_cells]) < THRESHOLDS[1]) ]

plt.figure()
ax = plt.gca()
ax.hist(np.log10(ad.obs['total_counts'][sample_cells]), 50)
ax.set_xlabel('log10(molecules per cell)')
ax.set_ylabel('Frequency')
sns.despine()<jupyter_output><empty_output><jupyter_text>Filter based on ATAC<jupyter_code>## CAUTION: CellRanger peak calling is unreliable and we only use it here to filter cells with low 
##          fraction of reads in peaks. This is ok since CellRanger peak caller dramatically over estimates 
##          the peak widths. We might incur some false negatives, but thats an acceptable error for filtering   
plt.figure()
ax = plt.gca()
ax.hist(bc_metrics['FRIP'][sample_cells], 50)
plt.title('Fraction of reads in peaks')
sns.despine()
THRESHOLD = 0.4
sample_cells = sample_cells[bc_metrics['FRIP'][sample_cells] > THRESHOLD]<jupyter_output><empty_output><jupyter_text>Filter based on mitochondria<jupyter_code>MITO_THRESHOLD = 25

x,y,dens = palantir.plot.density_2d(ad.obs['total_counts'][sample_cells], 
               ad.obs['pct_counts_mt'][sample_cells])

plt.figure()
ax = plt.gca()
ax.scatter(x, y, c=dens, s=5)
ax.hlines(MITO_THRESHOLD, ax.get_xlim()[0], ax.get_xlim()[1], color='black', linestyle='--')

plt.xlabel('Total molecule count')
plt.ylabel('Mitochondrial fraction')
sns.despine()

sample_cells = sample_cells[ad.obs['pct_counts_mt'][sample_cells] < MITO_THRESHOLD]
ad = ad[sample_cells, :]
ad<jupyter_output><empty_output><jupyter_text>Analysis<jupyter_code># Filter genes 
sc.pp.filter_genes(ad, min_cells=50)
# Normalize and log transform
sc.pp.normalize_per_cell(ad)
log_transform(ad)
# Highly variable genes [Num. of genes is a paramter, typically 1.5-2.5k genes work well]
sc.pp.highly_variable_genes(ad, flavor='cell_ranger', n_top_genes=2500)
ad
# PCA
sc.pp.pca(ad, use_highly_variable=True, n_comps=50)
ad
# 50 comps is an approxmiation here - choose number of components by knee point or fraction of explained variance
# UMAP and Leiden - requires the computation of nearest neighbors
sc.pp.neighbors(ad, use_rep='X_pca')
sc.tl.umap(ad)
sc.tl.leiden(ad)
ad
# Diffusion maps 
dm_res = palantir.utils.run_diffusion_maps(pd.DataFrame(ad.obsm['X_pca'], index=ad.obs_names))
ad.obsp['DM_Kernel'] = dm_res['kernel']
ad.obsm['DM_EigenVectors'] = dm_res['EigenVectors'].values
ad.uns['DMEigenValues'] = dm_res['EigenValues'].values
ad
# Phenograph 
ad.obs['phenograph'], _, _ = phenograph.cluster(ad.obsm['X_pca'])
ad.obs['phenograph'] = ad.obs['phenograph'].astype(str)
# Imputation
ad.obsm['MAGIC_imputed_data'] = impute_data(dm_res, ad)<jupyter_output><empty_output><jupyter_text>Visualize results<jupyter_code>ad.obs['log_n_counts'] = np.log10(ad.obs['total_counts'])
sc.pl.scatter(ad, basis='umap', color='log_n_counts', frameon=False)
sc.pl.scatter(ad, basis='umap', color=['leiden'], frameon=False)<jupyter_output><empty_output><jupyter_text>Celltype annotation<jupyter_code>genes = pd.Series([
    "CD34",  #Early Progenitor
    "GATA1", #Erythroid
    "PAX5", "MS4A1", "MME", #B-Cell Trajectory
    "CD14",  #Monocytes
    "CD3D", "CD8A", 'CD4', #TCells,
    'FCGR3A', #'CD16', 
     'LEF1', 'NKG7', 'TREM1', 'LYZ', 'CD19',  'IRF8', 'IRF4',
    'ITGA2B', # CD41,
    'NCAM1', # CD56,
    'CCR7' # Memory T cell
])
genes = genes[genes.isin(ad.var_names)]
imp_df = pd.DataFrame(ad.obsm['MAGIC_imputed_data'], 
                                                index=ad.obs_names, columns=ad.var_names)


palantir.plot.plot_gene_expression(imp_df,
                                   pd.DataFrame(ad.obsm['X_umap'], 
                                               index=ad.obs_names, columns=['x', 'y']), 
                                   genes)<jupyter_output><empty_output><jupyter_text>Cell type annotation can be performed using the same markers listed in https://github.com/settylab/single-cell-primers/blob/main/notebooks/PBMC-RNA-standalone.ipynb  Save<jupyter_code>temp = comb_ad[ad.obs_names, :]
temp.var_names_make_unique()
temp = temp[:, ad.var_names]
# Attach raw counts
ad.raw = sc.AnnData(temp.X)
ad.write(data_dir + 'pbmc_multiome_rna.h5ad')
# Export cell names for ATAC - only this subset of cells will be use for ATAC analysis
pd.DataFrame(ad.obs_names).to_csv(data_dir + 'pbmc_multiome_cells.csv')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14129018245004343, 'fraction_numerical': 0.015204170286707211, 'mean_word_length': 5.508127208480565, 'pattern_counts': {'":': 0, '<': 55, '<?xml version=': 0, '>': 54, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '407011', 'license': 'MIT', 'n_tokens_mistral': 3746, 'n_tokens_neox': 3565, 'n_words': 587, 'path': 'notebooks/PBMC-multiome-RNA.ipynb', 'repo_name': 'settylab/single-cell-primers', 'size': '9208'}