Tn5bias

To accurately eliminate the insertion bias of the Tn5 transposase enzyme, it is important to have a good bias model specific to the data and condition being explored. It is recommended to generate a separate bias model for each of the conditions being explored. This meant that within the scope of this research, the following models were generated:

• Tn5 bias model for mm10 adult scATAC-seq
• Tn5 bias model for mm10 neonatal scATAC-seq
• One Tn5 bias for each bulk cell line investigated (e.g., rn6 GH3, TaT1, At-T20, Lbt2)

All bias models generated within this research are available on HuggingFace.

Because the parameters are highly dependent on the condition and cell line of the model being explored, an iterative function was not created, and the chrombpnet bias pipeline was the function used, as described in the original ChromBPNet repository wiki. It is recommended to also look through their wiki page as well as this one when generating and evaluating a model.

First make a directory for the Tn5 bias models:

mkdir ~/ChromBPNet/Bias_Models

Model Training

chrombpnet bias pipeline \
        -ifrag ~/atac_seq_data/mouse_merged/all_cells_merged_sorted.tsv.gz \
        -d "ATAC" \
        -g ~/genome/mm10/mm10.fa \
        -c ~/genome/mm10/mm10.chrom.sizes \
        -p ~/ChromBPNet/data/mm10/peaks/peaks_no_blacklist.bed \
        -n ~/ChromBPNet/data/mm10/peaks/output_fold_0_negatives.bed \
        -fl ~/ChromBPNet/splits/mm10/fold_0.json \
        -b 0.5 \
        -o ~/ChromBPNet/Bias_Models/Tn5_bias_model_fold_0_mm10_adult_0.5

Input Format

As taken from the ChromBPNet repository

-ifrag     # Fragment file of the merged cohort (-ibam if using a BAM file instead)
-d         # Data type 
-g         # Genome fasta file
-c         # Chromosome sizes (see Note 1 below)
-p         # Peaks file with no overlap with the blacklist (generated in PreprocessPeaks)
-n         # Non-peaks (generated in PreprocessPeaks)
-fl        # Fold split used - fold 0 used continuously in this research (generated in PreprocessPeaks)
-b         # Bias threshold factor (see Note 2 below)
-o         # Output path - does not have to exist prior to running command

Note 1: The chromosome size file, when downloaded, will not include any scaffold regions within the genome. If the fragment file includes these regions, the function will return an error and not run. It is therefore important to check first, and if the regions are included they must first be filtered out.

cut -f1 ~/genome/mm10/mm10.chrom.sizes > valid_chroms.txt

zcat all_cells_merged_sorted.tsv.gz | \
  awk -v OFS="\t" 'NR==FNR{valid[$1]=1; next} ($1 in valid)' valid_chroms.txt - | \
  gzip > all_cells_merged_sorted.filtered.tsv.gz

Note 2: The bias threshold value is calculated as a threshold applied on maximum count of non-peak region for training bias model, which is set as this threshold x min(count over peak regions). The recommended start value for ATAC-seq data is 0.5.

The large-scale cohorts of this study used this threshold set at 0.5, but for the bulk models it required an increase to 0.9, because the peak calling was less stringent, potentially including non-peak regions and therefore requiring a larger value.

It is important to run the bias models iteratively and assess the outcomes on a case-by-case basis until an accurate model has been generated.

Output Format

As taken from the ChromBPNet repository

models\
	bias.h5    # Final bias model
logs\
	bias.log (loss per epoch)
	bias.log.batch (loss per batch per epoch)
	(..other hyperparameters used in training)
	
intermediates\
	...

evaluation\
        overall_report.html      # Summary of the training and motifs learnt by the bias model with quality check guidance.
        overall_report.pdf       # Same summary in PDF format
	pwm_from_input.png
        k562_epoch_loss.png 
	bias_metrics.json
	bias_only_peaks.counts_pearsonr.png
	bias_only_peaks.profile_jsd.png
	bias_only_nonpeaks.counts_pearsonr.png
	bias_only_nonpeaks.profile_jsd.png
        bias_predictions.h5
	bias_profile.pdf
	bias_counts.pdf
	...

Evaluating the models

Following training, DeepLIFT is automatically applied by the ChromBPNet pipeline to compute per-nucleotide sequence contribution scores, and TF-MoDISco used to derive sequence motifs from these scores. Bias model outputs were evaluated iteratively. If biologically interpretable TF motifs were detected in either the profile or counts contribution score outputs, the bias_threshold_factor was reduced and the model retrained, as such motifs indicate incomplete exclusion of TF binding signal. Iterative refinement was continued until all models yielded Tn5 insertion sequence signatures consistent with known GC-enriched Tn5 bias.

Total signal performance was quantified as the Pearson correlation coefficient (r) between predicted and observed log-transformed read counts across all evaluated genomic regions (Fig. A). For non-peak regions, high correlation is expected, as the bias model should accurately capture the Tn5 insertion preference (Fig. 1B). For peak regions, a low correlation is desired, indicating that the model does not model biological accessibility signals, only captures Tn5 bias (Fig. 1C). However, values below -0.3 risk capturing AT-content bias and were considered unsuitable for downstream ChromBPNet training, with bias models yielding correlations below -0.5 being incompatible with the ChromBPNet pipeline.

Profile shape performance was assessed using the Jensen-Shannon Divergence (JSD) between observed and predicted per-position accessibility profiles (Fig. 1A). Lower JSD values reflect greater agreement between predicted and observed profiles. However, as JSD is sensitive to sequencing depth, it was used as a qualitative rather than threshold-based metric.

All trained bias models returned Tn5 insertion sequence motifs characterised by GC-enriched sequence content (Fig. 1D), consistent with the established biochemical preference of the Tn5 transposase and confirming appropriate bias model behaviour.