deTASSELed

deTASSELed provides a Docker environment designed to run the TASSEL GBS pipeline reproducibly. It resolves dependency incompatibilities and version changes in newer Ubuntu environments by packaging all necessary dependencies, scripts, and binaries (including Beagle and PLINK) into a portable container.

This container is pre-configured with all necessary pipeline dependencies. To ensure maximum reproducibility, the software versions locked into the image are as follows:

• OS: Ubuntu 24.04.4 LTS (Noble Numbat)
• Java: OpenJDK 11
• Python: 3.12+ (standard library)
• TASSEL: 5.2.96 (March 31, 2025)
• BWA: 0.7.17-r1188
• Samtools: 1.19.2
• BCFtools: 1.19
• VCFtools: 0.1.16
• fastp: 0.23.4
• seqkit: 2.6.1
• PLINK2: v2.0.0-a.7.0LM (25 Apr 2026)
• Beagle: v5.5 (beagle.27Feb25.75f.jar)

This pipeline performs the following four major steps natively:

1. Discovery: Runs TASSEL GBSSeqToTagDBPlugin, TagExportToFastqPlugin, BWA alignment, SAMToGBSdbPlugin, and DiscoverySNPCallerPluginV2 to identify SNPs from FASTQ files.

2. Production: Uses the generated database to call genotypes across the taxa using ProductionSNPCallerPluginV2 and automatically reformats/indexes the resulting VCF.

3. Filtering: Runs a parallelized bash script using vcftools and bcftools to filter the raw VCF based on minor allele frequency, missingness, and heterozygosity at the taxa and SNP levels. Can also perform optional LD thinning.

4. REF/ALT Correction (Optional): Automatically invokes tassel_vcf_fixer.sh on the small filtered VCF to restore stripped TASSEL chromosome names and synchronize Major/Minor alleles into proper REF/ALT format before imputation.

5. Imputation: Uses Beagle to impute missing genotypes, handles re-headering the VCF file, and performs a final round of filtering on the imputed data.

The entire pipeline is operated by environment variables. By defining these variables at runtime, the container is directed to the correct data and settings.

• ENV_WORKDIR: The directory where the pipeline will store its output (e.g., /data/GBS-Runs).
• ENV_STUDY: A string representing the name of the study/project (e.g., virus_gwas_2025).
• ENV_KEYFILE: Path to the TASSEL keyfile.
• ENV_FASTQ: Directory containing the raw FASTQ sequence files.
• ENV_REF: Path to the reference genome FASTA (must be BWA-indexed).

By default, the pipeline orchestrator runs all steps sequentially. Users can toggle specific steps off if they only want to run a portion of the analysis:

• ENV_DISCOVERY="true" (Set to false to skip Discovery)
• ENV_PRODUCTION="true" (Set to false to skip Production)
• ENV_FILTER="true" (Set to false to skip Filtration)
• ENV_FIXREF="false" (Set to true to invoke REF/ALT Correction)
• ENV_IMPUTE="true" (Set to false to skip Imputation)

• ENV_DATABASE: Path to an existing discovery database (required if running Production without running Discovery first).
• ENV_MAP: Path to a Beagle map file. If left blank, Beagle will run with default settings without a genetic map.
• ENV_RAM: Memory allocated to Java processes (Default: 25g).
• ENV_NCORES: Number of CPU threads to utilize for parallel steps like BWA and Beagle (Default: 2).
• ENV_ENZYMES: Restriction enzymes used to generate the GBS library (Default: PstI-MspI).
• ENV_TAGLENGTH: The maximum length of tags used in the TASSEL database (Default: 64).
• ENV_MINQ: Minimum quality score within a tag (Default: 0).

Users can fine-tune the strictness of the filtering script with the following variables:

