Workflow for scATAC-Seq data#

If you haven’t already then we highly recommend that you please go the Workflow for scRNA-Seq analysis and also go through this vignette describing data organization and design principles of Scarf.

Here we will go through some of the basic steps of analysis of single-cell ATAC-seq (Assay for Transposase-Accessible Chromatin using sequencing) data.

%load_ext autotime

import scarf
scarf.__version__

'0.23.5'

time: 1.64 s (started: 2023-05-24 14:13:14 +00:00)

1) Fetch and convert data#

We will use 10x Genomics’s singel-cell ATAC-Seq data from peripheral blood mononuclear cells. Like single-cell RNA-Seq, Scarf only needs a count matrix to start the analysis of scATAC-Seq data. We can use fetch_dataset to download the data in 10x’s HDF5 format.

scarf.fetch_dataset(
    dataset_name='tenx_10K_pbmc-v1_atacseq',
    save_path='scarf_datasets'
)

time: 13.6 s (started: 2023-05-24 14:13:16 +00:00)

The CrH5Reader class provides access to the HDF5 file. We can load the file and quickly check the number of features, and also verify that Scarf identified the assay as an ATAC assay.

reader = scarf.CrH5Reader(
    'scarf_datasets/tenx_10K_pbmc-v1_atacseq/data.h5'
)
reader.assayFeats
ATAC
type Peaks
start 0
end 89796
nFeatures 89796 
time: 141 ms (started: 2023-05-24 14:13:30 +00:00)

writer = scarf.CrToZarr(
    reader,
    zarr_fn=f'scarf_datasets/tenx_10K_pbmc-v1_atacseq/data.zarr',
    chunk_size=(1000, 2000)
)
writer.dump(batch_size=1000)

time: 21.3 s (started: 2023-05-24 14:13:30 +00:00)

2) Create DataStore and filter cells#

We load the Zarr file on using DataStore class. The obtained DataStore object will be our single point on interaction for rest of this analysis. When loaded, Scarf will automatically calculate the number of cells where each peak is present and number of peaks that are accessible in cell (nFeatures). Scarf will also calculate the total number of fragments/cut sites within each cell.

ds = scarf.DataStore(
    'scarf_datasets/tenx_10K_pbmc-v1_atacseq/data.zarr', 
    nthreads=4
)

time: 14.6 s (started: 2023-05-24 14:13:51 +00:00)

We will use auto_filter_cells method which automatically remove outliers from the data. To identify outliers we generate a normal distribution using sample mean and variance. Using this normal distribution Scarf estimates the values with probability less 0.01 (default value) on both ends of the distribution and flags them for removal.

ds.auto_filter_cells()

INFO: 247 cells flagged for filtering out using attribute ATAC_nCounts

INFO: 294 cells flagged for filtering out using attribute ATAC_nFeatures

Matplotlib is building the font cache; this may take a moment.
../_images/342f160069b6a230e5cb04809d71a597da3258674eb3b47733f558ebab847284.png ../_images/70d2b3a81a3b72626c74653dce167856870a9a50b06fc2005a8bc41710683f66.png 
time: 17.7 s (started: 2023-05-24 14:14:06 +00:00)

3) Feature selection#

For scATAC-Seq data, the features are ranked by their TF-IDF normalized values, summed across all cells. The top n features are marked as prevalent_peaks and are used for downstream steps. Here we used top 25000 peaks, which makes for more than quarter of all the peaks which inline with the what has been suggested in other scATAC-Seq analysis protocols.

ds.mark_prevalent_peaks(top_n=25000)

time: 12.1 s (started: 2023-05-24 14:14:24 +00:00)

ds.ATAC.feats.head()
I ids names I__prevalent_peaks dropOuts nCells stats_I_prevalence
0 True chr1:565107-565550 chr1:565107-565550 False 8645 83 0.206181
1 True chr1:569174-569639 chr1:569174-569639 False 8529 199 0.350725
2 True chr1:713460-714823 chr1:713460-714823 True 6242 2486 1.993573
3 True chr1:752422-753038 chr1:752422-753038 True 8296 432 0.680701
4 True chr1:762106-763359 chr1:762106-763359 True 6850 1878 1.570790 
time: 59 ms (started: 2023-05-24 14:14:36 +00:00)

4) KNN graph creation#

For scATAC-Seq datasets, Scarf uses TF-IDF normalization. The normalization is automatically performed during the graph building step. The selected features, marked as prevalent_peaks in feature metadata, are used for graph creation. For the dimension reduction step, LSI (latent semantic indexing) is used rather than PCA. The rest of the steps are same as for scRNA-Seq data.

LSI reduction of scATAC-Seq is known to capture the sequencing depth of cells in the first LSI dimension. Hence, by default, the lsi_skip_first parameter is True but users can override it.

ds.make_graph(
    feat_key='prevalent_peaks',
    k=21,
    dims=50,
    n_centroids=100,
    lsi_skip_first=True,
)

