TeloTales: Telomere Content Analysis Pipeline

Reproducible bioinformatics pipeline for estimating telomere content from whole-genome sequencing (WGS) data from the 1000 Genomes Project.

This repository presents a Snakemake pipeline for automated retrieval of BAM/BAI files, extraction of telomeric repeat signals using qmotif, and downstream analysis of telomere content across samples.

Data are accessed directly from the 1000 Genomes Project via NCBI, and are not stored within this repository.

Table of Contents

• Requirements
• How the Pipeline Works
• Performance
• File Descriptions
• Repository Structure
• Setup
• Usage
• Outputs
• Citation
• Contact

Requirements

To run this pipeline, you'll need:

• Snakemake: to run the pipeline

• Java: Required to run qmotif

• qmotif: Download and install qmotif (Documentation here)

• A note on qmotif scaling: qmotif normalizes raw motif counts against a BAM file's total read count, scaled to 1 billion reads, so that samples with different sequencing depths can be compared on the same scale. This scaling treats every read equally and doesn't account for factors like unmapped reads or copy number changes (e.g., whole-arm amplifications), since there's no single "correct" way to adjust for those in tumor samples. Despite this simplicity, qmotif's scaled scores have been shown to correlate well with wet-lab telomere length measurements. See the qmotif documentation for more detail.

Supplementary Figure 1. Workflow for Telomere Content Variation Pipeline. Schematic representation of the qmotif-based pipeline used to estimate telomere content variation across samples from Phase 3 of the 1000 Genomes Project. The workflow involves extracting and quantifying telomeric reads from whole-genome sequencing data using qmotif v1.0. The tool operates through a two-stage matching system: in Stage 1, a simple string match is used to identify canonical telomeric repeats, while in Stage 2, a more complex regular expression is applied to detect variant telomeric sequences. At the end of Stage 2, a tally of all identified motifs is done, and the final number is recorded. [Figure created with BioRender.com].

Additional files to prepare:

