Translator: A Transfer Learning Approach to Facilitate Single-Cell ATAC-Seq Data Analysis from Reference Dataset

2.1.1. scATAC-seq data simulation with ground truth labels

Due to the challenge in obtaining gold-standard datasets, we first used a public software SCAN-ATAC-sim to generate labeled scATAC-seq data by efficiently down sampling bulk ATAC-seq profiles from ENCODE cell lines (Chen et al., 2021c). By tuning two key parameters: SNR and average fragment number per cell ($\mu$), as suggested by the SCAN-ATAC-sim website, we simulated two kinds of scATAC-seq data to mimic the HQ reference and the low-quality (LQ) target data scenarios. We repeated this procedure for five different cell types, MONO, NEU, CMP, MEGA, and ERY, from the ENCODE data portal (Consortium, 2012; Davis et al., 2018) and generated the signal tracks using the SCAN-ATAC-sim preprocessor.

To test Translator's robustness, we evaluated four possibilities: fixed depth, varying depth, mismatched cells with new cell types, and mismatched cells with missing cell types (Table 1). In addition, to evaluate the performance of eliminating confounding factors of Translator's invariant learning, we applied two experiments: one using one sample ($\mu =5000,\sigma =1.5,\rho =0.4$) to test the effectiveness of removing sequencing depths; another using two samples (${\mu }_1=2500,{\sigma }_1=1.5,{\rho }_1=0.4,{\mu }_2=5000,{\sigma }_2=1.5,{\rho }_2=0.5$) to test the performance under two samples.

2.1.2. PBMC scATAC-seq data preprocessing

To test Translator's performance on real single-cell datasets, we downloaded the fragment files of the PBMC multiome dataset publicly available from 10x Genomics. We used the ATAC modality for performance benchmarking and the RNA modality for more accurate labeling due to its higher quality. To mimic the scenario of datasets with different qualities, we took the PBMC multiome 10k dataset as the reference and the down-sampled PBMC multiome 3k dataset as the target dataset. The sequencing depth and other quality control (QC) parameters are provided in Supplementary Figure S1 and S2.

Specifically, we down-sampled the target dataset using the following algorithm. Given the sequencing depth of the raw target dataset d_r and the desired sequencing depth $d<d_r$, we calculated $r=\frac{\left(d_r-d\right)}{d_r}$, which is the proportion of the peaks needed to be deactivated. To avoid samples with extremely low sequencing depth, we set 1000 fragments as the lower bound depth. For each cell i, where n_i is the number of open chromatin peaks, we randomly sampled $min\left(r{n}_i,n_i-1000\right)$ peaks to be deactivated by setting their corresponding positions in the peak matrix as zero. Further iterations were conducted until the sequencing depth reached the desired sequencing depth $d_r\pm 50$.

We then used ArchR (Cao et al., 2021; Granja et al., 2021) (version 1.0.1) to preprocess the fragment file of the reference data using the default parameters. Specifically, we created an ArchR object and removed barcodes with a transcription starting site (TSS) enrichment score less than 4 or whose number of fragments was less than 1000 (as suggested by ArchR Tutorial). We further removed doublets using ArchR's default filters, kept peaks within the autosome chromosomes, and then generated the reference data cell-by-peak matrix. Finally, we binarized the reference cell-by-peak matrix to set any value greater than 1 as 1. For the target data, Translator imported the reference peak files to create the target cell-by-peak matrix for the training process. To create benchmarking metrics for the target dataset, we used the latent semantic indexing (LSI) commonly utilized in natural language processing models as well as in scATAC-seq analyses (Cao et al., 2021; Stuart et al., 2021). Given the binarized target data, we utilized the RunLSI function from the R package Signac (Hao et al., 2021) (version 4.0.4) with default parameters and n = 20 to keep its latent dimension the same as Translator's.

2.1.3. Preprocessing scATAC-seq data in the human prefrontal cortex

In addition to blood cells, we also looked at the prefrontal cortex (PFC) to further evaluate our model (post-mortem tissues were used; no IRB applicable). Previously, we deeply sequenced scATAC-seq data from three frozen PFC tissues MS0169WW, MS0177EE, and MS0181II. We conducted the following protocol to generate scATAC-seq data from the prefrontal cortex tissues: single-nucleus suspensions were isolated from 25 mg of frozen human dorsolateral prefrontal cortex (Brodmann Area 9/46). Tissue was initially lysed and nuclei released by using a sucrose-based solution and NP-40 detergent (Sigma) in Dounce homogenizer followed by centrifugation at 1000 g for 10 minutes.