/home/docs/checkouts/readthedocs.org/user_builds/scarf/envs/0.23.5/lib/python3.8/site-packages/sklearn/cluster/_kmeans.py:870: FutureWarning: The default value of `n_init` will change from 3 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  warnings.warn(

/home/docs/checkouts/readthedocs.org/user_builds/scarf/envs/0.23.5/lib/python3.8/site-packages/umap/distances.py:1063: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  @numba.jit()
/home/docs/checkouts/readthedocs.org/user_builds/scarf/envs/0.23.5/lib/python3.8/site-packages/umap/distances.py:1071: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  @numba.jit()
/home/docs/checkouts/readthedocs.org/user_builds/scarf/envs/0.23.5/lib/python3.8/site-packages/umap/distances.py:1086: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  @numba.jit()

/home/docs/checkouts/readthedocs.org/user_builds/scarf/envs/0.23.5/lib/python3.8/site-packages/umap/umap_.py:660: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  @numba.jit()

INFO: ANN recall: 99.41%

time: 4min 4s (started: 2023-05-24 14:14:36 +00:00)

5) UMAP reduction and clustering#

Non-linear dimension reduction using UMAP and tSNE are performed in the same way as for scRNA-Seq data. Because, in Scarf the core UMAP step are run directly on the neighbourhood graph, the scATAC-Seq data is handled similar to any other data.

ds.run_umap(
    n_epochs=500,
    min_dist=0.1, 
    spread=1, 
    parallel=True
)

time: 1min 11s (started: 2023-05-24 14:18:40 +00:00)

/home/docs/checkouts/readthedocs.org/user_builds/scarf/envs/0.23.5/lib/python3.8/site-packages/scarf/metadata.py:304: RuntimeWarning: overflow encountered in cast
  a = np.empty(self.N).astype(type(values[0]))

Same goes for clustering as well. The leiden clustering acts on the neighbourhood graph directly.

ds.run_leiden_clustering(resolution=0.6)

time: 514 ms (started: 2023-05-24 14:19:52 +00:00)

/home/docs/checkouts/readthedocs.org/user_builds/scarf/envs/0.23.5/lib/python3.8/site-packages/scarf/metadata.py:304: RuntimeWarning: invalid value encountered in cast
  a = np.empty(self.N).astype(type(values[0]))

Results from both UMAP and Leiden clustering are stored in cell attribute table. Here, because the UMAP and leiden clustering columns have ‘ATAC’ prefixes because they were executed on the default ‘ATAC’ assay. You can also see that the cells that were filtered out (marked by ‘False’ value in column ‘I’) have NaN value for UMAP axes and -1 for clustering

ds.cells.head()
I ids names ATAC_UMAP1 ATAC_UMAP2 ATAC_leiden_cluster ATAC_nCounts ATAC_nFeatures
0 True AAACGAAAGAGCGAAA-1 AAACGAAAGAGCGAAA-1 0.490977 13.089506 1 13187.0 5546.0
1 True AAACGAAAGAGTTTGA-1 AAACGAAAGAGTTTGA-1 3.260576 13.478303 1 16074.0 6586.0
2 True AAACGAAAGCGAGCTA-1 AAACGAAAGCGAGCTA-1 3.693170 -7.778605 2 28009.0 9214.0
3 False AAACGAAAGGCTTCGC-1 AAACGAAAGGCTTCGC-1 NaN NaN -1 221747.0 30325.0
4 True AAACGAAAGTGCTGAG-1 AAACGAAAGTGCTGAG-1 0.099497 13.407598 1 11440.0 4922.0 
time: 35.7 ms (started: 2023-05-24 14:19:52 +00:00)

We can visualize the UMAP embedding and the clusters of cells on the embedding

ds.plot_layout(
    layout_key='ATAC_UMAP', 
    color_by='ATAC_leiden_cluster'
)
../_images/47ca3cd1d15c4258c8f98ee53126f07b88e26b4f1625cd7e5667efa7f2f7fd53.png 
time: 842 ms (started: 2023-05-24 14:19:52 +00:00)

Those familiar with PBMC datasets might already be able to identify different cell types in the UMAP plot.

6) Calculating gene scores#

The features in snATAC-Seq in form of peak coordinates are hard to interpret by themselves. Hence, distilling the open chromatin regions in terms of accessible genes can help in identification of cell types. GeneScores are simply the summation of all the peak fragments that are present within gene bodies and their corresponding promoter region. Since, marker genes are often better understood than the non-coding regions the GeneScores can be used to annotate cell types.

In Scarf, the users can provide a BED file containing gene annotations. This BED should have no header, should be tab separated and must have the columns in following order:

  1. chromosome identifier

  2. start coordinate

  3. end coordinate

  4. Gene ID

  5. Gene Name

  6. Strand (Optional)

The start/end coordinate can extend through transcription start site (TSS) to include a portion of promoter.

