Input files (before launching ChromSCape)

Various existing technologies allow to produce single-cell genome-wide epigenomic datasets :

• scChIP-seq (Grosselin et al., 2019) scCUT&Tag (Kaya-Okur et al., 2019), scChIL-seq (Harada et al., 2019), scChIC-seq (Ku et al.,) reveal the distribution of histone marks (H3K27me3, H3K4me3) or transcription factors (RNA Polymerase 2,…) with single-cell resolution.
• scATAC-seq (Buenorostro et al., 2015) or sciATAC-seq (Cusanovitch et al.,

2016. assess regions of open chromatin in single-cells

• scDNA methylation profiling, scRRBS, scWGBS …
  

After sequencing a single-cell epigenomic experiment, the raw sequencing files (.fastq) are demultiplexed, aligned against a reference genome to output different kind of data depending on the technology. ChromSCape allows user to input a variety of different format. Depending on the output of the data-engineering/pre-processing pipeline used, the signal can be already summarized into features (Count matrix, Peak-Index-Barcode files) or stored directly in raw format (single-cell BAM or single-cell BED files).

Anyhow the format, ChromSCape needs signal to be summarized into features. If inputting raw signal (scBAM or scBED), the application lets user summarize signal of each cells into various features:

• Genomic features (extended region around TSS of genes, enhancers)
• Peaks called on bulk or single-cell signal (BED file must be provided by the user)
• Genomic bins (windows of constant length, e.g. 100kbp, 50kbp, 5kbp…)

Note that summarizing will take longer if using BAM files than BED files, and if the number of features is important (e.g. 5kbp bins, enhancers…).

Create & import a dataset

Once you have launched the application, you arrive on the “Select & Import”. Here you can selet an output directory, where your analyses will be saved. ChromSCape will automatically save the folder’s location so that you don’t have to select it each time you connect.

You must then name your analysis. The name shouldn’t contain special characters (except ’_‘). Choose either the Human (hg38) or Mouse (mm10) genome. This is used to annotate your features with the closest genes TSS for the Gene Set Analysis. Browse your computer for one or multiple matrices that will be analyzed together. To select multiple matrices, they must be placed in the same folder and then the user can select multiples matrices with (Shift + Click) or (Ctrl + Click). Finally, clicking on ’Create Analysis’ will create an analysis & import the count matrices in this analysis. This will create a folder named ‘ChromSCape_analyses’ in the output directory you specified, inside which another folder ‘Your_analysis_name’ is nested. If you create another analysis, it will also be created under ‘ChromSCape_analyses’. If you already created an analysis in a previous section, selecting the same output directory as you chose in the previous session will enable you to load your analysis.

Select & Import Tab: Select an output directory, Create an analysis or load previous analyses, import count matrix into your new analysis.

Select an type of format input, and eventually a type of feature to summarize signal.

2. Filter, Normalize & Reduce Dimensionality

In order to efficiently remove outlier cells from the analysis, e.g. cells with excessively high or low coverage, the user sets a threshold on a minimum read count per cell and the upper percentile of cells to remove. The latter could correspond to doublets, e.g. two cells in one droplet, while lowly covered cells are not informative enough or may correspond to barcodes ligated to contaminant DNA or library artifacts. Regions not supported by a minimum user-defined percentage of cells that have a coverage greater than 1,000 reads are filtered out.

Defaults parameters were chosen based on the analysis of multiple scChIP-seq datasets:

• A minimum coverage of 1,600 unique reads per cell
• Filtering out the cells with the top 5% coverage
• Keeping regions detected in at least 1 % of cells.

Post quality control filtering, the matrices are normalized by total read count and region size (CPM normalization. At this step, the user can provide a list of genomic regions, in BED format, to exclude from the subsequent analysis, in cases of known copy number variation regions between cells for example.

To reduce the dimensions of the normalized matrix for further analysis, PCA is applied to the matrix, with centering, and the 50 first PCs are kept for further analysis.

If you want to run t-SNE for vizualisation in addition to PC and UMAP, tick the ‘Run T-SNE’ checkbox. This might take some time to run, depending on the size of the dataset.

For some large datasets (e.g. more than 15,000 cells) or unbalance datasets (i.e. number of cells per sample is very different), it might be a good option to first run with the subsampling option. To do so, tick the ‘Perform subsampling’ and select the maximum number of cells to subsample from each sample. The sampling is done without replacement.

Users might want to remove specific loci from the analysis, e.g. repeat regions, known Copy Number Altered regions that might bias the analysis. To do so, tick the ‘Exclude specific genomic regions’ and browse for a BED file containing the regions to exclude.

A batch correction option using mutual nearest neighbors ‘FastMNN’ function from batchelor is implemented to remove any known batch effect in the reduced feature space. To do so, tick the ‘Perform batch correction’, select the number of batches and add each sample to a batch.

Filter, Normalize & Reduce Tab: Adjust paramters to filter out low covered cells & features, as well as too highly covered cells that might be potential doublets. Normalize & Reduce dimensionality using PCA.

Visualize cells in reduced dimensions

The user can visualize the data after quality control in the PCA , t-SNE or UMAP reduced dimensional space. The PCA and t-SNE plots are a convenient way to check if cells form clusters in a way that was expected before any clustering method is applied. For instance, the user should verify whether the QC filtering steps and normalization procedure were efficient by checking the distribution of cells in PC1 and PC2 space. Cells should group independently of normalized coverage.

Depending on wether you have ticked the ‘Run T-SNE’ checkbox in the previous Tab, you will see the T-SNE two dimensional plot.

Note: You can change the colors of each samples, batch or cluster and the saved colors will be saved in the Analysis.

Visualize Cells: Visualize cells in reduced dimensions: PCA, UMAP & TSNE. Color cells by sample, batch, or total count to make sure cells are not separating due to library size only.

A brief understanding into what features drives the main PCs is shown in the ‘Contribution to PCA’ box.

• In the barplot, you will get the top loci / closest genes associated to the top loci contributing to the Component of your choice. The top contributing features are taken by sorting the values of the features in the componenent of intereset (e.g. PC1, PC2…). If you observe a huge gap between the top 2-3 features and the rest, this means that very few features are driving most of the dimensionality reduction, and that other less contributing features might be hidden. This also gives you interesting features to look at in the Peak Calling & Coverage or Differential Analysis tabs.

• In the pie chart, the chromosome repartition of the top 100 most contributing features to the componenent of interest (e.g. PC1, PC2…) are displayed. This is used to see if there is any strong imbalance in the “chromosome contribution”, e.g. if a chromosome is contributing more than it should be. Usually, if a single chromosome has > 30% representation, this might indicate that CNVs or genetic events are happening on this chromosome, and you might want to look into more details at these events, or exclude them from the analysis (see Filter & Normalize tab)

Most contributing features

Cluster cells

Using the first 50 first PCs of computed PCA as input, hierarchical clustering is performed, taking 1-Pearson’s correlation score as distance metric. The cell to cell correlation heatmap and dendogram describing the hierarchical clustering is represented in the left panel.

To improve the stability of the clustering approaches and to remove from the analysis isolated cells that do not belong to any subgroup, cells displaying a Pearson’s pairwise correlation score below a threshold t with at least p % of cells can be filtered out. The correlation threshold t is calculated as a user-defined percentile of Pearson’s pairwise correlation scores for a randomized dataset (percentile is recommended to be set as the 99th percentile). To do so, open the “Filter lowly correlated cells” panel and click on “Filter & Save”. This will filter any cells that are too rare to correlate with enough other cells and form stable clusters.

To select a number of clusters, looking at the correlation heatmap and dendogram, choose a number in the ‘Select Cluster Number’ panel and click on ‘Choose Cluster’.

If you are unsure of how many cluster to choose, you can run consensus clustering using Bioconductor’s package ConsensusClusterPlus package to determine what is the appropriate k-partition of the dataset into k clusters. To do so, it evaluates the stability of the clusters and computes item consensus score for each cell for each possible partition from k=2 to n. For each k, consensus partitions of the dataset are done on the basis of 1,000 resampling iterations (80% of cells sampled at each iteration) of hierarchical clustering, with Pearson’s dissimilarity as the distance metric and Ward’s method for linkage analysis. The optimal number of clusters is then chosen by the user; one option is to maximize intra-cluster correlation scores based on the graphics displayed on the ‘Consensus Clustering’ Tab after processing.

If you ran consensus clustering, re-pick the cluster numbers in the ‘Select Cluster Number’ panel making sure you check the ‘Use consensus’ checkbox and re-click on ‘Choose Cluster’.

Cell clusters can then be visualized in two dimensions with the t-SNE or UMAP plot.

Note: The correlation filtering and consensus clustering steps are optional as they need long runtime. You might want to proceed to Differential Analysis and Gene Set Analysis directly if you are confident in your choice of clusters.

Cluster Cells: Choose a number of cluster if you have clear clusters of cells, or run consensus clustering to test multiple clustering and find the optimal cluster number. Optionally filter out lowly correlated cells that might impact the clustering performance.

## Inter & Intra correlation violin plots

In this box, the heterogeneity of cells in various groups (e.g. sample or cluster) are displayed with two plots:

• In the intra-correlation plot, the distribution of cell-to-cell Pearson correlation scores within each cluster or sample is displayed. Heterogeneous groups will have lower correlation than homogeneous groups.

• In the inter-correlation plot, the distribution of cell-to-cell Pearson correlation scores between each cluster or sample is displayed. This allows to understand which groups are closer to each other than other groups.

For both plots, you can add a layer of information by adding the single-cell as jitter points on the plots. This is interesting when combining samples and clusters to see the repartition of conditions into the clusters.

Inter/Intra correlation

Peak Calling to refine bins to gene annotation

This step of the analysis is optional, but recommended for Gene Set Analysis if working on genomic bins in order to refine the peak annotation prior to enrichment analysis. To be able to run this module, some additional command line tools are required such as SAMtools, BEDTools and MACS2. Therefore, this step is only available for Linux and MacOS systems.

You need to input BAM files for the samples (one separate BAM file per sample), with each read being labeled with the barcode ID in the ‘XB’ tag. ChromSCape merges all files and splits them again according to the previously determined clusters of cells (one separate BAM file per cluster). To do so, place all your BAM files in a folder and click on the ‘Choose a folder containing the BAM files’ button to browse for your folder. If some of the BAM files in the folder are not part of your analyis, remove them by pressing the write box and deleting them.

You can select significance threshold for peak detection (passed to MACS2) and merging distance for peaks (defaults to p-value=0.05 and peak merge distance to 5,000). This allows to identify peaks in close proximity (<1,000bp) to a gene transcription start site (TSS); these genes will be later used as input for the enrichment analysis.

Peak Calling to refine bins to gene annotation: Call peaks on ‘pseudo-bulk’ cell clusters, and refine bins to genes annotation.

Differential Analysis

To identify differentially enriched regions across single-cells for a given cluster, ChromSCape can perform 1. A non-parametric two-sided Wilcoxon rank sum test comparing normalized counts from individual cells from one cluster versus all other cells, or cluster of choice, 2. Or a parametric test comparing raw counts from individual cells, using the Bioconductor package edgeR, based on the assumption that the data follows a negative-binomial distribution. We test for the null hypothesis that the distribution of normalized counts from the two compared groups have the same median, with a confidence interval 0.95. The calculated p-values are then corrected by the Benjamini-Hocheberg procedure. You can modify the log2 fold-change threshold and corrected p-value threshold for regions to be considered as significantly differentially enriched (default settings are a p-value and log2 fold-change thresholds respectively of 0.01 and 1).

If you have specified batches to correct for batch effect, the differential analysis is done using the ‘pairwiseWilcox’ function from the Bioconductor package scran , setting the batch of origin as a ‘blocking level’ for each cell.

Differential Analysis: Identify most differential features between the various cell clusters

Gene Sets Analysis

In this tab, using the refined annotation of peaks if available, the final step is to look for enriched gene sets of the MSigDB v5 database within differentially enriched regions (either enriched or depleted regions in the studied histone mark). Click on ‘Start Enrichment Analysis’ in order to apply hypergeometric tests to identify gene sets from the MSigDB v5 database over-represented within differentially enriched regions, correcting for multiple testing with the Benjamini-Hochberg procedure. Users can then visualize most significantly enriched or depleted gene sets corresponding to the epigenetic signatures of each cluster and download gene sets enrichment tables.

The top 100 most significant differential regions, single-cell Epigenomic enrichment levels can be visualized overlaying read counts for each cell at selected genes onto a t-SNE plot (right panel).

Gene Set Analysis: Look for pathways activated or depleted in most differential regions between clusters