• ENV_TAXASUB: Text file indicating taxa to keep, one per line (Default: NA, performs no subsetting).
• ENV_TAXAMISS: Maximum missing data allowed per taxa (Default: 1).
• ENV_TAXAHET: Maximum heterozygosity allowed per taxa (Default: 1).
• ENV_SNPMISS: Maximum missing data allowed per SNP (Default: 1).
• ENV_SNPHET: Maximum heterozygosity allowed per SNP (Default: 1).
• ENV_MAF: Minimum minor allele frequency (Default: 0).
• ENV_MINDEP: Minimum read depth per SNP (Default: 0).
• ENV_MAXDEP: Maximum read depth per SNP (Default: 1e6).
• ENV_MINBP: Minimum physical distance between any two SNPs. This acts as a blunt physical-distance based pruning method (Default: 0).
• ENV_MAXLD: Maximum pairwise LD between SNPs. This performs smart statistical LD-based pruning over 50,000 SNP-wide windows (Default: NA).
• ENV_REMOVECHR: Chromosomes to remove in filtration (Default: NA).

System Requirements: Executing this build script requires a Unix-like operating system (Linux, macOS, or Windows via WSL) and requires docker to be installed and available in the system PATH.

To compile the reproducible environment, open a terminal at the root of this repository and execute the provided build script:

bash docker_build.sh

(Optional): If you intend to run this on an HPC cluster via Apptainer, you can pass the -B (or --build-tar) flag to automatically export the image to a .tar file after compiling:

./docker_build.sh -B

This script builds the detasseled-env image using Ubuntu 24.04 as a base, installs all offline .deb packages, and links all provided binaries to the system PATH.

The primary entry point inside the container is run_pipeline.sh. The method of execution depends on the deployment environment (Docker vs. Apptainer).

If running on a server with Docker privileges, you first need to load the exported tarball into the server's Docker registry. Transfer the detasseled-env.tar file to the server and run:

docker load -i detasseled-env.tar

Once the image is loaded, the pipeline can be executed directly using docker run. The -v flag must be used to "bind" the data directory into the container so it can read and write files.

A Docker template is provided in this repository as example_docker_run.sh. Fill in the paths at the top of the file and execute it:

bash example_docker_run.sh

HPC clusters like SCINet Atlas do not grant root privileges, so deployment must use Apptainer (formerly Singularity). To avoid issues with broken or missing system modules on HPC environments, we highly recommend installing your own version of Apptainer using Conda.

1. Setup Apptainer via Miniconda (Once): If you do not already have Miniconda installed (or if you are not using the module version), please install it into your default home directory (e.g. ~/miniconda3). Then, create the isolated Apptainer environment using the provided YAML file:

conda env create -f apptainer_env.yml

2. Transfer and Build the Image: The image cannot be built directly on the cluster using the Dockerfile. First, build the image on a local machine using the tarball flag:

bash docker_build.sh -B

Upload the resulting detasseled-env.tar file to the cluster. Then, activate your Conda environment and convert it into an Apptainer .sif file:

conda activate detasseled_apptainer
apptainer build detasseled-env.sif docker-archive://detasseled-env.tar

3. Submit via SLURM: Instead of passing -e flags in one massive command like Docker, HPC environments require users to export the ENV_ variables inside a .slurm batch script. The script should activate the Conda environment and run apptainer exec detasseled-env.sif run_pipeline.sh.

A SLURM template configured for SCINet Atlas (with Conda activation built-in) is provided in this repository as example_scinet_submission.slurm.

These scripts have been used to assist in TASSEL pipeline tasks in the past. We may add more in the future. Below is the definition of each and what they do.

This utility sanitizes corrupted FASTQ data and removes short reads that are known to fatally crash the TASSEL GBSSeqToTagDBPlugin (e.g., reads shorter than the ENV_TAGLENGTH). It is highly recommended to run this on raw data before starting the pipeline.

It performs a two-step process:

1. Sanitization: Uses seqkit sana to repair improperly formatted or truncated FASTQ records.
2. Filtering: Uses fastp to filter out reads below a user-defined length threshold.

Usage Arguments:

• -i Input directory of raw .fastq.gz files
• -o Output directory for the clean .fastq.gz files
• -t Threads (Default: 4)
• -l Minimum length (Default: 75)

Example Usage (Docker):