• files_list.txt: contains relative paths to BAM and corresponding BAI (BAM index) files from the 1000 Genomes Project, one per line (the base URL they're downloaded from is set separately in config.yaml)
• chrnames.txt: contains chromosome names for tallying telomeric reads (autosomes only — chr1–chr22; sex chromosomes are not included)
• qmotif.ini: qmotif configuration file used for all analyses in the paper (included in workflow/scripts/)

How the Pipeline Works

For each sample, the Snakefile runs these steps:

1. Download: fetches the BAM/BAI listed in files_list.txt from the 1000 Genomes FTP.
2. qmotif: scans every read using qmotif.ini:
   • Stage 1: exact match to TTAGGGTTAGGGTTAGGG (a fast string prefilter for canonical repeats)
   • Stage 2: regex (...GGG){2,}|(CCC...){2,} (captures variant repeats on both strands)
   • revcomp=true, window_size=10000, cutoff_size=5
   • INCLUDES regions: the terminal windows of each chromosome end, defined on hg19 (see [INCLUDES] in qmotif.ini)
3. Per-chromosome tally: counts stage 2 telomeric reads aligning to the terminal window of each autosome (chr1–22, listed in chrnames.txt). These are chromosome-end-associated read counts, not chromosome-specific telomere lengths.
4. Aggregation: combines all samples into stage2coverage.txt and ScaledGenomicOutput.txt.

Normalization

scaled count = raw count × 1,000,000,000 / totalReadCount

totalReadCount is every record in the input BAM, including unmapped and duplicate-flagged reads, before any qmotif filtering. We verified this on HG00096: qmotif reported 145,063,589 reads, identical to samtools view -c.

Performance

We benchmarked qmotif v1.0 against TelSeq, Telomerecat and TelomereHunter on five 1000 Genomes Phase 3 samples (four GBR, one FIN; low-coverage Illumina BAMs of 14–36 GB aligned to hg19). Each tool ran on one node with up to 8 threads (TelSeq is single-threaded) and a 20 GB memory limit. Times are wall-clock per sample.

Tool	Mean time (range)	Peak memory
qmotif v1.0	1.5 min (0.8–2.5)	6.7–7.4 GB
Telomerecat	3.3 min (2.1–4.8)	1.6–1.7 GB
TelSeq	11.2 min (6.2–17.1)	~5 MB
TelomereHunter	17.0 min (9.5–26.8)	~76 MB

qmotif was 2.3–11.6× faster than the other tools on average, at the cost of higher memory use. Peak usage was about 7 GB, well under the 20 GB Java heap limit we set (-Xmx20g). All runs completed without errors.

Disk: the downloaded BAM dominates storage (14–36 GB per sample here). qmotif's own outputs are small: the XML and extracted-reads BAM are under 0.3 MB each and the index about 1 MB. The 42 MB per-sample log comes from running at DEBUG log level; use --loglevel INFO to reduce it.

File Descriptions

• Snakefile: defines the full pipeline as a set of Snakemake rules, run per sample
• config.yaml: pipeline configuration: input/output directories, the FTP base URL for downloads, and the path to chrnames.txt
• files_list.txt: list of relative BAM/BAI paths for the samples to process
• chrnames.txt: list of chromosome names used for per-chromosome coverage tallying
• qmotif.ini: qmotif configuration file with the motif strings and regex, window size, and the telomeric INCLUDES regions used in the paper (hg19 coordinates)

Note: This pipeline previously ran as a set of standalone scripts (download_files.sh, runqmotif.py, stage2.py, realcoverage.sh, scaledgenomic.sh) executed manually in sequence. It has since been converted to a Snakemake pipeline: the Snakefile handles downloading, running qmotif, and aggregating results automatically, with per-sample dependency tracking and the ability to resume after failures.

Repository Structure

TeloTales/
├── docs/
│   ├── telotales.png
│   ├── consent_form_reference.pdf
│   ├── README_supplementary.md
│   └── examples/            (example reference outputs, e.g. HG00143)
├── metadata/
│   └── 1000genomes_sample_info.xlsx
├── notebooks/
│   ├── chapter1_final_across_chromosomes_kbps.rmd
│   └── Telomere_1k_genomes_project.rmd
├── results/
│   └── tables/               (analysis output CSVs)
├── workflow/
│   └── scripts/
│       ├── Snakefile
│       ├── config.yaml
│       ├── files_list.txt
│       ├── chrnames.txt
│       └── qmotif.ini         
├── LICENSE
└── README.md

Setup

1. Clone the repository

git clone https://github.com/paribytes/TeloTales.git
cd TeloTales/workflow/scripts

2. Create a conda environment for the pipeline

conda create -n telotales_env -c bioconda -c conda-forge snakemake
conda activate telotales_env

3. Check qmotif.ini The configuration used in the paper is included in workflow/scripts/. Edit it only if you want to change parameters or regions (e.g. a different reference build).

4. Confirm config.yaml config.yaml should point to your input/output directories and chrnames.txt:

input_dir: "data/bam"
output_dir: "results"
ftp_base: "https://ftp-trace.ncbi.nih.gov/1000genomes/ftp/phase3/data"
chromosomes_file: "chrnames.txt"

5. Provide files_list.txt

This should list relative BAM/BAI paths, one per line, matching the structure hosted on the 1000 Genomes FTP server, e.g.:

HG00096/alignment/HG00096.mapped.ILLUMINA.bwa.GBR.low_coverage.20120522.bam
HG00096/alignment/HG00096.mapped.ILLUMINA.bwa.GBR.low_coverage.20120522.bam.bai

Usage

Do a dry run first to check the pipeline is wired up correctly without downloading or computing anything:

snakemake --cores 1 -n

Then run for real:

snakemake --cores N

Replace N with the number of cores to use in parallel.

macOS note: if a shell: rule seems to hang indefinitely with no CPU activity, it may be a zsh-specific non-interactive shell quirk rather than an actual bug. Try forcing bash instead:

SHELL=/bin/bash snakemake --cores N

• Download Recommendations: For optimal performance, we recommend using an HPC or a Linux-based system, given the large file sizes involved (individual BAM files can be 15GB+). Ensure sufficient disk space and a stable internet connection before running this on a large sample list.

• For long-running jobs: on Linux or an HPC cluster, submit this as a background/queued job — see documentation on tmux and nohup for keeping long processes running across disconnected sessions. On macOS, prevent your machine from sleeping mid-run with:

caffeinate -i snakemake --cores N

If the pipeline stops partway through (e.g., due to a network interruption), just rerun the same snakemake command — completed steps are cached and won't be redone.

Outputs

• data/bam/{sample}.bam and .bam.bai - downloaded BAM/BAI files

• data/logs/{sample}_log.txt - qmotif log output

• data/qmotif/{sample}_output.txt - qmotif output, including scaled telomeric read counts

• results/stage2/{sample}_stage2_coverage.txt - chromosome-specific telomeric read counts

• results/stage2coverage.txt - combined telomeric read counts for each chromosome across all samples

• results/ScaledGenomicOutput.txt - scaled telomeric reads for all samples

• Example output files generated by the qmotif tool are included in the docs/examples directory. These files serve as reference outputs to understand the results produced by the qmotif analysis process.

Citation

• If you use TeloTales, please cite our paper

Shah, P., & Sethuraman, A. (2025). A Comprehensive Catalog of Telomere Content Variation across Human Populations. https://doi.org/10.1101/2025.11.03.686324 

Contact

• If you’d like to discuss this project or get in touch for other inquiries, please email me at priyanshishah213@gmail.com or connect with me on LinkedIn.