For convenience we have generated such BED files for human and mouse assemblies using the annotation information from GENCODE project. We downloaded the GFF3 format primary chromosome annotations and used Scarf’s GFFReader to convert the files into BED and add promoter offset of 2KB. These BED files containing gene annotations can be downloaded using fetch_dataset command and passing ‘annotations’ parameter.

scarf.fetch_dataset(
    dataset_name='annotations', 
    save_path='scarf_datasets'
)

time: 15.5 s (started: 2023-05-24 14:19:53 +00:00)

Now we have the annotations and are ready to calculate the ‘GeneScores’. The add_melded_assay is the DataStore method that will be used for this purpose. The add_melded_assay method is actually designed to map any arbitrary genomic coordinate information (not just gene annotations) to the ATAC peaks. Here, we use the term ‘melding’ for the process wherein for a given loci the values from all the overlapping peaks are merged/melded into that feature. The peaks values are TF-IDF normalized before melding.

Now we explain the parameters that are usually passed to the add_melded_assay.

  • from_assay: Name of the assay to be acted on. You can generally skip if you have only one assay. We use this parameter here only for demonstration puspose

  • external_bed_fn: This is the annotation file. Here we pass the annotation for human GRCh37/hg19 assembly based GENCODE v38 annotations.

  • peak_col: This is the column in ds.ATAC.feats table that contains the peak coordinates. In this case it is ‘ids’, but could be anyt other column. Please note that the coordinate information in this column should be in format ‘chr:start-end’ (please note the positions of colon and hyphen)

  • renormalization: By overriding the default value of True to False here, we turned off ‘renormalization’ step. The renormalization step makes sure that all the feature values for each cells in the melded assay sum up to the same value. Here we turned this off because we will have ‘GeneScores’ as an ‘RNAassay’ which uses library size normalization that has the same effect as ‘renormalization’.

  • assay_label: The label/name of the melded/output assay. Because we are using gene bodies as our input, ‘GeneScores’ is sensible choice.

  • assay_type: Here we set the type of assay as ‘RNA’ which means that the new melded assay will be treated as if it was an scRNA-Seq assay. Alternatively, we could have set it as a generic ‘Assay’

ds.add_melded_assay(
    from_assay='ATAC',
    external_bed_fn='scarf_datasets/annotations/human_GRCh37_gencode_v38_gene_body.bed.gz',
    peaks_col='ids',
    renormalization=False,
    assay_label='GeneScores',
    assay_type='RNA'
)

INFO: 33898/62409 features did not overlap with any peak

WARNING: Minimum cell count (6.539499236477921) is lower than size factor multiplier (1000)

time: 5min 30s (started: 2023-05-24 14:20:09 +00:00)

We can now print out the DataStore and see that the ‘GeneScores’ assay has indeed been added. The ‘add_melded_assay’ also printed a useful bit of information that almost half of features(genes bodies) did not overlap with even a single peak. The melded assay anyway contains all the genes but sets these ‘empty’ genes as invalid.

ds

DataStore has 8426 (8728) cells with 2 assays: ATAC GeneScores
   Cell metadata:
            'I', 'ids', 'names', 'ATAC_UMAP1', 'ATAC_UMAP2', 
            'ATAC_leiden_cluster', 'ATAC_nCounts', 'ATAC_nFeatures', 'GeneScores_nCounts', 'GeneScores_nFeatures', 
            'GeneScores_percentMito', 'GeneScores_percentRibo'
   ATAC assay has 86049 (89796) features and following metadata:
            'I', 'ids', 'names', 'I__prevalent_peaks', 'dropOuts', 
            'nCells'
   GeneScores assay has 28357 (62409) features and following metadata:
            'I', 'ids', 'names', 'dropOuts', 'nCells', 
          

time: 17.4 ms (started: 2023-05-24 14:25:39 +00:00)

Let’s now visualize some ‘GeneScores’ for some of the known marker genes for PBMCs on the UMAP plot calculated on the ATAC assay

ds.plot_layout(
    layout_key='ATAC_UMAP', from_assay='GeneScores', 
    color_by=['CD3D', 'MS4A1', 'LEF1', 'NKG7', 'TREM1', 'LYZ'],
    clip_fraction=0.01, 
    n_columns=3,
    width=3,
    height=3,
    point_size=5,
    scatter_kwargs={'lw': 0.01},
)
../_images/c2c052d48aac0693db2660dac0465f5909eaa8ed31a193d0f623f44b76d48e10.png 
time: 13 s (started: 2023-05-24 14:25:39 +00:00)

We stop here in this vignette. But will soon add other vignettes that show how we can other kinds of melded assays like motifs, enhancers, etc. We will also explore how we can use the ‘GeneScores’ to integrate scATAC-Seq datasets with scRNA-Seq datasets from same population of cells (but not the exact same set of cells)

That is all for this vignette.