Zubair Khalid

Virologist/Molecular Biologist | Veterinarian | Bioinformatician

Conventional & Molecular Virology • Vaccine Development • Computational Biology

Dr. Zubair Khalid is a veterinarian and virologist specializing in conventional and molecular virology, vaccine development, and computational biology. Dedicated to advancing animal health through innovative research and multi-omics approaches.

Dr. Zubair Khalid - Veterinarian, Virologist, and Vaccine Development Researcher specializing in Computational Biology, Multi-omics, Animal Health, and Infectious Disease Research

Category: Guides

Sequencing Pipeline

A sequencing pipeline is a defined computational workflow that converts raw sequencing reads into biological insights. This guide is for researchers, bioinformatics trainees, and lab scientists who need a practical framework to design, run, and troubleshoot their own pipelines. Sources like NCBI Bookshelf provide foundational references for the concepts discussed here.

The core goal of any pipeline is to accurately process high throughput sequencing data while controlling for errors and biases. The Galaxy Training Network offers hands on tutorials that mirror the steps outlined in this guide.

At a Glance

Stage Description Common Tools Critical Output
Raw data quality control Assess base quality, adapter content, and read length distributions FastQC, MultiQC Quality report and trimming decisions
Read alignment Map reads to a reference genome or transcriptome BWA, STAR, Bowtie2 SAM/BAM alignment files
Post alignment processing Sort, mark duplicates, index, and recalibrate base scores SAMtools, Picard, GATK Cleaned, sorted BAM files
Variant calling or quantification Identify variants or count features (genes, transcripts) GATK HaplotypeCaller, featureCounts, HTSeq VCF files or count matrices
Annotation and interpretation Assign biological meaning to variants or expression changes ANNOVAR, SnpEff, Ensembl VEP Annotated tables, reports

Core Concepts and Decision Points

A sequencing pipeline must match the experimental design. The first decision is the type of sequencing: whole genome shotgun, whole exome, targeted amplicon, RNA seq, ChIP seq, or bisulfite sequencing, among others. Each requires different downstream logic. The EMBL EBI Training resources break down these categories with clear learning paths.

Key decision criteria include:

  • Reference versus de novo assembly. If a high quality reference exists for the species, alignment based pipelines are typical. For non model organisms or metagenomic samples, de novo assembly may be necessary.
  • Read length and pairing. Illumina short reads (paired end) dominate most pipelines today. Long read platforms (PacBio, Oxford Nanopore) require specialized base calling and alignment tools, often with higher error rates that must be corrected.
  • Library preparation strand specificity. For RNA seq pipelines, strand specific protocols change how quantification is performed. Ignoring strandedness leads to misassignment of antisense reads.
  • Duplex versus single unique molecular identifier (UMI). UMIs allow accurate duplicate removal and error correction in amplicon or low frequency variant detection.

The pipeline is not a single universal script. You must adapt the order of steps and parameter settings based on the biological question and data source.

Practical Workflow

The following sequence outlines a typical pipeline for Illumina short read variant calling. Adapt steps for other assay types as needed.

1. Obtain raw sequencing data. Raw reads are stored in FASTQ format. Public data can be downloaded from the NCBI Sequence Read Archive. Ensure you have the correct SRA accession. Use fastq dump or fasterq dump to convert SRA files to FASTQ.

2. Quality control with FastQC. Run FastQC on each FASTQ file to examine per base quality scores, GC content, sequence duplication levels, adapter contamination, and overrepresented sequences. MultiQC can aggregate multiple FastQC reports into a single summary. If adapters are detected, plan to trim them.

3. Trimming and filtering. Use tools like Trimmomatic, cutadapt, or fastp to remove adapter sequences, low quality bases (often Q score < 20), and short reads (e.g., < 30 bases). Record trimming parameters for reproducibility.

4. Alignment to reference genome. Choose an aligner suited to your data. For short reads mapping to a linear genome, BWA mem or Bowtie2 are standard. For RNA seq with splice junctions, STAR is preferred. Aligners output SAM format, which should be converted to sorted BAM using SAMtools.

5. Post alignment processing. Mark duplicate reads using Picard MarkDuplicates (for PCR based libraries) or GATK MarkDuplicates. For WGS and WES pipelines, base quality score recalibration (BQSR) in GATK improves variant calling accuracy. Index the final BAM file with SAMtools index.

6. Variant calling. Use a caller that matches the study design. GATK HaplotypeCaller is common for germline single nucleotide variants and small indels. For RNA seq variant calling, consider GATK with appropriate splicing aware parameters. Output is a VCF (Variant Call Format) file.

7. Annotation and filtering. Annotate variants with functional impact using ANNOVAR, SnpEff, or Ensembl VEP. Add population frequency filters (e.g., gnomAD), predicted pathogenicity scores (CADD, SIFT, PolyPhen), and clinical significance databases. Filter to focus on rare, high impact variants.

For RNA seq quantification pipelines, replace steps 6 and 7 with feature counting (featureCounts or HTSeq) or transcript abundance estimation (Salmon, Kallisto). The Bioconductor project provides extensive R packages for differential expression analysis after counting.

