Preprocess Raw Sequences

Sequences that come off of a sequencer are typically not ready for analytical use. If the data hasn't already been separated out by sample, you need to do that first. If you're using linked-read data, you then need to identify the linked-read barcodes and pull them out of the sequence. Below describes the workflows for the original haplotagging protocol published by Meier et al. (2021), along with the one we're developed at Cornell University (Iqbal et al., in prep).

gih

harpy preprocess gih OPTIONS... INPUTS

harpy preprocess gih --threads 20 data/*R{1,2}.fastq.gz

Running Options

In addition to the common runtime options , the preprocess meier2021 module is configured using these command-line arguments:

The ME Sequence

The ME sequence is very important for preprocessing these data correctly, as identifying the position it starts will allow the algorithm to determine how much to pad the sequence (to offset the variable stagger). The GIH at Cornell University uses the AGATGTGTATAAGAGACAG, which is why it is set to the default.

The default --mismatch value tolerate 2 mismatches in the ME sequence before determining the ME sequence is absent. The N nucleotide counts as a partial mismatch (0.3), whereas the standard ATCG bases count as a single mismatch. This means 1 ATCG mismatch and 2 Ns would still pass with --mismatch 2:

\text{-m} = AGATGTGTATAAGAGACAG

\text{sequence} = \hat GGATGTGTATAAGA\hat NACA\hat N

\text{mismatch} = 1_{ATCG} + 0.3_{N} + 0.3_{N} = 1.6

Workflow

graph LR
    A[multiplexed FASTQ]:::clean-->B([pad the stagger]):::clean
    B-->C([pheniqs demultiplex]):::clean
    C-->D([convert to standard]):::clean
    D-->E([quality metrics]):::clean
    classDef clean fill:#f5f6f9,stroke:#b7c9ef,stroke-width:2px

Preprocess/
├── Sample1.F.fq.gz
├── Sample1.R.fq.gz
├── Sample2.F.fq.gz
├── Sample2.R.fq.gz
└── reports
    ├── performance.QC.ipynb
    └── preprocess.QC.html

meier2021

harpy preprocess meier2021 OPTIONS... R1_FQ R2_FQ I1_FQ I2_FQ

harpy preprocess meier2021 --threads 20 --schema demux.schema Plate_1_S001_R*.fastq.gz Plate_1_S001_I*.fastq.gz

Running Options

In addition to the common runtime options , the preprocess meier2021 module is configured using these command-line arguments:

Keeping Unknown Samples

It's not uncommon that some sequences cannot be demultiplexed due to sequencing errors at the ID location. Use --keep-unknown-samples/-u to have Harpy still separate those reads from the original multiplex. Those reads will be labelled _unknown_sample.R*.fq.gz

Keeping Unknown Barcodes

It's likewise not uncommon that sequencing errors make it so that the sequences don't match the list of known barcode segments. Use --keep-unknown-barcodes/-b to have Harpy separate those reads out from the original multiplex as _unknown_barcodes.R*.fq.gz.

Keep QX and RX Tags

Using --qx-rx, you can opt-in to retain the QX:Z (barcode PHRED scores) and RX:Z (nucleotide barcode) tags in the sequence headers. These tags aren't used by any subsequent analyses, but may be useful for your own diagnostics.

Demultiplexing Schema

Generation I haplotags typically use a unique Cxx barcode per sample-- that's the barcode segment that will be used to identify sequences by sample. However, any of the 4 segments (A,B,C,D) are valid, so long as the schema only features a single segment. You will need to provide a simple text file to --schema (-s) with two columns, the first being the sample name, the second being the identifying segment barcode (e.g., C19). This file is to be tab or space delimited and must have no column names.

Sample01    C01
Sample02    C02
Sample03    C03
Sample04    C04

This will result in splitting the multiplexed reads into individual file pairs Sample01.F.fq.gz, Sample01.R.fq.gz, Sample02.F.fq.gz, etc. A sample can have multiple barcodes, but a barcode cannot have multiple samples:

Sample01    D01
Sample02    D02
Sample03    D03
Sample03    D21

Sample01    C01
Sample02    C02
Sample03    C02

Sample01    C01
Sample02    D02
Sample03    C03

Workflow

graph LR
    subgraph Inputs
        direction TB
        A[multiplexed FASTQ]:::clean---BX
        BX[index reads FASTQ]:::clean---SCH
        SCH[Sample Schema]:::clean
    end
    Inputs-->B([demultiplex samples]):::clean
    B-->D([quality metrics]):::clean
    style Inputs fill:#f0f0f0,stroke:#e8e8e8,stroke-width:2px,rx:10px,ry:10px
    classDef clean fill:#f5f6f9,stroke:#b7c9ef,stroke-width:2px

Preprocess/
├── Sample1.F.fq.gz
├── Sample1.R.fq.gz
├── Sample2.F.fq.gz
├── Sample2.R.fq.gz
└── reports
    └── preprocess.QC.html

The power of dmox!

Harpy v2 introduced a new demultiplexer under the hood called dmox, which is singificantly faster, lighter on memory, and has better maintenance than the previous solution. Iago Bonnici of Montpellier Bioinformatics Biodiversity (MBB) saw the need for better demultiplexing performance and took it upon themselves to donate their time to write a brand-new purpose-built demultiplexer for the Meier/Chan haplotagging bead design. Beyond just being way more performant, this new demultiplexer has more features, has more output options, and is flexible for haplotagging bead designs where the sample ID is not the C-segment. If you're happy with the performance of the new demultiplexing workflow, please let Iago/MBB know!