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
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
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.
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
...
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.