Quality Checks and Verification

Quality control is not a single step but a continuous process throughout the pipeline.

  • Base quality distribution. After trimming, re run FastQC on cleaned reads. The median per base quality should be above Q30 for most positions.
  • Alignment rates. For human WGS, at least 95% of reads should map to the reference. Low alignment may indicate contamination, adapter issues, or wrong reference.
  • Insert size metrics (paired end). Evaluate the insert size distribution from the BAM file using Picard CollectInsertSizeMetrics. Consistent library size supports proper alignment.
  • Coverage depth. Use SAMtools depth or GATK to compute mean coverage and the fraction of target bases covered at 10x, 20x, 30x thresholds. For clinical variant detection, minimum coverage thresholds often apply.
  • Variant transition transversion ratio. For whole genome germline data, a Ti/Tv ratio around 2.0 to 2.1 is expected. Deviations suggest systematic errors.

Verify that computational environment does not introduce batch effects. Use principal component analysis on variant or count data to check for separation by sequencing run or library batch.

Common Mistakes to Avoid

  • Skipping adapter trimming. Adapter sequences cause alignment artifacts and false positive variant calls. Always run a trimmer even if FastQC shows no adapters. Some protocols use indexed adapters that may not appear as standard overrepresented sequences.
  • Using the wrong reference genome. Different builds (hg19 vs GRCh38) have different coordinates and alt loci variants. Mismatch leads to unreproducible results. Check assembly accession against the source data descriptor.
  • Ignoring read strandedness in RNA seq. Strand information is required for accurate quantification of antisense transcripts. Use tools that infer or accept an explicit library type (e.g., RF, FR).
  • **Not removing PCR duplicates. Duplicates inflate coverage of specific molecules and bias allele frequency estimates. Deduplication is mandatory for variant calling in amplicon or low input libraries.
  • Applying hard filters without context. For example, removing all variants with low quality scores may eliminate true positives. Use variant quality score recalibration (VQSR) in GATK with a training set of known variants rather than arbitrary thresholds.
  • Merging technical replicates without checking concordance. Replicates should be compared for correlation before merging. High discordance indicates a batch effect or library problem.

A study of parental postzygotic mutations across rare disease trios Landscape of parental postzygotic mutations highlights how sensitive variant detection is to read depth and duplicate removal. Pipelines that skip these checks miss or overcall de novo mutations.

Limits and Uncertainty

Every sequencing pipeline introduces uncertainty that must be understood before biological interpretation.

Coverage and variant detection sensitivity. Low coverage regions lead to false negative variant calls. Even at high mean coverage, some genomic regions remain inaccessible due to GC bias, repetitive elements, or structural variation. Structural variant callers still have high false positive rates across different algorithms.

Mapping ambiguity. Multimapping reads (reads that align equally well to multiple genomic locations) are often discarded or assigned arbitrarily. This can distort quantification for repetitive genes or transposable elements.

Base quality score calibration. Sanger type quality scores from Illumina machines may be systematically over or under confident. Recalibration (BQSR) can reduce but not eliminate errors, especially near indels or homopolymer runs.

Validation requirements. Computational variant calls remain predictions. Clinical and functional follow up, such as Sanger sequencing or orthogonal genotyping, is required for high stakes applications. The field of artificial intelligence in clinical genetics Artificial Intelligence in Clinical Genetics shows that even advanced machine learning models depend on high quality variant calling from traditional pipelines.

Interpretation of non coding variants. Most pipelines focus on coding regions. Non coding variants are abundant but their functional prediction is far less standardized. Methods from evolutionary fusion analyses Evolutionary innovation through fusion may offer new ways to assess non coding elements, but they remain research grade.

Always report the exact software versions, parameters, and reference database versions used. Pipelines should be containerized (e.g., Docker, Singularity) or use workflow managers (Nextflow, Snakemake) to ensure reproducibility across environments.

Frequently Asked Questions

What is the most common file format in a sequencing pipeline? The FASTQ format is universal for raw reads. After alignment, BAM (binary SAM) becomes the standard. Variant calls are stored in VCF. All these formats are documented on the NCBI Bookshelf.

How long does a typical pipeline take to run? Time depends on data volume, read length, and compute resources. Human whole genome sequencing with 30x coverage can take 24 to 48 hours on a server with 16 cores and 64 GB RAM. Cloud or cluster backends speed up parallel steps.

Do I need to use a container or virtual environment? Yes. Software version differences are a leading cause of irreproducibility. Use Conda environments, Docker images, or a pipeline manager like Nextflow with built in container support. The Galaxy Training Network runs entire pipelines in controlled environments.

Can I reuse the same pipeline for different species? The conceptual steps are similar, but reference genome quality, aligner parameters, and variant calling databases differ. For non model species, you may need to build a de novo assembly first or use a closely related reference. For bacterial genomes, speed optimized mappers like Minimap2 are more appropriate.

References and Further Reading

Related Articles