align with linked-read information

Arachne is the successor to lariat, the linked-read aware aligner originally developed by 10X Genomics. It incorporates linked-read information to better place alignments, which tends to show improvement when aligning over repetitive regions. The workflow first preprocesses FASTQ files into the format expected by Arachne (sorted by barcode, reads with invalid barcodes removed), then aligns the processed linked-reads with Arachne, while everything else gets aligned using minibwa.

Once sequences have been trimmed and passed through other QC filters, they will need to be aligned to a reference genome. This module within Harpy expects filtered reads as input, such as those derived using harpy qc . You can map reads onto a genome assembly with Harpy using the align arachne module:

harpy align arachne OPTIONS... REFERENCE INPUTS...

harpy align arachne genome.fasta Sequences/ 

Running Options

In addition to the common runtime options , the align bwa / align strobe modules are configured using these command-line arguments:

Output format

Regardless of the input linked-read format, the align workflows will standardize the output alignment records such that the barcode is contained in the BX:Z tag and barcode validation is in the VX:i tag.

Molecule distance

The --molecule-distance option is used during the alignment workflow in two places. First, it is used by Arachne directly when considering barcode information for alignment placement. Second, it uses the barcode information after aligning to deconvolute alignments with the same barcode that might not have originated from the same DNA molecule based on the distance threshold you specify. This happens during the linked-read stats step to internally split molecules based on this value, but it doesn't modify the barcodes in the output. Set this value to 0 to skip distance-based deconvolution during the this reporting step. Ignored if using --skip-reports.

Quality filtering

The --min-quality argument filters out alignments below a given MQ threshold. The default, 30, keeps alignments that are at least 99.9% likely correctly mapped. Set this value to 1 if you only want alignments removed with MQ = 0 (0% likely correct). You may also set it to 0 to keep all alignments for diagnostic purposes. The plot below shows the relationship between MQ score and the likelihood the alignment is correct and will serve to help you decide on a value you may want to use. It is common to remove alignments with MQ <30 (<99.9% chance correct) or MQ <40 (<99.99% chance correct).

Every alignment in a BAM file has an associated mapping quality score (MQ) that informs you of the likelihood that the alignment is accurate. This score can range from 0-40, where higher numbers mean the alignment is more likely correct. The math governing the MQ score actually calculates the percent chance the alignment is incorrect:

\%\ chance\ incorrect = 10^\frac{-MQ}{10} \times 100\\
\text{where }0\le MQ\le 40

You can simply subtract it from 100 to determine the percent chance the alignment is correct:

\%\ chance\ correct = 100 - \%\ chance\ incorrect\\
\text{or} \\
\%\ chance\ correct = (1 - 10^\frac{-MQ}{10}) \times 100

Marking PCR duplicates

Arachne marks duplicates internally, but reads with invalid barcodes get processed with samtools markdup to mark putative PCR duplicates. The read name is also parsed to determine if the sequencing platform was HiSeq/NovaSeq to distinguish between PCR and optical duplicates. Duplicate marking also uses the -S option to mark supplementary (chimeric) alignments as duplicates if the primary alignment was marked as a duplicate. Duplicates get marked but are not removed.

Arachne Workflow

graph LR
    A([index genome]):::clean --> B([align to genome]):::clean
    B-->C([sort alignments]):::clean
    C-->XX([standardize barcodes]):::clean
    XX-->D([mark duplicates]):::clean
    D-->E([assign molecules]):::clean
    E-->F([alignment metrics]):::clean
    D-->G([barcode stats]):::clean
    G-->F
    subgraph aln [Inputs]
        Z[FASTQ files]:::clean---genome[genome]:::clean
    end
    aln-->B & A
    subgraph markdp [mark duplicates via `samtools`]
        direction LR
        collate:::clean-->fixmate:::clean
        fixmate-->sort:::clean
        sort-->markdup:::clean
    end
    style markdp fill:#f0f0f0,stroke:#e8e8e8,stroke-width:2px,rx:10px,ry:10px
    style aln fill:#f0f0f0,stroke:#e8e8e8,stroke-width:2px,rx:10px,ry:10px
    classDef clean fill:#f5f6f9,stroke:#b7c9ef,stroke-width:2px

Align/{aligner}
├── Sample1.bam
├── Sample1.bam.bai
├── logs
│   ├── sample1.arachne.log
│   ├── sample1.markdup.log
│   │── sample1.sort.log
└── reports
    ├── barcodes.summary.ipynb
    ├── arachne.stats.ipynb
    ├── Sample1.ipynb
    └── data
        ├── lrstats
        │   └── Sample1.lrstats.gz
        └── coverage
            ├── Sample1.molcov.gz
            └── Sample1.cov.gz

minibwa map -y -x sr -R "@RG\tID:samplename\tSM:samplename"

strobealign [--use-index -r ...] -t THREADS -U -C --rg-id={sample} --rg=SM:{sample} {input.genome} {input.fastq}

minimap2 -t {threads} -a --MD -y -x map-{technology} -R "@RG\tID:samplename\tSM:samplename" {input.genome} {input.fastq}