docker run --rm \
  -v /path/to/data:/path/to/data \
  detasseled-env fastq_sanitizer.sh -i /path/to/data/raw_fastq -o /path/to/data/clean_fastq -t 16 -l 75

Example Usage (Apptainer / SLURM):

apptainer exec detasseled-env.sif fastq_sanitizer.sh -i /path/to/data/raw_fastq -o /path/to/data/clean_fastq -t 16 -l 75

When downloading raw FASTQ data from sequencing facilities, files often arrive with long names. To ensure TASSEL and downstream tools parse these names correctly, users can utilize this interactive renaming utility.

To use it, run ./scripts/fastq_renamer_v1.03.sh on a local terminal (outside the container). The script will prompt the user for the working directory and the target file.

Failsafes:

1. Dynamic Header Parsing: It reads the gzipped sequence file (zcat | head) to extract the Flowcell ID and Lane Number embedded in the sequencing metadata, ensuring the new name is accurate.
2. Interactive Validation: It requires [y/n] user confirmation before modifying files.
3. Collision Prevention: Before performing the mv command, it checks the directory for existing files with the proposed name. If a duplicate is found, the script aborts to prevent overwriting sequencing data.

Example Usage (Local Terminal):

bash ./scripts/fastq_renamer_v1.03.sh

Example Usage (Docker): (Note: requires the -it flag for interactive prompts)

docker run -it --rm \
  -v /path/to/your/actual/data:/path/to/your/actual/data \
  detasseled-env fastq_renamer_v1.03.sh

Example Usage (Apptainer):

apptainer exec detasseled-env.sif fastq_renamer_v1.03.sh

TASSEL outputs VCF files with Major/Minor alleles instead of the standard Reference/Alternate (REF/ALT) format. Furthermore, TASSEL truncates chromosome names (e.g. converting Chr1A into 1).

This script resolves both issues. It expects ENV_VCFIN, ENV_REF, and ENV_VCFOUT to be set.

1. Dynamic Mapping (AWK): It scans the truncated chromosome names in the VCF, compares them to the headers in the FASTA file, replicates TASSEL's stripping logic (ignoring case, chr prefix, and suffixes), and builds a chromosome map.
2. Annotation: It runs bcftools annotate --rename-chrs to restore the chromosome names.
3. REF/ALT Flipping: It pipes the restored VCF through bcftools norm --check-ref ws -f <ref.fa>, checking every SNP against the FASTA and swapping Major/Minor into REF/ALT orientation.

Example Usage (Docker):

docker run --rm \
  -v /path/to/actual/data:/path/to/actual/data \
  -e ENV_VCFIN="/path/to/actual/data/TASSEL_output.vcf.gz" \
  -e ENV_REF="/path/to/actual/data/ref.fa" \
  -e ENV_VCFOUT="/path/to/actual/data/Final_Fixed.vcf.gz" \
  detasseled-env tassel_vcf_fixer.sh

Example Usage (Apptainer / SLURM):

Inside a .slurm submission script, export the required variables and run the script instead of the main pipeline orchestrator:

export ENV_VCFIN="/path/to/actual/data/TASSEL_output.vcf.gz"
export ENV_REF="/path/to/actual/data/ref.fa"
export ENV_VCFOUT="/path/to/actual/data/Final_Fixed.vcf.gz"

apptainer exec detasseled-env.sif tassel_vcf_fixer.sh

This utility extracts flanking genomic sequences around SNPs directly from a reference genome FASTA file for every variant present in a REF/ALT-corrected VCF file (e.g. following tassel_vcf_fixer.sh or post-imputation). It outputs a tab-delimited file formatted for PCR and KASP molecular marker design with the SNP locus centered in bracket notation (e.g. ATTAG[C/T]GTACG).

Key Features:

1. High Throughput: Streams extraction via bcftools query and bedtools getfasta, processing tens of thousands of loci in seconds.
2. Reference Validation: Verifies that the VCF REF allele strictly matches the base present in the reference FASTA at the specified position.
3. Chromosome Compatibility Check: Ensures chromosome naming between the VCF and FASTA match, flagging warnings if alleles require tassel_vcf_fixer.sh first.

