Preprocessing and clustering scATAC PBMCs¶
The data consist in ~3000 cells of human PBMCs. Cells were FAC sorted thus providing a ground truth for cell type identification. The data was produces by Buenrostro et al. 2018 The original count matrix (using bulk peaks as features) from Buenrostro et al. is available here: GSE96772
On a unix system, you can uncomment and run the following to download the count matri in its anndata format.
In [1]:
# !mkdir data
# !wget https://www.dropbox.com/s/cwlaaxl70t27tb2/data_tutorial_buenrostro.tar.gz?dl=0 -O data/data_tutorial_buenrostro.tar.gz
# !cd data; tar -xzf data_tutorial_buenrostro.tar.gz
# !mkdir write
If episcanpy isn't installed already. Or if some of the functions aren't working. Install the following version:
In [2]:
# !pip install git+https://github.com/colomemaria/epiScanpy
Load libraries¶
In [3]:
import scanpy as sc
import anndata as ad
import numpy as np
import pandas as pd
import episcanpy.api as epi
In [4]:
sc.settings.set_figure_params(dpi=80, color_map='gist_earth')
In [5]:
results_file = 'write/buenrostro_pbmc.h5ad' # the file that will store the analysis results
Load the data¶
Read in the count matrix into an AnnData object, which holds many slots for annotations and different representations of the data. It also comes with its own HDF5 file format: .h5ad.
In [6]:
adata = ad.read('data/data_tutorial_buenrostro/all_buenrostro_bulk_peaks.h5ad')
adata
Out[6]:
In [7]:
# display information stored in adata.obs
adata.obs
Out[7]:
In [8]:
# checking the variable names (bulk peaks coordinates)
print(adata.var_names)
Extract the FACs information from the file names¶
In [9]:
adata.obs['facs_label'] = ['MEP' if 'MEP' in line else line.split('.bam')[0].lstrip('singles-').split('BM')[-1].split('-')[1] for line in adata.obs['cell_name'].tolist()]
adata
Out[9]:
Load the additional metadata¶
In [10]:
!head 'data/data_tutorial_buenrostro/metadata.tsv'
In [11]:
# format the cell names to match the metadata file
adata.obs_names = [x.split('/')[-1].split('.dedup.st.bam')[0] for x in adata.obs_names.tolist()]
adata.obs_names
Out[11]:
In [12]:
epi.pp.load_metadata(adata,
'data/data_tutorial_buenrostro/metadata.tsv',
separator='\t')
adata
Out[12]:
In [13]:
adata.obs['cell_type'] = [x.rstrip('\n') for x in adata.obs['cell_type\n']]
del adata.obs['cell_type\n']
In [14]:
adata
Out[14]:
In [15]:
# check the 2 different annotation obtained
pd.crosstab(adata.obs['cell_type'], adata.obs['facs_label'])
Out[15]:
Load gene/transcript annotation¶
Download annotation file (the data are aligned on hg19)
In [16]:
# !wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz -O data/data_tutorial_buenrostro/gencode.v19.annotation.gtf.gz
# !cd data/data_tutorial_buenrostro; gunzip gencode.v19.annotation.gtf
In [17]:
epi.tl.find_genes(adata,
gtf_file='data/data_tutorial_buenrostro/gencode.v19.annotation.gtf',
key_added='transcript_annotation',
upstream=2000,
feature_type='transcript',
annotation='HAVANA',
raw=False)
In [18]:
adata
Out[18]:
Preprocessing¶
Check if the data matrix is binary - if not, binarize the data matrix¶
In [19]:
print(np.max(adata.X))
In [20]:
# Here, the data is already binary. If not, run the following commands
# epi.pp.binarize(adata)
# print(np.max(adata.X))
Proceed to basic filtering.
In [21]:
adata
Out[21]:
In [22]:
# remove any potential empty features or barcodes
epi.pp.filter_cells(adata, min_features=1)
epi.pp.filter_features(adata, min_cells=1)
In [23]:
# filtered out 24210 peaks
adata
Out[23]:
Quality controls¶
In [24]:
adata.obs['log_nb_features'] = [np.log10(x) for x in adata.obs['nb_features']]
adata
Out[24]:
In [25]:
epi.pl.violin(adata, ['nb_features'])
epi.pl.violin(adata, ['log_nb_features'])
In [26]:
# set a minimum number of cells to keep
min_features = 1000
epi.pp.coverage_cells(adata, binary=True, log=False, bins=50,
threshold=min_features, save='Buenrostro_bulk_peaks_coverage_cells.png')
epi.pp.coverage_cells(adata, binary=True, log=10, bins=50,
threshold=min_features, save='Buenrostro_bulk_peaks_coverage_cells_log10.png')
In [27]:
# minimum number of cells sharing a feature
min_cells = 5
epi.pp.coverage_features(adata, binary=True, log=False,
threshold=min_cells, save='Buenrostro_bulk_peaks_coverage_peaks.png')
epi.pp.coverage_features(adata, binary=True, log=True,
threshold=min_cells, save='Buenrostro_bulk_peaks_coverage_peaks_log10.png')
Actually proceed to filter the cells and peaks based on the QC plots¶
In [28]:
min_features = 1000
epi.pp.filter_cells(adata, min_features=min_features)
adata
Out[28]:
In [29]:
min_cells = 5
epi.pp.filter_features(adata, min_cells=min_cells)
adata
Out[29]:
Looking at the QC plots after filtering¶
In [30]:
epi.pp.coverage_cells(adata, binary=True, log='log10', bins=50, threshold=min_features)
In [31]:
epi.pp.coverage_features(adata, binary=True, log='log10', bins=50, threshold=min_cells)
Identifying the most variable features¶
We aim to select a cuttof after the elbow.
In [32]:
epi.pp.cal_var(adata)
In [33]:
min_score_value = 0.515
nb_feature_selected = 120000
epi.pl.variability_features(adata,log=None,
min_score=min_score_value, nb_features=nb_feature_selected,
save='variability_features_plot_bonemarrow_peakmatrix.png')
epi.pl.variability_features(adata,log='log10',
min_score=min_score_value, nb_features=nb_feature_selected,
save='variability_features_plot_bonemarrow_peakmatrix_log10.png')
In [34]:
# save the current matrix in the raw layer
adata.raw = adata
In [35]:
# create a new AnnData containing only the most variable features
adata = epi.pp.select_var_feature(adata,
nb_features=nb_feature_selected,
show=False,
copy=True)
In [36]:
adata
Out[36]:
In [37]:
epi.pl.violin(adata, ['nb_features'])
epi.pl.violin(adata, ['log_nb_features'])