The nuclear pellet was resuspended and further purified using an Iodixanol (Optiprep, AxisShield) solution gradient by an additional centrifugation (3000 g for 30 minutes) and collection of nuclei in the interphase (30%/40%). Nuclei were washed and resuspended before running on the 10X Chromium Single Cell ATAC platform (#PN-1000110; 10x Genomics). Library quantification, quality checking, and sequencing (250 million reads) were done according to the manufacturer's recommendations by using the Illumina NovaSeq6000 (Illumina) flow cell S4 (Illumina) at the Yale Center for Genome Analysis.

Following the data generation, we used cell-ranger ATAC (version 2.0.0) to preprocess the raw reads with default parameters (Satpathy et al., 2019). Then, for quality control, cells with TSS enrichment score less than 4 or number of fragments less than 1000 were filtered out using ArchR. In addition, Harmony was used to integrate samples (Korsunsky et al., 2019). The sequencing depth and other QC parameters are provided in Supplementary Figures S3–S5. We performed dimensionality reduction, clustering, and cell type identification using default parameters in ArchR. Cells from these three samples were clearly clustered into biologically relevant cell types that matched known markers for neuronal and non-neuronal cell types, validating their high sequencing quality.

Additionally, a lower-quality PFC scATAC data was also collected in the same cohort, with moderate sequencing depth ($2743.5$), lower cell number (6285), and lower TSS enrichment than that of the reference dataset ($4.35$). Supplementary Figure S6 indicates the sample's moderate quality. Similar to the reference dataset, default parameters were used in cell-ranger (version 6.0.1) and ArchR to preprocess the lower-quality data.

For all PFC data, after the QC step, we added the reference peak set, which was the union of all three references' cell type-specific peak sets sorted by chromosomes. We generated the count matrix using the addPeakMatrix function from ArchR. Finally, we saved the matrix using getMatrixFromProject with useMatrix set to PeakMatrix and binarize set to True so that the exported matrix was binarized. Using the writeMM function in R allowed us to save the exported sparse matrix as a mtx file and then exported as a npz file usable in Python.

2.1.4. Data preprocessing of the Share-seq mouse skin scATAC-seq data

In addition, we further applied Translator on the mouse skin dataset from Share-seq dataset (Ma et al., 2020). We extracted the count matrix from its raw data using the default peak set provided with Signac. To construct the reference and the target data, we randomly sampled 1000 cells as the target data. The rest of them would be the reference data. We then conducted a feature selection workflow similar to SCALE (Xiong et al., 2019) by filtering out all peaks with a total signal less than 500 across the entire dataset. Then, we chose the top 80,000 peaks as inputs with the TF-IDF algorithm.

2.2.1. Deep VAE with invariant representation learning

To capture complex feature interactions and remove various confounding factors (e.g., age, gender, sequencing platforms, and batch), Translator first used a deep VAE with an invariant representation learning scheme to learn robust cell embeddings free from bias (Kingma and Welling, 2013). Specifically, we encouraged cell representations $z\in {ℛ}^d$ to only reflect intrinsic biological states but be independent of confounding factors c. Specifically, a penalty term was added in the VAE to minimize the mutual information $I\left(z,c\right)$. Given a peak profile $x\in 0,1^n,$ we used a VAE to handle complex peak interactions and maximize evidence lower bound of the log-likelihood

In Equation (1), the encoder probability, denoted by $q_{\theta }\left(z|x\right)$, is the probability of the latent embedding z after the reparameterization trick given the input data x and the encoder parameter $\theta$. The decoder probability $p_{\varphi }\left(x|z,c\right)$ is the likelihood of the reconstructed input data x given the latent embedding z and other confounding factors c (see Table 3 for detailed parameterizations). Since z could potentially depend on c, we added an upper bounded penalty $I\left(z,c\right)$ to encourage the independence of z and c, where

where the $H\left(x|c\right)$ is a constant and can be removed from Equation (2). Additionally, in VAE training, a common challenge was Kullback-Leibler (KL) annealing, lowering $D_{K}L\left(\cdot \right)$ with respect to training epochs. Therefore, we applied the cyclic annealing schedule (Fu et al., 2019) by adding a term $w=min\left(1,\frac{2\left(nmodN\right)}{N}\right)$, where n is the current epoch and N is the number of epochs in a cycle. Therefore, the final loss function we aimed to minimize was

2.2.2. Neural network architecture

A typical scATAC-seq data will generate at least 20k peaks as the input, so the first fully connected layer in the VAE model contributes the largest to the number of parameters used and usually requires a large number of cells for a stable training process. To reduce over-parameterization on data with a moderate cell number, we only allowed intrachromosomal peak-neuron connections in the first layer of the VAE with a straightforward intuition: CRE interactions usually happen within a single chromosome instead of across different chromosomes. As shown in Figure 1, each chromosome had its own fully connected network with 64–32 units. This way, the number of parameters in the first two neural network layers could be dramatically reduced to less than $10\%$ of those in the fully connected layer, improving training efficiency and reducing the number of parameters trained to avoid potential overfitting. The 32-unit layers across all chromosomes were then concatenated, followed by a fully connected layer of 20 units, indicating the multivariate Gaussian mean and variance. The decoder was symmetric to the encoder.

2.2.3. Training procedure

Given a set of single-cell ATAC-seq data, we first sorted and grouped the peaks by their chromosome. Then, we fed the vertically stacked reference and target data into the encoder together in batches to obtain the learned mean and variance. The reparameterization step took the multivariate mean and variance from the encoder, generated a random sample, and fed them into the decoder. After the loss was computed, we conducted back propagation to update the parameters. Adam optimizer was used throughout the process. The main hyperparameters included learning rate lr, weight decay $\alpha$, and the $\lambda$ term in the loss function. To tune these parameters, for each training task, only the warmup procedure was done using parameters $lr\in \left\{1E-4,5E-4,1E-3,2E-3\right\},\alpha \in \left\{1E-5,1E-4,5E-4\right\},\lambda \in \left\{0.1,0.2,0.5,1,2,5,10\right\}.$ The set of hyperparameters that led to the best loss after the warmup was chosen to train the final model. All trainings were conducted using NVIDIA Tesla K80 GPUs.

2.2.4. Training from scratch using VAE

To evaluate the performance of Translator with and without the reference dataset, we adopted a baseline training scheme without the reference data, namely VAE-SCRATCH. In this protocol, only the target data were used in the training process. The training procedure was identical to the Translator training procedure. After training, we directly generated the latent embedding using the trained encoder with the target data.

2.3.1. Motif analysis for PBMC dataset

To better understand the cells and different clusters of the PBMC dataset, we generated the motifs for all cell types using Translator's transferred latent embeddings. Since we borrowed the cell type labels from the RNA modality, we could compare and contrast generated motifs and potentially find new subtypes that enables further downstream analysis. Using the R package Signac (Stuart et al., 2020) (version 1.3.0), we first created a ChromatinAssay object, and then used the getMatrixSet function to get a set of motifs from the JASPAR2020 dataset (Fornes et al., 2020) with collection being CORE and tax group being vertebrates. Then, we used the AddMotifs function to add both the reference genome, UCSC hg38 (Schneider et al., 2017), and the motif to the object. Using the function runChromVAR in signac, we received a matrix of motif enrichment scores for each motif in the set and for each cell (Schep et al., 2017). Using the LEIDEN clusters calculated based on latent embeddings, we plotted and generated statistical metrics of the distribution of certain cell type-specific motifs for each cluster, allowing qualitative and quantitative comparisons.

2.3.2. Gene enrichment analysis and cell type annotation for the postmortem brain dataset

We also conducted gene enrichment analysis for the PFC scATAC-seq dataset to evaluate the clustering accuracy of our model. In the analysis, we clustered Translator's latent embedding with LEIDEN. Then, using the list of cell-type marker genes identified in previous work (Lake et al., 2016; Zonouski et al., 2019), we curated a gene score matrix with a score for each cell and each gene. More specifically, we used the following marker genes to determine the cell-type annotations—astrocytes: GLUL, SOX9, AQP4, GJA1, NDRG2, and GFAP; microglia: MRC1, TMEM119, and CX3CR1; endothelial: CLDN5; OPCs: PDGFRA, PCDH15, and OLIG2; oligodendrocytes: PLP1, MAG, and MOG; excitatory: SATB2 and SLC17A7; and inhibitory: GAD1 and SLC32A1. Finally, for each gene in every cluster, we inspected the gene enrichment score distribution to assign cell types. As a result, qualitative evaluation was done using overlapped UMAP visualization between the reference and target datasets. In addition, to evaluate the performance of Translator under similar or distinct data, we train the PFC target data with PBMC reference and conduct a UMAP and ARI analysis as well as comparison.

2.3.3. Reference and target analysis of the Share-seq dataset

We conducted multiple experiments to evaluate our model on Share-seq. First, we train Translator, Signac, ArchR, and SCALE (Xiong et al., 2019; Cao et al., 2021; Stuart et al., 2021) on the same preprocessed data as aforementioned, all with 20 dimensions. Note that because we have done feature selection, we disabled the feature selection function on SCALE. As a quantitative measurement, UMAP was generated using the same methods for both reference and target data. Quantitative measurements include silhouette scores and ARI, calculated consistent with the other datasets.