Usage Arguments / Environment Variables:

• -i, --vcfin / ENV_VCFIN: Input VCF or BCF file (required)
• -r, --ref / ENV_REF: Reference genome FASTA file (required)
• -o, --outfile / ENV_OUTFILE: Output TSV file path (Default: <workdir>/<study>_snp_flanking_<length>bp.tsv)
• -l, --length / ENV_LENGTH: Flanking sequence length in bp on each side (Default: 100)
• -m, --allow-mismatch / ENV_ALLOW_MISMATCH: Retain variants even if REF base mismatches FASTA (Default: false)

Example Usage (Docker):

docker run --rm \
  -v /path/to/actual/data:/path/to/actual/data \
  -e ENV_VCFIN="/path/to/actual/data/Final_Fixed.vcf.gz" \
  -e ENV_REF="/path/to/actual/data/ref.fa" \
  -e ENV_LENGTH=100 \
  -e ENV_OUTFILE="/path/to/actual/data/kasp_flanking_sequences.tsv" \
  detasseled-env vcf_snp_sequence_puller.sh

Example Usage (Apptainer / SLURM):

export ENV_VCFIN="/path/to/actual/data/Final_Fixed.vcf.gz"
export ENV_REF="/path/to/actual/data/ref.fa"
export ENV_LENGTH=100
export ENV_OUTFILE="/path/to/actual/data/kasp_flanking_sequences.tsv"

apptainer exec detasseled-env.sif vcf_snp_sequence_puller.sh

Converts LGC KlusterCaller KASP genotyping cluster calls (X:X, Y:Y, X:Y, No Call) into a standardized, sorted, bgzipped, and indexed VCF (.vcf.gz) using a high-performance Python engine.

Key Features:

1. Explicit True Reference Assignment (ref_allele): Eliminates arbitrary allele assignment by letting users denote whether probe X or Y (or a specific nucleotide base) corresponds to the true reference allele.
2. Strand Orientation & Complementation (strand): Supports tracking leading (+) vs lagging (-) assay design strands, automatically reverse-complementing probe alleles to the reference genome's forward strand when needed.
3. Unmapped Marker Handling: Automatically assigns sequential pseudo-positions for markers with missing coordinates (NA), flagging them with UNMAPPED_PSEUDO_POS=1 so downstream tools do not crash.
4. CSI Indexing: Automatically generates .csi indexes for compatibility with large crop genomes (like wheat chromosomes >512 Mb).

Keyfile Format (Tab-Delimited):

• marker: Marker identifier (matches KlusterCaller column header)
• chr: Chromosome identifier (e.g. Chr1B, Chr4A)
• position: 1-based genomic coordinate (or NA)
• X: Nucleotide for Probe X / FAM (e.g. C)
• Y: Nucleotide for Probe Y / VIC/HEX (e.g. T)
• ref_allele: (Optional) X or Y (Default: X)
• strand: (Optional) + (leading/forward) or - (lagging/reverse) (Default: +)
• No Call: (Optional) Character for missing calls (Default: N)

Example Usage (Docker):

docker run --rm \
  -v /path/to/actual/data:/path/to/actual/data \
  -e ENV_KLUSTER="/path/to/actual/data/klustercaller_calls.txt" \
  -e ENV_KEYFILE="/path/to/actual/data/kasp_keyfile.tsv" \
  -e ENV_VCFOUT="/path/to/actual/data/kasp_calls.vcf.gz" \
  detasseled-env klustercaller_to_vcf.sh

Example Usage (Apptainer / SLURM):

export ENV_KLUSTER="/path/to/actual/data/klustercaller_calls.txt"
export ENV_KEYFILE="/path/to/actual/data/kasp_keyfile.tsv"
export ENV_VCFOUT="/path/to/actual/data/kasp_calls.vcf.gz"

apptainer exec detasseled-env.sif klustercaller_to_vcf.sh