diff --git a/.github/filters.yml b/.github/filters.yml index c999fd5f8..f900df24f 100644 --- a/.github/filters.yml +++ b/.github/filters.yml @@ -20,8 +20,8 @@ validate: &validate - 'test/bam/**' - 'harpy/utils/check_bam.py' - 'harpy/utils/check_fastq.py' - - 'harpy/report/notebooks/validate_fastq.ipynb' - - 'harpy/report/notebooks/validate_bam.ipynb' + - 'harpy/notebooks/validate_fastq.ipynb' + - 'harpy/notebooks/validate_bam.ipynb' deconvolve: &deconvolve - *common - *environments @@ -44,15 +44,20 @@ qc: &qc - 'harpy/utils/bx_stats_fq.py' - 'harpy/validation/fastq.py' - 'test/fastq/**' + - 'harpy/notebooks/fastp_qc.ipynb' + - 'harpy/notebooks/qc_bx_stats.ipynb' align: &align - *common - *environments + - 'harpy/utils/standardize/standardize.go' - 'harpy/commands/align.py' - 'harpy/snakefiles/align_bwa.smk' - 'harpy/snakefiles/align_strobe.smk' - - 'harpy/report/notebooks/align_stats.ipynb' - - 'harpy/report/notebooks/align_lrstats.ipynb' - - 'harpy/report/notebooks/samtools_stats.ipynb' + - 'harpy/snakefiles/align_minimap.smk' + - 'harpy/snakefiles/align.smk' + - 'harpy/notebooks/align_stats.ipynb' + - 'harpy/notebooks/align_lrstats.ipynb' + - 'harpy/notebooks/samtools_stats.ipynb' - 'harpy/utils/bx_stats_sam.py' - 'harpy/utils/molecule_coverage.py' - 'harpy/validation/fasta.py' @@ -64,7 +69,7 @@ snp: &snp - 'harpy/commands/snp.py' - 'harpy/snakefiles/snp_mpileup.smk' - 'harpy/snakefiles/snp_freebayes.smk' - - 'harpy/report/notebooks/bcftools_stats.ipynb' + - 'harpy/notebooks/bcftools_stats.ipynb' - 'harpy/validation/fasta.py' - 'harpy/validation/xam.py' - 'test/bam/**' @@ -73,8 +78,8 @@ impute: &impute - *environments - 'harpy/commands/impute.py' - 'harpy/snakefiles/impute.smk' - - 'harpy/report/notebooks/impute.ipynb' - - 'harpy/report/notebooks/stitch_collate.ipynb' + - 'harpy/notebooks/impute.ipynb' + - 'harpy/notebooks/stitch_collate.ipynb' - 'harpy/validation/fasta.py' - 'harpy/validation/xam.py' - 'test/bam/**' @@ -84,7 +89,7 @@ leviathan: &leviathan - *environments - 'harpy/commands/sv.py' - 'harpy/snakefiles/sv_leviathan.smk' - - 'harpy/report/notebooks/sv.ipynb' + - 'harpy/notebooks/sv.ipynb' - 'harpy/validation/fasta.py' - 'harpy/validation/xam.py' - 'test/bam/**' @@ -93,7 +98,7 @@ naibr: &naibr - *environments - 'harpy/commands/sv.py' - 'harpy/snakefiles/sv_naibr.smk' - - 'harpy/report/notebooks/sv.ipynb' + - 'harpy/notebooks/sv.ipynb' - 'harpy/utils/infer_sv.py' - 'harpy/validation/fasta.py' - 'harpy/validation/xam.py' @@ -109,7 +114,7 @@ phase: &phase - 'test/bam/**' - 'test/vcf/test.bcf' - 'test/vcf/test.phased.bcf' - - 'harpy/report/notebooks/hapcut.ipynb' + - 'harpy/notebooks/hapcut.ipynb' - 'harpy/utils/parse_phaseblocks.py' assembly: &assembly - *common diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index d76a436f4..3cebb33c4 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -190,6 +190,10 @@ jobs: run: | harpy align bwa --quiet 2 -x "-A 2" test/genome/genome.fasta.gz test/fastq && \ ls -lh Align/bwa + - name: test minimap + run: | + harpy align minimap --quiet 2 test/genome/genome.fasta.gz test/fastq && \ + ls -lh Align/minimap snp: needs: changes diff --git a/CHANGELOG.md b/CHANGELOG.md index 80e43abe2..a8feab23b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,26 @@ # New ## QC - replace multiqc report with native harpy report +## Align +- BWAMEM2 has been replaced with minibwa. Long live BWA! + - it's much faster, and takes much less time to index a reference +- coverage depth added to aggregate report for processed alignments +- minimap2 added back in (with extra perks) for long-read compatability + - called with `harpy align minimap` + +## Report +- harpy reports can be converted to less-nice but functional standalone HTML files + - this feature is accessed using `harpy report static` + - to accomodate this, `harpy report` (live report website) is now `harpy report live` + +# Changes +## Align +- `-d` (molecule distance threshold) has its default restored to 50kb since this value is used exclusively for reporting and does not alter the data + +## misc +- removed FASTA format validation because it can be dreadfully slow with existing tools + + ## Report - harpy reports can be converted to less-nice but functional standalone HTML files @@ -14,4 +34,9 @@ ## misc - constrain CASAVA regex in FASTQ file validation so it doesn't trigger false positives when new CASAVA appears in unexpected places - [internal] notebooks no longer a submodule/subdirectory of `harpy.report` -- utility `check_fastq.py` no longer employs globals, instead uses a sensible class system \ No newline at end of file +- utility `check_fastq.py` no longer employs globals, instead uses a sensible class system +- add multithreading to pre-workflow VCF and XAM file validation and parsing +- [internal-ish] the bwa, strobealign, and minimap2 workflows are nearly identical except for the reference preprocessing and alignment, so to minimize redundancy and duplication, those workflows have a single `align.smk` that imports a second snakefile `align_{aligner}.smk` that handles just the preprocessing and direct alignment for those aligners, then hands off to `align.smk` for all the downstream things (dedup, sorting, reports, etc) + +# Documentation +- the pages for bwa, strobealign, and minimap have been consolidated into a single page bc they are nearly identical diff --git a/docs/Commands/align/align.md b/docs/Commands/align/align.md index af8073275..2974aca09 100644 --- a/docs/Commands/align/align.md +++ b/docs/Commands/align/align.md @@ -9,18 +9,18 @@ will need to align them to a reference genome before you can call variants. Harpy offers several aligners for this purpose: {.compact .clean .whitespace-nowrap} -| aligner | speed | repository | publication | -|:-------------------------|:-------------:|:------------------------------------------------:|:---------------------------------------------------:| -| [BWA](bwa.md) | fast ⚡ | [github](https://github.com/lh3/bwa) | [paper](http://arxiv.org/abs/1303.3997) | -| [strobealign](strobe.md) | super fast ⚡⚡ | [github](https://github.com/ksahlin/strobealign) | [paper](https://doi.org/10.1186/s13059-022-02831-7) | +| command | aligner | best for | repository | publication | +| :------ | :---------- | :---------- | -----------------------------------------------: | :----------------------------------------------------: | +| bwa | minibwa | general use | [github](https://github.com/lh3/minibwa) | [preprint](https://github.com/lh3/minibwa) | +| minimap | minimap2 | long reads | [github](https://github.com/lh3/minimap2) | [paper](https://doi.org/10.1093/bioinformatics/bty191) | +| strobe | strobealign | speed | [github](https://github.com/ksahlin/strobealign) | [paper](https://doi.org/10.1186/s13059-022-02831-7) | -Neither of these are linked-read aware aligners, but Harpy transfers the barcode information from the sequence headers into the alignments and will -assign molecule identifiers (`MI:i` SAM tags) based on these barcodes and the [molecule distance threshold](../../Getting_Started/linked_read_data.md#barcode-thresholds). +Neither of these are linked-read aware aligners, but Harpy transfers the barcode information from the sequence headers into the alignments. ## Non linked-read WGS data Starting with Harpy `v2.x`, you can skip the workflow routines that do things specific to linked reads, meaning you can comfortably use -[!badge corners="pill" text="harpy align bwa"](bwa.md) and [!badge corners="pill" text="harpy align strobe"](strobe.md) to align your WGS sequence data. +[!badge corners="pill" text="harpy align bwa"](standard.md) and [!badge corners="pill" text="harpy align strobe"](standard.md) to align your WGS sequence data. - version `2.0-2.7` : `--ignore-bx` - version `>2.7` : `--lr-type none` - version `>=3.0`: autodetected or forced with `--unlinked` diff --git a/docs/Commands/align/arachne.md b/docs/Commands/align/arachne.md new file mode 100644 index 000000000..b11623e57 --- /dev/null +++ b/docs/Commands/align/arachne.md @@ -0,0 +1,327 @@ +--- +label: arachne +description: Align sequences using linked-read information with arachne +category: [linked-read] +tags: [linked-read] +icon: dot +order: 5 +hidden: true +--- + +# :icon-quote: 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. + + +=== :icon-checklist: You will need +- at least 4 cores/threads available +- a genome assembly in FASTA format: [!badge variant="success" text=".fasta"] [!badge variant="success" text=".fa"] [!badge variant="success" text=".fasta.gz"] [!badge variant="success" text=".fa.gz"] [!badge variant="secondary" text="case insensitive"] +- paired-end fastq sequence files [!badge variant="secondary" icon=":heart:" text="gzipped recommended"] + - **sample name**: [!badge variant="success" text="a-z"] [!badge variant="success" text="0-9"] [!badge variant="success" text="."] [!badge variant="success" text="_"] [!badge variant="success" text="-"] [!badge variant="secondary" text="case insensitive"] + - **forward**: [!badge variant="success" text="_F"] [!badge variant="success" text=".F"] [!badge variant="success" text=".1"] [!badge variant="success" text="_1"] [!badge variant="success" text="_R1_001"] [!badge variant="success" text=".R1_001"] [!badge variant="success" text="_R1"] [!badge variant="success" text=".R1"] + - **reverse**: [!badge variant="success" text="_R"] [!badge variant="success" text=".R"] [!badge variant="success" text=".2"] [!badge variant="success" text="_2"] [!badge variant="success" text="_R2_001"] [!badge variant="success" text=".R2_001"] [!badge variant="success" text="_R2"] [!badge variant="success" text=".R2"] + - **fastq extension**: [!badge variant="success" text=".fq"] [!badge variant="success" text=".fastq"] [!badge variant="secondary" text="case insensitive"] +=== + +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 [!badge corners="pill" text="harpy qc"](../qc.md). You can map reads onto a genome assembly with Harpy using the [!badge corners="pill" text="align arachne"] module: + +```bash usage +harpy align arachne OPTIONS... REFERENCE INPUTS... +``` +```bash example +harpy align arachne genome.fasta Sequences/ +``` + +## :icon-terminal: Running Options +In addition to the [!badge variant="info" corners="pill" text="common runtime options"](/Getting_Started/common_options.md), the [!badge corners="pill" text="align bwa"]/[!badge corners="pill" text="align strobe"] modules are configured using these command-line arguments: + +{.compact .clean} +| argument {.whitespace-nowrap} | default {.whitespace-nowrap} | description | +| :------------------------------- | :--------------------------: | :---------------------------------------------------------------------------------------------------------------------------------------------- | +| `REFERENCE` | | [!badge variant="info" text="required"] Reference assembly for read mapping | +| `INPUTS` | | [!badge variant="info" text="required"] Files or directories containing [input FASTQ files](/Getting_Started/common_options.md#input-arguments) | +| `-depth-window` `-w` | `50000` | Interval size (in bp) for depth stats | +| `--extra-params` `-x` | | Additional BWA arguments, in quotes | +| `--keep-unmapped` `-u` | false | Output unmapped sequences too | +| `--molecule-distance` `-d` | `0` | Base-pair distance threshold to separate molecules given as base pairs, disabled with `0` | +| `--min-quality` `-q` | `30` | Minimum `MQ` (SAM mapping quality) to pass filtering | + + +### 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](/Getting_Started/linked_read_data.md#barcode-thresholds) +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). + +==- What is the $MQ$ score? +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 +$$ + +![A visual explanation of MQ Score](/static/MQscore.png) +=== + +## 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**. + +---- + +## :icon-git-pull-request: Arachne Workflow + ++++ :icon-git-pull-request: workflow +```mermaid +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 +``` + ++++ :icon-file-directory: output +The default output directory is `Align/arachne` with the folder structure below. +`Sample1` is a generic sample name for demonstration purposes. The resulting folder also includes a `workflow` directory +(not shown) with workflow-relevant runtime files and information. +``` +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 +``` +{.compact} +| item {.whitespace-nowrap} | description | +| :---------------------------------- | :------------------------------------------------------------------------------------- | +| `*.bam` | sequence alignments for each sample | +| `*.bai` | sequence alignment indexes for each sample | +| `logs/*arachne.log` | output of the aligner during run | +| `logs/*markdup.log` | stats provided by `samtools markdup` | +| `logs/*sort.log` | output of `samtools sort` | +| `reports/` | various counts/statistics/reports relating to sequence alignment | +| `reports/barcodes.summary.ipynb` | report summarizing barcode-specific metrics across all samples | +| `reports/arachne.stats.ipynb` | report summarizing `samtools stats` of raw and processed alignments across all samples | +| `reports/Sample1.ipynb` | report summarizing BX tag metrics and alignment coverage | +| `reports/data/coverage/*.cov.gz` | output from mosdepth, used for reports | +| `reports/data/coverage/*.molcov.gz` | molecular coverage stats, used for reports | +| `reports/data/lrstats` | tabular data containing the information used to generate the BX stats in reports | ++++ + ++++ :icon-git-merge: details +- incorporates barcode information + +The [arachne](https://github.com/pdimens/arachne) workflow maps all reads against the reference genome. +Duplicates for valid-barcoded reads are marked internally by arachne, while invalid-barcoded reads use `samtools markdup`. +The `-m` threshold is used for alignment molecule assignment during arachne aligning and when calculating statistics. + ++++ :icon-code-square: minbwa parameters +By default, Harpy runs `arachne` with these parameters (excluding inputs and outputs): +```bash +minibwa map -y -x sr -R "@RG\tID:samplename\tSM:samplename" +``` + +Below is a list of all `minibwa map` command line arguments, excluding those Harpy already uses or those made redundant by Harpy's implementation of BWA. + +{.compact .clean} +| argument {.whitespace-nowrap} | category {.whitespace-nowrap} | description | +| :---------------------------- | :-----------------------------: | :--------------------------------------------------------------------- | +| `-l` | common | treat reads 80% of the best hit, output them to XA [5] | +| `-Y` | file IO | use soft clipping for supplementary alignments | +| `-H` | file IO | if STR starts with @, insert to header; or insert lines in file STR [] | +| `-5` | file IO | take the alignment with the smallest query position as primary | +| `-K` | file IO | process NUM1-NUM2 bp of query sequences in a batch [100m,1g] | +| `--mmap[=lite]` | file IO | load the index via memory mapped files (slower mapping) [] | + ++++ + +==- strobealign ++++ :icon-git-merge: details +- ignores barcode information (but retains in output) +- ultra-fast +- [as-good-or-better accuracy](https://github.com/ksahlin/strobealign/blob/main/evaluation.md) to BWA MEM for sequences greater than 100bp + - accuracy may be lower for sequences shorter than 100bp + +The [strobealign](https://github.com/lksahlinh3/strobealign) workflow is nearly identical to the BWA workflow, +the only real difference being how the input genome is indexed and that alignment is performed with +`strobealign` instead of BWA. Duplicates are marked using `samtools markdup`. +The `BX:Z` tags in the read headers are still added to the alignment headers, even though barcodes +are not used to inform mapping. The `-m` threshold is used for alignment molecule assignment. + + ++++ :icon-code-square: strobealign parameters +By default, Harpy runs `strobealign` with these parameters (excluding inputs and outputs): +```bash +strobealign [--use-index -r ...] -t THREADS -U -C --rg-id={sample} --rg=SM:{sample} {input.genome} {input.fastq} +``` + +Below is a list of all `strobealign` command line arguments, excluding those Harpy already uses or those made redundant by Harpy's implementation of it. + +{.compact .clean} +| argument {.whitespace-nowrap} | type {.whitespace-nowrap} | description | +| :---------------------------- | :-------------------------: | :-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | +| `-v` | toggle | Verbose output | +| `--aemb` | toggle | Output the estimated abundance value of contigs, the format of output file is: contig_id abundance_value | +| `--eqx` | toggle | Emit =/X instead of M CIGAR operations | +| `--no-PG` | toggle | Do not output PG header | +| `--details` | toggle | Add debugging details to SAM records | +| `--rg=` | [TAG:VALUE...] | Add read group metadata to SAM header (can be specified multiple times). Example: SM:samplename | +| `-N` | integer | Retain at most INT secondary alignments (is upper bounded by -M and depends on -S) [0] | +| `-m` | integer | Maximum seed length. Defaults to r - 50. For reasonable values on -l and -u, the seed length distribution is usually determined by parameters l and u. Then, this parameter is only active in regions where syncmers are very sparse. | +| `-k` | integer | Strobe length, has to be below 32. [20] | +| `-l` | integer | Lower syncmer offset from k/(k-s+1). Start sample second syncmer k/(k-s+1) + l syncmers downstream [0] | +| `-u` | integer | Upper syncmer offset from k/(k-s+1). End sample second syncmer k/(k-s+1) + u syncmers downstream [7] | +| `-c` | integer | Bitcount length between 2 and 63. [8] | +| `-s` | integer | Submer size used for creating syncmers [k-4]. Only even numbers on k-s allowed. A value of s=k-4 roughly represents w=10 as minimizer window [k-4]. It is recommended not to change this parameter unless you have a good understanding of syncmers as it will drastically change the memory usage and results with non default values. | +| `-b` | integer | No. of top bits of hash to use as bucket indices (8-31)[determined from reference size] | +| `-A` | integer | Matching score [2] | +| `-B` | integer | Mismatch penalty [8] | +| `-O` | integer | Gap open penalty [12] | +| `-E` | integer | Gap extension penalty [1] | +| `-L` | integer | Soft clipping penalty [10] | +| `-f` | float | Top fraction of repetitive strobemers to filter out from sampling [0.0002] | +| `-S` | float | Try candidate sites with mapping score at least S of maximum mapping score [0.5] | +| `-M` | integer | Maximum number of mapping sites to try [20] | +| `-R` | integer | Rescue level. Perform additional search for reads with many repetitive seeds filtered out. This search includes seeds of R*repetitive_seed_size_filter (default: R=2). Higher R than default makes strobealign significantly slower but more accurate. R <= 1 deactivates rescue and is the fastest. | + ++++ + +==- minimap2 ++++ :icon-git-merge: details +- ignores barcode information (but retains in output) +- ultra-fast +- highly tuned for long-read data + +[Minimap2](https://github.com/lh3/minimap2) is a versatile sequence alignment program that aligns DNA or mRNA sequences against a large reference database. +For ~10kb noisy reads sequences, minimap2 is tens of times faster than mainstream long-read mappers such as BLASR, BWA-MEM, NGMLR and GMAP. +The `BX:Z` tags in the read headers are still added to the alignment headers, even though barcodes +are not used to inform mapping. The `-m` threshold is used for alignment molecule assignment. + + ++++ :icon-code-square: minimap2 parameters +By default, Harpy runs `minimap2` with these parameters (excluding inputs and outputs): +```bash +minimap2 -t {threads} -a --MD -y -x map-{technology} -R "@RG\tID:samplename\tSM:samplename" {input.genome} {input.fastq} +``` + +Below is a list of all `minimap2` command line arguments, excluding those Harpy already uses or those made redundant by Harpy's implementation of it. +Values in `[brackets]` are the software's default. + +{.compact .clean} +| argument {.whitespace-nowrap} | category {.whitespace-nowrap} | description | +| :---------------------------- | :-----------------------------: | :------------------------------------------------------------------------------ | +| `-H` | indexing | use homopolymer-compressed k-mer (preferrable for PacBio) | +| `-k` | indexing | k-mer size (no larger than 28) [15] | +| `-w` | indexing | minimizer window size [10] | +| `-I` | indexing | split index for every ~NUM input bases [8G] | +| `-f` | mapping | filter out top FLOAT fraction of repetitive minimizers [0.0002] | +| `-g` | mapping | stop chain enlongation if there are no minimizers in INT-bp [5000] | +| `-G` | mapping | max intron length (effective with -xsplice; changing -r) [200k] | +| `-F` | mapping | max fragment length (effective with -xsr or in the fragment mode) [800] | +| `-r` | mapping | chaining/alignment bandwidth and long-join bandwidth [500,20000] | +| `-n` | mapping | minimal number of minimizers on a chain [3] | +| `-m` | mapping | minimal chaining score (matching bases minus log gap penalty) [40] | +| `-X` | mapping | skip self and dual mappings (for the all-vs-all mode) | +| `-p` | mapping | min secondary-to-primary score ratio [0.8] | +| `-N` | mapping | retain at most INT secondary alignments [5] | +| `-A` | alignment | matching score [2] | +| `-B` | alignment | mismatch penalty (larger value for lower divergence) [4] | +| `-O` | alignment | gap open penalty [4,24] | +| `-E` | alignment | gap extension penalty; a k-long gap costs min{O1+k*E1,O2+k*E2} [2,1] | +| `-z` | alignment | Z-drop score and inversion Z-drop score [400,200] | +| `-s` | alignment | minimal peak DP alignment score [80] | +| `-u` | alignment | how to find GT-AG. f:transcript strand, b:both strands, n:don't match GT-AG [n] | +| `-J` | alignment | splice mode. 0: original minimap2 model; 1: miniprot model [1] | +| `-j` | alignment | junctions in BED12 to extend *short* RNA-seq alignment [] | +| `-L` | file IO | write CIGAR with >65535 ops at the CG tag | +| `--cs[=STR]` | file IO | output the cs tag; STR is 'short' (if absent) or 'long' [none] | +| `--ds` | file IO | output the ds tag, which is an extension to cs | +| `--eqx` | file IO | write =/X CIGAR operators | +| `-Y` | file IO | use soft clipping for supplementary alignments | +| `-K` | file IO | minibatch size for mapping [500M] | ++++ + +=== \ No newline at end of file diff --git a/docs/Commands/align/bwa.md b/docs/Commands/align/bwa.md deleted file mode 100644 index 6abf48358..000000000 --- a/docs/Commands/align/bwa.md +++ /dev/null @@ -1,209 +0,0 @@ ---- -label: bwa -description: Align sequences with BWA-MEM2 -category: [linked-read, wgs] -tags: [linked-read, wgs] -icon: dot -order: 5 ---- - -# :icon-quote: align bwa - -=== :icon-checklist: You will need -- at least 4 cores/threads available -- a genome assembly in FASTA format: [!badge variant="success" text=".fasta"] [!badge variant="success" text=".fa"] [!badge variant="success" text=".fasta.gz"] [!badge variant="success" text=".fa.gz"] [!badge variant="secondary" text="case insensitive"] -- paired-end fastq sequence files [!badge variant="secondary" icon=":heart:" text="gzipped recommended"] - - **sample name**: [!badge variant="success" text="a-z"] [!badge variant="success" text="0-9"] [!badge variant="success" text="."] [!badge variant="success" text="_"] [!badge variant="success" text="-"] [!badge variant="secondary" text="case insensitive"] - - **forward**: [!badge variant="success" text="_F"] [!badge variant="success" text=".F"] [!badge variant="success" text=".1"] [!badge variant="success" text="_1"] [!badge variant="success" text="_R1_001"] [!badge variant="success" text=".R1_001"] [!badge variant="success" text="_R1"] [!badge variant="success" text=".R1"] - - **reverse**: [!badge variant="success" text="_R"] [!badge variant="success" text=".R"] [!badge variant="success" text=".2"] [!badge variant="success" text="_2"] [!badge variant="success" text="_R2_001"] [!badge variant="success" text=".R2_001"] [!badge variant="success" text="_R2"] [!badge variant="success" text=".R2"] - - **fastq extension**: [!badge variant="success" text=".fq"] [!badge variant="success" text=".fastq"] [!badge variant="secondary" text="case insensitive"] -=== - -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 [!badge corners="pill" text="harpy qc"](../qc.md). You can map reads onto a genome assembly with Harpy -using the [!badge corners="pill" text="align bwa"] module: - -```bash usage -harpy align bwa OPTIONS... REFERENCE INPUTS... -``` -```bash example -harpy align bwa genome.fasta Sequences/ -``` - -## :icon-terminal: Running Options -In addition to the [!badge variant="info" corners="pill" text="common runtime options"](/Getting_Started/common_options.md), the [!badge corners="pill" text="align bwa"] module is configured using these command-line arguments: - -{.compact .clean} -| argument {.whitespace-nowrap} | default {.whitespace-nowrap} | description | -| :------------------------------- | :--------------------------: | :---------------------------------------------------------------------------------------------------------------------------------------------- | -| `REFERENCE` | | [!badge variant="info" text="required"] Reference assembly for read mapping | -| `INPUTS` | | [!badge variant="info" text="required"] Files or directories containing [input FASTQ files](/Getting_Started/common_options.md#input-arguments) | -| `-depth-window` `-w` | `50000` | Interval size (in bp) for depth stats | -| `--extra-params` `-x` | | Additional BWA arguments, in quotes | -| `--keep-unmapped` `-u` | false | Output unmapped sequences too | -| `--molecule-distance` `-d` | `0` | Base-pair distance threshold to separate molecules given as base pairs, disabled with `0` | -| `--min-quality` `-q` | `30` | Minimum `MQ` (SAM mapping quality) to pass filtering | - -### 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 -to deconvolute alignments with the same barcode that might not have originated -from the same DNA molecule based on the [distance threshold](/Getting_Started/linked_read_data.md#barcode-thresholds) -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). - -==- What is the $MQ$ score? -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 -$$ - -![A visual explanation of MQ Score](/static/MQscore.png) -=== - -## Marking PCR duplicates -Harpy uses `samtools markdup` to mark putative PCR duplicates by using both the `BX` tag -as a UMI (unique molecule identified) for more accurate duplicate detection. 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**. - ----- - -## :icon-git-pull-request: BWA workflow -+++ :icon-git-merge: details -- ignores (but retains) barcode information -- fast - -The [BWA MEM](https://github.com/bwa-mem2/bwa-mem2) workflow maps all reads against the reference genome. Duplicates are marked using `samtools markdup`. -The `BX:Z` tags in the read headers are still added to the alignment headers, even though barcodes -are not used to inform mapping. The `-m` threshold is used for alignment molecule assignment. - -```mermaid -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:::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 - -``` -+++ :icon-file-directory: BWA output -The default output directory is `Align/bwa` with the folder structure below. `Sample1` is a generic sample name for demonstration purposes. -The resulting folder also includes a `workflow` directory (not shown) with workflow-relevant runtime files and information. -``` -Align/bwa -├── Sample1.bam -├── Sample1.bam.bai -├── logs -│ ├── sample1.bwa.log -│ ├── sample1.markdup.log -│ │── sample1.sort.log -└── reports - ├── barcodes.summary.html - ├── bwa.stats.html - ├── Sample1.html - └── data -    ├── bxstats - │ └── Sample1.bxstats.gz - └── coverage - └── Sample1.cov.gz -``` -{.compact} -| item {.whitespace-nowrap} | description | -| :------------------------------- | :------------------------------------------------------------------------------- | -| `*.bam` | sequence alignments for each sample | -| `*.bai` | sequence alignment indexes for each sample | -| `logs/*bwa.log` | output of BWA during run | -| `logs/*markdup.log` | stats provided by `samtools markdup` | -| `logs/*sort.log` | output of `samtools sort` | -| `reports/` | various counts/statistics/reports relating to sequence alignment | -| `reports/barcodes.summary.ipynb` | report summarizing barcode-specific metrics across all samples | -| `reports/bwa.summary.ipynb` | report summarizing `samtools stats` of raw alignments across all samples | -| `reports/Sample1.ipynb` | report summarizing BX tag metrics and alignment coverage | -| `reports/data/coverage/*.cov.gz` | output from mosdepth, used for reports | -| `reports/data/lrstats` | tabular data containing the information used to generate the BX stats in reports | - -+++ :icon-code-square: BWA parameters -By default, Harpy runs `bwa` with these parameters (excluding inputs and outputs): -```bash -bwa-mem2 mem -v 2 -T 10 -m 10 -C -R "@RG\tID:samplename\tSM:samplename" -``` - -Below is a list of all `bwa-mem2 mem` command line arguments, excluding those Harpy already uses or those made redundant by Harpy's implementation of BWA. -These are taken directly from the [BWA documentation](https://bio-bwa.sourceforge.net/bwa.shtml). -```bwa arguments - Algorithm options: - -k INT minimum seed length [19] - -w INT band width for banded alignment [100] - -d INT off-diagonal X-dropoff [100] - -r FLOAT look for internal seeds inside a seed longer than {-k} * FLOAT [1.5] - -y INT seed occurrence for the 3rd round seeding [20] - -c INT skip seeds with more than INT occurrences [500] - -D FLOAT drop chains shorter than FLOAT fraction of the longest overlapping chain [0.50] - -W INT discard a chain if seeded bases shorter than INT [0] - -S skip mate rescue - -P skip pairing; mate rescue performed unless -S also in use -Scoring options: - -A INT score for a sequence match, which scales options -TdBOELU unless overridden [1] - -B INT penalty for a mismatch [4] - -O INT[,INT] gap open penalties for deletions and insertions [6,6] - -E INT[,INT] gap extension penalty; a gap of size k cost '{-O} + {-E}*k' [1,1] - -L INT[,INT] penalty for 5'- and 3'-end clipping [5,5] - -U INT penalty for an unpaired read pair [17] -Input/output options: - -p smart pairing (ignoring in2.fq) - -H STR/FILE insert STR to header if it starts with @; or insert lines in FILE [null] - -j treat ALT contigs as part of the primary assembly (i.e. ignore .alt file) - -5 for split alignment, take the alignment with the smallest coordinate as primary - -q don't modify mapQ of supplementary alignments - -K INT process INT input bases in each batch regardless of nThreads (for reproducibility) [] - -h INT[,INT] if there are 80% of the max score, output all in XA [5,200] - -a output all alignments for SE or unpaired PE - -V output the reference FASTA header in the XR tag - -Y use soft clipping for supplementary alignments - -M mark shorter split hits as secondary - -I FLOAT[,FLOAT[,INT[,INT]]] - specify the mean, standard deviation (10% of the mean if absent), max - (4 sigma from the mean if absent) and min of the insert size distribution. - FR orientation only. [inferred] -``` - -+++ diff --git a/docs/Commands/align/standard.md b/docs/Commands/align/standard.md new file mode 100644 index 000000000..ad59fb6a9 --- /dev/null +++ b/docs/Commands/align/standard.md @@ -0,0 +1,336 @@ +--- +label: bwa/strobe/minimap +description: Align sequences with minibwa / strobealign /minimap2 +category: [linked-read, wgs] +tags: [linked-read, wgs] +icon: dot +order: 5 +--- + +# :icon-quote: align without linked-read information +The process and commands for using minibwa, strobealign, and minimap2 are identical and consolidated here. + +=== :icon-checklist: You will need +- at least 4 cores/threads available +- a genome assembly in FASTA format: [!badge variant="success" text=".fasta"] [!badge variant="success" text=".fa"] [!badge variant="success" text=".fasta.gz"] [!badge variant="success" text=".fa.gz"] [!badge variant="secondary" text="case insensitive"] +- paired-end fastq sequence files [!badge variant="secondary" icon=":heart:" text="gzipped recommended"] + - **sample name**: [!badge variant="success" text="a-z"] [!badge variant="success" text="0-9"] [!badge variant="success" text="."] [!badge variant="success" text="_"] [!badge variant="success" text="-"] [!badge variant="secondary" text="case insensitive"] + - **forward**: [!badge variant="success" text="_F"] [!badge variant="success" text=".F"] [!badge variant="success" text=".1"] [!badge variant="success" text="_1"] [!badge variant="success" text="_R1_001"] [!badge variant="success" text=".R1_001"] [!badge variant="success" text="_R1"] [!badge variant="success" text=".R1"] + - **reverse**: [!badge variant="success" text="_R"] [!badge variant="success" text=".R"] [!badge variant="success" text=".2"] [!badge variant="success" text="_2"] [!badge variant="success" text="_R2_001"] [!badge variant="success" text=".R2_001"] [!badge variant="success" text="_R2"] [!badge variant="success" text=".R2"] + - **fastq extension**: [!badge variant="success" text=".fq"] [!badge variant="success" text=".fastq"] [!badge variant="secondary" text="case insensitive"] +=== + +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 [!badge corners="pill" text="harpy qc"](../qc.md). You can map reads onto a genome assembly with Harpy +using the [!badge corners="pill" text="align bwa"] module: + +```bash usage +harpy align bwa|stobe|minimap OPTIONS... REFERENCE INPUTS... +``` +```bash example +harpy align bwa genome.fasta Sequences/ +``` + +## :icon-terminal: Running Options +In addition to the [!badge variant="info" corners="pill" text="common runtime options"](/Getting_Started/common_options.md), the [!badge corners="pill" text="align bwa"]/[!badge corners="pill" text="align strobe"] modules are configured using these command-line arguments: + +{.compact .clean} +| argument {.whitespace-nowrap} | default {.whitespace-nowrap} | description | +| :------------------------------- | :--------------------------: | :---------------------------------------------------------------------------------------------------------------------------------------------- | +| `REFERENCE` | | [!badge variant="info" text="required"] Reference assembly for read mapping | +| `INPUTS` | | [!badge variant="info" text="required"] Files or directories containing [input FASTQ files](/Getting_Started/common_options.md#input-arguments) | +| `-depth-window` `-w` | `50000` | Interval size (in bp) for depth stats | +| `--extra-params` `-x` | | Additional aligner-specific arguments, in quotes | +| `--keep-unmapped` `-u` | false | Output unmapped sequences too | +| `--molecule-distance` `-d` | `0` | Base-pair distance threshold to separate molecules given as base pairs, disabled with `0` | +| `--min-quality` `-q` | `30` | Minimum `MQ` (SAM mapping quality) to pass filtering | +| `--technology` `-t` | `sr` | [!badge variant="secondary" text="minimap only"] Sequence type preset [`sr`, `pb`, `hifi`, `ont`, `iclr`] | + +### Technology +This option is specific to [!badge corners="pill" text="align minimap"] and ignored otherwise. It configures the aligner for +specific data types, where: + +{.compact .clean} +| technology | description | +| :--------- | :----------------------------------------------------- | +| `sr` | short-read, typically Illumina, data (<600bp, default) | +| `pb` | PacBio data that is not HiFi | +| `hifi` | PacBio HiFi data | +| `ont` | Oxford Nanopore data | +| `iclr` | Illumina Complete Long Reads data | + +### 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 +to deconvolute alignments with the same barcode that might not have originated +from the same DNA molecule based on the [distance threshold](/Getting_Started/linked_read_data.md#barcode-thresholds) +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). + +==- What is the $MQ$ score? +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 +$$ + +![A visual explanation of MQ Score](/static/MQscore.png) +=== + +## Marking PCR duplicates +Harpy uses `samtools markdup` to mark putative PCR duplicates by using both the `BX` tag +as a UMI (unique molecule identified) for more accurate duplicate detection. 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**. + +---- + +## :icon-git-pull-request: Workflows +Regardless of which aligner you choose, the workflows are mostly identical: + ++++ :icon-git-pull-request: workflow +```mermaid +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 +``` + ++++ :icon-file-directory: output +The default output directory is `Align/{aligner}` with the folder structure below. +`Sample1` is a generic sample name for demonstration purposes. The resulting folder also includes a `workflow` directory +(not shown) with workflow-relevant runtime files and information. +``` +Align/{aligner} +├── Sample1.bam +├── Sample1.bam.bai +├── logs +│ ├── sample1.arachne.log +│ ├── sample1.markdup.log +│ │── sample1.sort.log +└── reports + ├── barcodes.summary.ipynb + ├── {aligner}.stats.ipynb + ├── Sample1.ipynb + └── data +    ├── lrstats + │ └── Sample1.lrstats.gz + └── coverage + ├── Sample1.molcov.gz + └── Sample1.cov.gz +``` +{.compact} +| item {.whitespace-nowrap} | description | +| :---------------------------------- | :------------------------------------------------------------------------------------- | +| `*.bam` | sequence alignments for each sample | +| `*.bai` | sequence alignment indexes for each sample | +| `logs/*{aligner}.log` | output of arachne during run | +| `logs/*markdup.log` | stats provided by `samtools markdup` _for invalid-barcoded reads_ | +| `logs/*sort.log` | output of `samtools sort` | +| `reports/` | various counts/statistics/reports relating to sequence alignment | +| `reports/barcodes.summary.ipynb` | report summarizing barcode-specific metrics across all samples | +| `reports/{aligner}.summary.ipynb` | report summarizing `samtools stats` of raw and processed alignments across all samples | +| `reports/Sample1.ipynb` | report summarizing BX tag metrics and alignment coverage | +| `reports/data/coverage/*.cov.gz` | output from mosdepth, used for reports | +| `reports/data/coverage/*.molcov.gz` | molecular coverage stats, used for reports | +| `reports/data/lrstats` | tabular data containing the information used to generate the BX stats in reports | ++++ + +==- minibwa ++++ :icon-git-merge: details +- ignores barcode information (but retains in output) +- fast and accurate + +The [minibwa](https://github.com/lh3/minibwa) workflow maps all reads against the reference genome. Duplicates are marked using `samtools markdup`. +The `BX:Z` tags in the read headers are still added to the alignment headers, even though barcodes +are not used to inform mapping. The `-m` threshold is used for calculating statistics. + ++++ :icon-code-square: minbwa parameters +By default, Harpy runs `minibwa` with these parameters (excluding inputs and outputs): +```bash +minibwa map -y -x sr -R "@RG\tID:samplename\tSM:samplename" +``` + +Below is a list of all `minibwa map` command line arguments, excluding those Harpy already uses or those made redundant by Harpy's implementation of BWA. + +{.compact .clean} +| argument {.whitespace-nowrap} | category {.whitespace-nowrap} | description | +| :---------------------------- | :-----------------------------: | :--------------------------------------------------------------------- | +| `-l` | common | treat reads 80% of the best hit, output them to XA [5] | +| `-Y` | file IO | use soft clipping for supplementary alignments | +| `-H` | file IO | if STR starts with @, insert to header; or insert lines in file STR [] | +| `-5` | file IO | take the alignment with the smallest query position as primary | +| `-K` | file IO | process NUM1-NUM2 bp of query sequences in a batch [100m,1g] | +| `--mmap[=lite]` | file IO | load the index via memory mapped files (slower mapping) [] | + ++++ + +==- strobealign ++++ :icon-git-merge: details +- ignores barcode information (but retains in output) +- ultra-fast +- [as-good-or-better accuracy](https://github.com/ksahlin/strobealign/blob/main/evaluation.md) to BWA MEM for sequences greater than 100bp + - accuracy may be lower for sequences shorter than 100bp + +The [strobealign](https://github.com/lksahlinh3/strobealign) workflow is nearly identical to the BWA workflow, +the only real difference being how the input genome is indexed and that alignment is performed with +`strobealign` instead of BWA. Duplicates are marked using `samtools markdup`. +The `BX:Z` tags in the read headers are still added to the alignment headers, even though barcodes +are not used to inform mapping. The `-m` threshold is used for alignment molecule assignment. + + ++++ :icon-code-square: strobealign parameters +By default, Harpy runs `strobealign` with these parameters (excluding inputs and outputs): +```bash +strobealign [--use-index -r ...] -t THREADS -U -C --rg-id={sample} --rg=SM:{sample} +``` + +Below is a list of all `strobealign` command line arguments, excluding those Harpy already uses or those made redundant by Harpy's implementation of it. + +{.compact .clean} +| argument {.whitespace-nowrap} | type {.whitespace-nowrap} | description | +| :---------------------------- | :-------------------------: | :-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | +| `-v` | toggle | Verbose output | +| `--aemb` | toggle | Output the estimated abundance value of contigs, the format of output file is: contig_id abundance_value | +| `--eqx` | toggle | Emit =/X instead of M CIGAR operations | +| `--no-PG` | toggle | Do not output PG header | +| `--details` | toggle | Add debugging details to SAM records | +| `--rg=` | [TAG:VALUE...] | Add read group metadata to SAM header (can be specified multiple times). Example: SM:samplename | +| `-N` | integer | Retain at most INT secondary alignments (is upper bounded by -M and depends on -S) [0] | +| `-m` | integer | Maximum seed length. Defaults to r - 50. For reasonable values on -l and -u, the seed length distribution is usually determined by parameters l and u. Then, this parameter is only active in regions where syncmers are very sparse. | +| `-k` | integer | Strobe length, has to be below 32. [20] | +| `-l` | integer | Lower syncmer offset from k/(k-s+1). Start sample second syncmer k/(k-s+1) + l syncmers downstream [0] | +| `-u` | integer | Upper syncmer offset from k/(k-s+1). End sample second syncmer k/(k-s+1) + u syncmers downstream [7] | +| `-c` | integer | Bitcount length between 2 and 63. [8] | +| `-s` | integer | Submer size used for creating syncmers [k-4]. Only even numbers on k-s allowed. A value of s=k-4 roughly represents w=10 as minimizer window [k-4]. It is recommended not to change this parameter unless you have a good understanding of syncmers as it will drastically change the memory usage and results with non default values. | +| `-b` | integer | No. of top bits of hash to use as bucket indices (8-31)[determined from reference size] | +| `-A` | integer | Matching score [2] | +| `-B` | integer | Mismatch penalty [8] | +| `-O` | integer | Gap open penalty [12] | +| `-E` | integer | Gap extension penalty [1] | +| `-L` | integer | Soft clipping penalty [10] | +| `-f` | float | Top fraction of repetitive strobemers to filter out from sampling [0.0002] | +| `-S` | float | Try candidate sites with mapping score at least S of maximum mapping score [0.5] | +| `-M` | integer | Maximum number of mapping sites to try [20] | +| `-R` | integer | Rescue level. Perform additional search for reads with many repetitive seeds filtered out. This search includes seeds of R*repetitive_seed_size_filter (default: R=2). Higher R than default makes strobealign significantly slower but more accurate. R <= 1 deactivates rescue and is the fastest. | + ++++ + +==- minimap2 ++++ :icon-git-merge: details +- ignores barcode information (but retains in output) +- ultra-fast +- highly tuned for long-read data + +[Minimap2](https://github.com/lh3/minimap2) is a versatile sequence alignment program that aligns DNA or mRNA sequences against a large reference database. +For ~10kb noisy reads sequences, minimap2 is tens of times faster than mainstream long-read mappers such as BLASR, BWA-MEM, NGMLR and GMAP. +The `BX:Z` tags in the read headers are still added to the alignment headers, even though barcodes +are not used to inform mapping. The `-m` threshold is used for alignment molecule assignment. + + ++++ :icon-code-square: minimap2 parameters +By default, Harpy runs `minimap2` with these parameters (excluding inputs and outputs): +```bash +minimap2 -t {threads} -a --MD -y -x map-{technology} -R "@RG\tID:samplename\tSM:samplename" +``` + +Below is a list of all `minimap2` command line arguments, excluding those Harpy already uses or those made redundant by Harpy's implementation of it. +Values in `[brackets]` are the software's default. + +{.compact .clean} +| argument {.whitespace-nowrap} | category {.whitespace-nowrap} | description | +| :---------------------------- | :-----------------------------: | :------------------------------------------------------------------------------ | +| `-H` | indexing | use homopolymer-compressed k-mer (preferrable for PacBio) | +| `-k` | indexing | k-mer size (no larger than 28) [15] | +| `-w` | indexing | minimizer window size [10] | +| `-I` | indexing | split index for every ~NUM input bases [8G] | +| `-f` | mapping | filter out top FLOAT fraction of repetitive minimizers [0.0002] | +| `-g` | mapping | stop chain enlongation if there are no minimizers in INT-bp [5000] | +| `-G` | mapping | max intron length (effective with -xsplice; changing -r) [200k] | +| `-F` | mapping | max fragment length (effective with -xsr or in the fragment mode) [800] | +| `-r` | mapping | chaining/alignment bandwidth and long-join bandwidth [500,20000] | +| `-n` | mapping | minimal number of minimizers on a chain [3] | +| `-m` | mapping | minimal chaining score (matching bases minus log gap penalty) [40] | +| `-X` | mapping | skip self and dual mappings (for the all-vs-all mode) | +| `-p` | mapping | min secondary-to-primary score ratio [0.8] | +| `-N` | mapping | retain at most INT secondary alignments [5] | +| `-A` | alignment | matching score [2] | +| `-B` | alignment | mismatch penalty (larger value for lower divergence) [4] | +| `-O` | alignment | gap open penalty [4,24] | +| `-E` | alignment | gap extension penalty; a k-long gap costs min{O1+k*E1,O2+k*E2} [2,1] | +| `-z` | alignment | Z-drop score and inversion Z-drop score [400,200] | +| `-s` | alignment | minimal peak DP alignment score [80] | +| `-u` | alignment | how to find GT-AG. f:transcript strand, b:both strands, n:don't match GT-AG [n] | +| `-J` | alignment | splice mode. 0: original minimap2 model; 1: miniprot model [1] | +| `-j` | alignment | junctions in BED12 to extend *short* RNA-seq alignment [] | +| `-L` | file IO | write CIGAR with >65535 ops at the CG tag | +| `--cs[=STR]` | file IO | output the cs tag; STR is 'short' (if absent) or 'long' [none] | +| `--ds` | file IO | output the ds tag, which is an extension to cs | +| `--eqx` | file IO | write =/X CIGAR operators | +| `-Y` | file IO | use soft clipping for supplementary alignments | +| `-K` | file IO | minibatch size for mapping [500M] | ++++ + +=== \ No newline at end of file diff --git a/docs/Commands/align/strobe.md b/docs/Commands/align/strobe.md deleted file mode 100644 index 4c8dd9fbb..000000000 --- a/docs/Commands/align/strobe.md +++ /dev/null @@ -1,198 +0,0 @@ ---- -label: strobe -description: Align sequences with strobealign -category: [linked-read, wgs] -tags: [linked-read, wgs] -icon: dot -order: 5 ---- - -# :icon-quote: align strobe - -=== :icon-checklist: You will need -- at least 4 cores/threads available -- a genome assembly in FASTA format: [!badge variant="success" text=".fasta"] [!badge variant="success" text=".fa"] [!badge variant="success" text=".fasta.gz"] [!badge variant="success" text=".fa.gz"] [!badge variant="secondary" text="case insensitive"] -- paired-end fastq sequence files [!badge variant="secondary" icon=":heart:" text="gzipped recommended"] - - **sample name**: [!badge variant="success" text="a-z"] [!badge variant="success" text="0-9"] [!badge variant="success" text="."] [!badge variant="success" text="_"] [!badge variant="success" text="-"] [!badge variant="secondary" text="case insensitive"] - - **forward**: [!badge variant="success" text="_F"] [!badge variant="success" text=".F"] [!badge variant="success" text=".1"] [!badge variant="success" text="_1"] [!badge variant="success" text="_R1_001"] [!badge variant="success" text=".R1_001"] [!badge variant="success" text="_R1"] [!badge variant="success" text=".R1"] - - **reverse**: [!badge variant="success" text="_R"] [!badge variant="success" text=".R"] [!badge variant="success" text=".2"] [!badge variant="success" text="_2"] [!badge variant="success" text="_R2_001"] [!badge variant="success" text=".R2_001"] [!badge variant="success" text="_R2"] [!badge variant="success" text=".R2"] - - **fastq extension**: [!badge variant="success" text=".fq"] [!badge variant="success" text=".fastq"] [!badge variant="secondary" text="case insensitive"] -=== - -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 [!badge corners="pill" text="harpy qc"](../qc.md). You can map reads onto a genome assembly with Harpy -using the [!badge corners="pill" text="align strobe"] module: - -```bash usage -harpy align strobe OPTIONS... REFERENCE INPUTS... -``` -```bash example -harpy align strobe genome.fasta Sequences/ -``` - -## :icon-terminal: Running Options -In addition to the [!badge variant="info" corners="pill" text="common runtime options"](/Getting_Started/common_options.md), the [!badge corners="pill" text="align strobe"] module is configured using these command-line arguments: - -{.compact .clean} -| argument {.whitespace-nowrap} | default {.whitespace-nowrap} | description | -| :-------------------------------- | :--------------------------: | :---------------------------------------------------------------------------------------------------------------------------------------------- | -| `REFERENCE` | | [!badge variant="info" text="required"] Reference genome for read mapping | -| `INPUTS` | | [!badge variant="info" text="required"] Files or directories containing [input FASTQ files](/Getting_Started/common_options.md#input-arguments) | -| `-depth-window` `-w` | `50000` | Interval size (in bp) for depth stats | -| `--extra-params` `-x` | | Additional stroebealign arguments, in quotes | -| `--keep-unmapped` `-u` | false | Output unmapped sequences too | -| `--min-quality` `-d` | `30` | Minimum `MQ` (SAM mapping quality) to pass filtering | -| `--molecule-distance` `-d` | `0` | Base-pair distance threshold to separate molecules given as base pairs, disabled with `0` | - -### Molecule distance -The `--molecule-distance` option is used during the alignment workflow -to deconvolute alignments with the same barcode that might not have originated -from the same DNA molecule based on the [distance threshold](/Getting_Started/linked_read_data.md#barcode-thresholds) -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). - -==- What is the $MQ$ score? -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 -$$ - -![A visual explanation of MQ Score](/static/MQscore.png) -=== - -## Marking PCR duplicates -Harpy uses `samtools markdup` to mark putative PCR duplicates by using both the `BX` tag -as a UMI (unique molecule identified) for more accurate duplicate detection. 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**. - ----- - -## :icon-git-pull-request: Strobealign workflow -+++ :icon-git-merge: details -- ignores (but retains) barcode information -- ultra-fast -- [as-good-or-better accuracy](https://github.com/ksahlin/strobealign/blob/main/evaluation.md) to BWA MEM for sequences greater than 100bp - - accuracy may be lower for sequences less than 100bp - -The [strobealign](https://github.com/lh3/strobealign) workflow is nearly identical to the BWA workflow, -the only real difference being how the input genome is indexed and that alignment is performed with -`strobealign` instead of BWA. Duplicates are marked using `samtools markdup`. -The `BX:Z` tags in the read headers are still added to the alignment headers, even though barcodes -are not used to inform mapping. The `-m` threshold is used for alignment molecule assignment. - -```mermaid -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 -``` -+++ :icon-file-directory: strobealign output -The default output directory is `Align/strobealign` with the folder structure below. `Sample1` is a generic sample name for demonstration purposes. -The resulting folder also includes a `workflow` directory (not shown) with workflow-relevant runtime files and information. -``` -Align/strobealign -├── Sample1.bam -├── Sample1.bam.bai -├── logs -│ ├── sample1.strobealign.log -│ ├── sample1.markdup.log -│ │── sample1.sort.log -└── reports - ├── barcodes.summary.html - ├── strobealign.stats.html - ├── Sample1.html - └── data -    ├── bxstats - │ └── Sample1.bxstats.gz - └── coverage - └── Sample1.cov.gz -``` -{.compact} -| item {.whitespace-nowrap} | description | -| :---------------------------------- | :------------------------------------------------------------------------------- | -| `*.bam` | sequence alignments for each sample | -| `*.bai` | sequence alignment indexes for each sample | -| `logs/*strobe.log` | output of strobealign during run | -| `logs/*markdup.log` | stats provided by `samtools markdup` | -| `logs/*sort.log` | output of `samtools sort` | -| `reports/` | various counts/statistics/reports relating to sequence alignment | -| `reports/barcodes.summary.ipynb` | report summarizing barcode-specific metrics across all samples | -| `reports/strobealign.summary.ipynb` | report summarizing `samtools stats` of raw alignments across all samples | -| `reports/Sample1.ipynb` | html report summarizing BX tag metrics and alignment coverage | -| `reports/data/coverage/*.cov.gz` | output from mosdepth, used for reports | -| `reports/data/lrstats` | tabular data containing the information used to generate the BX stats in reports | - -+++ :icon-code-square: strobealign parameters -By default, Harpy runs `strobealign` with these parameters (excluding inputs and outputs): -```bash -strobealign [--use-index -r ...] -t THREADS -U -C --rg-id={sample} --rg=SM:{sample} {input.genome} {input.fastq} -``` - -Below is a list of all `strobealign` command line arguments, excluding those Harpy already uses or those made redundant by Harpy's implementation of it. - -{.compact} -| argument {.whitespace-nowrap} | type {.whitespace-nowrap} | description | -| :---------------------------- | :-------------------------: | :-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | -| -v | toggle | Verbose output | -| --aemb | toggle | Output the estimated abundance value of contigs, the format of output file is: contig_id abundance_value | -| --eqx | toggle | Emit =/X instead of M CIGAR operations | -| --no-PG | toggle | Do not output PG header | -| --details | toggle | Add debugging details to SAM records | -| --rg= | [TAG:VALUE...] | Add read group metadata to SAM header (can be specified multiple times). Example: SM:samplename | -| -N | integer | Retain at most INT secondary alignments (is upper bounded by -M and depends on -S) [0] | -| -m | integer | Maximum seed length. Defaults to r - 50. For reasonable values on -l and -u, the seed length distribution is usually determined by parameters l and u. Then, this parameter is only active in regions where syncmers are very sparse. | -| -k | integer | Strobe length, has to be below 32. [20] | -| -l | integer | Lower syncmer offset from k/(k-s+1). Start sample second syncmer k/(k-s+1) + l syncmers downstream [0] | -| -u | integer | Upper syncmer offset from k/(k-s+1). End sample second syncmer k/(k-s+1) + u syncmers downstream [7] | -| -c | integer | Bitcount length between 2 and 63. [8] | -| -s | integer | Submer size used for creating syncmers [k-4]. Only even numbers on k-s allowed. A value of s=k-4 roughly represents w=10 as minimizer window [k-4]. It is recommended not to change this parameter unless you have a good understanding of syncmers as it will drastically change the memory usage and results with non default values. | -| -b | integer | No. of top bits of hash to use as bucket indices (8-31)[determined from reference size] | -| -A | integer | Matching score [2] | -| -B | integer | Mismatch penalty [8] | -| -O | integer | Gap open penalty [12] | -| -E | integer | Gap extension penalty [1] | -| -L | integer | Soft clipping penalty [10] | -| -f | float | Top fraction of repetitive strobemers to filter out from sampling [0.0002] | -| -S | float | Try candidate sites with mapping score at least S of maximum mapping score [0.5] | -| -M | integer | Maximum number of mapping sites to try [20] | -| -R | integer | Rescue level. Perform additional search for reads with many repetitive seeds filtered out. This search includes seeds of R*repetitive_seed_size_filter (default: R=2). Higher R than default makes strobealign significantly slower but more accurate. R <= 1 deactivates rescue and is the fastest. | - -+++ diff --git a/docs/Commands/phase/phase_snp.md b/docs/Commands/phase/phase_snp.md index d651b91b7..af79505c0 100644 --- a/docs/Commands/phase/phase_snp.md +++ b/docs/Commands/phase/phase_snp.md @@ -20,7 +20,7 @@ order: 6 You may want to phase your genotypes into haplotypes, as haplotypes tend to be more informative than unphased genotypes (higher polymorphism, captures relationship between genotypes). Phasing -genotypes into haplotypes requires alignment files, such as those produced by [!badge corners="pill" text="align bwa"](../Align/bwa.md) +genotypes into haplotypes requires alignment files, such as those produced by [!badge corners="pill" text="align bwa"](../Align/standard.md) and a variant call file, such as one produced by [!badge corners="pill" text="snp freebayes"](../snp.md) or [!badge corners="pill" text="impute"](../impute.md). **Phasing only works on SNP/indel data**, and will not work for structural variants produced by [!badge corners="pill" text="sv leviathan"](../SV/leviathan.md) diff --git a/docs/Commands/sv/leviathan.md b/docs/Commands/sv/leviathan.md index 3a4020958..5ee2833b0 100644 --- a/docs/Commands/sv/leviathan.md +++ b/docs/Commands/sv/leviathan.md @@ -43,7 +43,7 @@ from the sample names. A simple fix would be to use underscores (`_`) to differe !!! === -After reads have been aligned, _e.g._ with [!badge corners="pill" text="align bwa"](../Align/bwa.md), you can use those alignment files +After reads have been aligned, _e.g._ with [!badge corners="pill" text="align bwa"](../Align/standard.md), you can use those alignment files (`.bam`) to call structural variants in your data using LEVIATHAN. To make sure your data will work seemlessly with LEVIATHAN, the alignments in the [input BAM files](/Getting_Started/common_options.md) should **end** with a `BX:Z` tag. Use [!badge corners="pill" text="validate bam"](../validate.md) if you want to double-check file diff --git a/docs/Commands/sv/naibr.md b/docs/Commands/sv/naibr.md index e0a4b6416..120493f2b 100644 --- a/docs/Commands/sv/naibr.md +++ b/docs/Commands/sv/naibr.md @@ -40,7 +40,7 @@ from the sample names. A simple fix would be to use underscores (`_`) to differe !!! === -After reads have been aligned, _e.g._ with [!badge corners="pill" text="align bwa"](../Align/bwa.md), you can use those alignment files +After reads have been aligned, _e.g._ with [!badge corners="pill" text="align bwa"](../Align/standard.md), you can use those alignment files (`.bam`) to call structural variants in your data using NAIBR. While our testing shows that NAIBR tends to find known inversions that LEVIATHAN misses, the program requires **phased** bam files as input. That means the alignments have a `PS` or `HP` tag that indicate which haplotype the read/alignment belongs to. diff --git a/docs/Getting_Started/Guides/wgs_data.md b/docs/Getting_Started/Guides/wgs_data.md index 3c4c03f4a..81208f7bd 100644 --- a/docs/Getting_Started/Guides/wgs_data.md +++ b/docs/Getting_Started/Guides/wgs_data.md @@ -33,8 +33,8 @@ harpy qc --unlinked --trim-adapters auto --min-length 50 data/WGS/sample_*.gz ``` ## Sequence Alignment -Setting `--unlinked` disables linked-read specific routines in [!badge corners="pill" text="harpy align bwa"](/Commands/Align/bwa.md) - and [!badge corners="pill" text="harpy align strobe"](/Commands/Align/strobe.md). Doing so also ignores `--molecule-distance`. +Setting `--unlinked` disables linked-read specific routines in [!badge corners="pill" text="harpy align bwa"](/Commands/Align/standard.md) + and [!badge corners="pill" text="harpy align strobe"](/Commands/Align/standard.md). Doing so also ignores `--molecule-distance`. ```bash align example harpy align bwa --unlinked --min-quality 25 genome.fasta data/WGS/trimmed diff --git a/docs/Getting_Started/Resources/software.md b/docs/Getting_Started/Resources/software.md index fa4a8e75d..1119ef5d5 100644 --- a/docs/Getting_Started/Resources/software.md +++ b/docs/Getting_Started/Resources/software.md @@ -18,7 +18,6 @@ Issues with specific tools might warrant a discussion with the authors/developer | athena | [github](https://github.com/abishara/athena_meta), [publication](https://doi.org/10.1038/nbt.4266) | | bash | [website](https://www.gnu.org/software/bash/) | | bcftools | [github](https://github.com/samtools/bcftools), [website](https://samtools.github.io/bcftools/bcftools.html) | -| bwa | [github](https://github.com/lh3/bwa), [publication](http://arxiv.org/abs/1303.3997) | | conda | [github](https://github.com/conda) | | dmox | [gitlab](https://gitlab.mbb.cnrs.fr/ibonnici/dmox) | | fastp | [github](https://github.com/OpenGene/fastp), [publication](https://doi.org/10.1093/bioinformatics/bty560) | @@ -28,8 +27,8 @@ Issues with specific tools might warrant a discussion with the authors/developer | LEVIATHAN | [github](https://github.com/morispi/LEVIATHAN) ,[publication](https://doi.org/10.1101/2021.03.25.437002) | | links | [github](https://github.com/bcgsc/links), [publication](https://gigascience.biomedcentral.com/articles/10.1186/s13742-015-0076-3) | | LRez | [github](https://github.com/morispi/LRez), [publication](https://academic.oup.com/bioinformaticsadvances/article/1/1/vbab022/6375438?login=false) | -| mamba | [github](https://github.com/mamba-org/mamba) | -| Mimick | [github](https://github.com/pdimens/mimick), [VISOR/XENIA](https://github.com/davidebolo1993/VISOR/tree/master/VISOR/XENIA) | +| minibwa | [github](https://github.com/lh3/minibwa), [preprint](https://arxiv.org/abs/2606.15357) | +| minimap2 | [github](https://github.com/lh3/minimap2), [publication](https://doi.org/10.1093/bioinformatics/bty191) | | NAIBR | [github](https://github.com/raphael-group/NAIBR), [github (fork)](https://github.com/pontushojer/NAIBR) ,[publication](https://doi.org/10.1093/bioinformatics/btx712) | | papermill | [github](https://github.com/nteract/papermill) | | python | [website](https://www.python.org/) | @@ -42,7 +41,6 @@ Issues with specific tools might warrant a discussion with the authors/developer | strobealign | [github](https://github.com/ksahlin/strobealign), [publication](https://doi.org/10.1186/s13059-022-02831-7) | | tigmint | [github](https://github.com/bcgsc/tigmint), [publication](https://doi.org/10.1186/s12859-018-2425-6) | | whatshap | [github](https://github.com/whatshap/whatshap), [publication](https://doi.org/10.1101/085050) | -| xeus-python | [github](https://github.com/jupyter-xeus/xeus-python) | ## Software Packages {.compact .clean} @@ -51,11 +49,11 @@ Issues with specific tools might warrant a discussion with the authors/developer | AG-Grid | JavaScript | [github](https://github.com/ag-grid/ag-grid) | | altair | Python/JavaScript | [github](https://github.com/vega/altair) | | click | Python | [github](https://github.com/pallets/click) | -| pandas | Python | [github](https://github.com/pandas-dev/pandas) | -| vega-lite | JavaScript | [github](https://github.com/vega/vega-lite), [publication](https://doi.org/10.1109/tvcg.2016.2599030) | -| pysam | Python | [github](https://github.com/pysam-developers/pysam) | +| IPython | Python | [github](https://github.com/ipython/ipython) | | numpy | Python | [github](https://github.com/numpy/numpy) | +| pandas | Python | [github](https://github.com/pandas-dev/pandas) | | polars | Python/Rust | [github](https://github.com/pola-rs/polars) | -| IPython | Python/Rust | [github](https://github.com/ipython/ipython) | +| pysam | Python | [github](https://github.com/pysam-developers/pysam) | | rich | Python | [github](https://github.com/Textualize/rich) | | rich-click | Python | [github](https://github.com/ewels/rich-click) | +| vega-lite | JavaScript | [github](https://github.com/vega/vega-lite), [publication](https://doi.org/10.1109/tvcg.2016.2599030) | diff --git a/docs/Getting_Started/common_options.md b/docs/Getting_Started/common_options.md index 4c5f4bb47..2236ff489 100644 --- a/docs/Getting_Started/common_options.md +++ b/docs/Getting_Started/common_options.md @@ -56,7 +56,7 @@ in the modules' help strings and can be configured using these arguments: | `--contigs` | file path or list | | Contigs to plot in the report(s) | | `--help` | | | Show the module docstring | | `--hpc` `-H` | | | Have snakemake submit all jobs to an HPC ([details](Resources/hpc.md)) | -| `--output` `-O` | string | varies | Name of output directory | +| `--output` `-O` | string | varies | Name of output directory | | `--quiet` `-Q` | [0,1,2] | 0 | `0` prints all progress information, `1` prints unified progress bar, `2` suppressess all console output except errors | | `--setup` `-N` | toggle [!badge variant="secondary" corners="pill" text="hidden"] | | Perform validations and setup workflow environment, but don't run anything | | `--skip-reports` `-R` | toggle | | Skip the processing and generation of HTML reports in a workflow | @@ -99,7 +99,7 @@ exceed that number. !!! ### example -You could call [!badge corners="pill" text="align strobe"](/Commands/Align/strobe.md) and specify 20 threads with no output to console: +You could call [!badge corners="pill" text="align strobe"](/Commands/Align/standard.md) and specify 20 threads with no output to console: ```bash harpy align strobe --threads 20 --quiet 2 genome.fasta samples/trimmedreads @@ -119,10 +119,12 @@ and the contents therein also allow you to rerun the workflow manually. The `wor {.compact .clean} | item {.whitespace-nowrap} | contents | utility | | :------------------------- | :------------------------------------------------------------------------------------------------------------- | :----------------------------------------------------- | -| `workflow.smk` | Snakefile with the full recipe of the workflow | understanding the entire workflow | +| `workflow.smk` | Snakefile with the full recipe of the workflow[^1] | understanding the entire workflow | | `profile.yaml` | Configuration file for Snakemake workflow dispatching | general bookkeeping, advanced runs | | `workflow.yaml` | Configuration file generated from command-line arguments and consumed by the Snakefile | general bookkeeping, advanced runs | | `envs/` | Configurations of the software environments required by the workflow | bookkeeping | | `reference/` | Folder with a link or copy to the FASTA file used as the reference for various workflows | necessary for concurrent workflows to avoid data races | | `*.ipynb` | Jupyter notebook templates used to execute and generate reports | seeing math behind plots/tables or borrow code from | | `*.summary` | Plain-text overview of the important parts of the workflow | bookkeeping and writing Methods in manuscripts | + +[^1]: The `align` workflows are an exception, as they have two snakefiles: `workflow.smk` which contains the bulk of the workflow, and another with the aligner-specific rules named `align_{aligner}.smk` \ No newline at end of file diff --git a/docs/Getting_Started/inputformat.md b/docs/Getting_Started/inputformat.md index 8deb03934..79fcc720b 100644 --- a/docs/Getting_Started/inputformat.md +++ b/docs/Getting_Started/inputformat.md @@ -23,7 +23,7 @@ Compressed files are expected be compressed with either `gzip` or `bgzip` and en Unfortunately, there are many different ways of naming FASTQ files, which makes it difficult to accomodate every wacky iteration currently in circulation. While Harpy tries its best to be flexible, there are limitations. -To that end, for the [!badge corners="pill" text="preprocess"](/Commands/preprocess.md), [!badge corners="pill" text="qc"](/Commands/qc.md), and [!badge corners="pill" text="align"](/Commands/Align/bwa.md) modules, the +To that end, for the [!badge corners="pill" text="preprocess"](/Commands/preprocess.md), [!badge corners="pill" text="qc"](/Commands/qc.md), and [!badge corners="pill" text="align"](/Commands/Align/standard.md) modules, the most common FASTQ naming styles are supported: - **sample names**: [!badge variant="success" text="a-z"] [!badge variant="success" text="0-9"] [!badge variant="success" text="."] [!badge variant="success" text="_"] [!badge variant="success" text="-"] [!badge variant="secondary" text="case insensitive"] - you can mix and match special characters, but that's bad practice and not recommended diff --git a/docs/Getting_Started/reports.md b/docs/Getting_Started/reports.md index cdd505ab2..99ee0279a 100644 --- a/docs/Getting_Started/reports.md +++ b/docs/Getting_Started/reports.md @@ -13,7 +13,7 @@ were reached, Harpy 4.0 introduced a complete overhaul of the reporting system u Jupyter comes several benefits: - Code and output are stored in the notebook - GitHub, JupyterLab, and VScode (and derivatives) natively render notebooks nicely -- Harpy can leverage [MyST](https://mystmd.org/) (via Jupyter Book) to render everything into a _cohesive_ report webiste +- Harpy can leverage [MyST](https://mystmd.org/) (via Jupyter Book) to render everything into a _cohesive_ [report webiste](https://pdimens.github.io/GIH-experiments/) !!! Complete overhaul is not an overstatement-- all the R code was ported into Python and reformatted for Jupyter format. diff --git a/docs/retype.yml b/docs/retype.yml index fb19ece32..b84f8cde3 100644 --- a/docs/retype.yml +++ b/docs/retype.yml @@ -21,7 +21,7 @@ lastUpdated: enabled: true source: committer branding: - label: v4.1 + label: v4.2 logo: static/logo_trans.png logoDark: static/logo_trans.png logoAlign: left @@ -33,17 +33,13 @@ edit: base: / branch: "docs" links: -- text: Therkildsen Lab - link: https://therkildsen-lab.github.io/ - icon: https://therkildsen-lab.github.io/images/TherkildsenLogo.png +- text: BLink-seq + link: https://blinkseq.github.io/ + icon: https://blinkseq.github.io/_astro/logo.DGThwx8R_1NHPA4.webp target: blank -- text: Innovation Lab - link: https://www.genomicsinnovation.cornell.edu/ - icon: static/biotech.png - target: blank -- text: LASTQ Standard - link: https://pdimens.github.io/lastq/ - icon: https://pdimens.github.io/lastq/img/logo.png +- text: LASTQ Format + link: https://blinkseq.github.io/lastq/ + icon: https://blinkseq.github.io/lastq/img/logo.png target: blank - text: Source link: https://github.com/pdimens/harpy diff --git a/harpy/__main__.py b/harpy/__main__.py index bfe38b510..2fff256e8 100644 --- a/harpy/__main__.py +++ b/harpy/__main__.py @@ -86,4 +86,3 @@ def cli(): cli.add_command(validate.validate) cli.add_command(view.view) cli.add_command(template.template) - diff --git a/harpy/commands/align.py b/harpy/commands/align.py index 2101c30c4..ff375cf53 100644 --- a/harpy/commands/align.py +++ b/harpy/commands/align.py @@ -5,7 +5,7 @@ import rich_click as click from harpy.common.cli_filetypes import FASTAfile, FASTQfile, HPCProfile -from harpy.common.cli_params import BwaParams, SnakemakeParams, StrobeAlignParams +from harpy.common.cli_params import BwaParams, MinimapParams, SnakemakeParams, StrobeAlignParams from harpy.common.system_ops import container_ok from harpy.common.workflow import Workflow from harpy.validation.fasta import FASTA @@ -23,12 +23,12 @@ def align(): is carried over to the alignment records. """ -@click.command(no_args_is_help = True, context_settings={"allow_interspersed_args" : False}, epilog= "Documentation: https://pdimens.github.io/harpy/workflows/align/bwa/") +@click.command(no_args_is_help = True, context_settings={"allow_interspersed_args" : False}, epilog= "Documentation: https://pdimens.github.io/harpy/workflows/align/standard/") @click.option('-w', '--depth-window', panel = "Parameters", default = 50000, show_default = True, type = click.IntRange(min = 50), help = 'Interval size (in bp) for depth stats') @click.option('-x', '--extra-params', panel = "Parameters", type = BwaParams(), help = 'Additional bwa mem parameters, in quotes') @click.option('-u', '--keep-unmapped', panel = "Parameters", is_flag = True, default = False, help = 'Include unmapped sequences in output') @click.option('-q', '--min-quality', panel = "Parameters", default = 30, show_default = True, type = click.IntRange(0, 40, clamp = True), help = 'Minimum mapping quality to output') -@click.option('-d', '--molecule-distance', panel = "Parameters", default = 0, show_default = True, type = click.IntRange(min = 0), help = 'Distance cutoff for molecule assignment (bp)') +@click.option('-d', '--molecule-distance', panel = "Parameters", default = 50000, show_default = True, type = click.IntRange(min = 0), help = 'Distance cutoff for molecule assignment (bp)') @click.option('-O', '--output', panel = "Workflow Options", type = click.Path(exists = False, resolve_path = True), default = "Align/bwa", show_default=True, help = 'Output directory name') @click.option('-@', '--threads', panel = "Workflow Options", default = 4, show_default = True, type = click.IntRange(4,999, clamp = True), help = 'Number of threads to use') @click.option('-U','--unlinked', panel = "Parameters", is_flag = True, default = False, help = "Treat input data as not linked reads") @@ -45,16 +45,16 @@ def align(): @click.argument('inputs', required=True, type=FASTQfile(), nargs=-1) def bwa(reference, inputs, output, depth_window, unlinked, threads, keep_unmapped, extra_params, min_quality, molecule_distance, snakemake, skip_reports, quiet, hpc, clean, container, no_temp, setup): """ - Align sequences to reference genome using BWA MEM2 + Align sequences to reference genome using minibwa Provide the reference fasta followed by input fastq files and/or directories at the end of the command as individual files/folders, using shell wildcards (e.g. `data/echidna*.fastq.gz`), or both. - BWA is a fast, robust, and reliable aligner that does not use barcodes when mapping. + minibwa is the official successor to BWA, a fast, robust, and reliable aligner that does not use barcodes when mapping. Presence and type of linked-read data is auto-detected, but can be deliberately ignored using `-U`. Setting `--molecule-distance` to `>0` activates alignment-distance based barcode deconvolution for reporting only (the barcodes remain unmodified). """ - workflow = Workflow("align_bwa", "align_bwa.smk", output, container, clean, quiet) + workflow = Workflow("align_bwa", "align.smk", output, container, clean, quiet) workflow.setup_snakemake(threads, hpc, snakemake, no_temp) workflow.notebook_files = ["align_stats.ipynb", "align_lrstats.ipynb", "samtools_stats.ipynb"] workflow.conda = ["align", "qc"] @@ -85,12 +85,12 @@ def bwa(reference, inputs, output, depth_window, unlinked, threads, keep_unmappe workflow.initialize(setup) -@click.command(no_args_is_help = True, context_settings={"allow_interspersed_args" : False}, epilog= "Documentation: https://pdimens.github.io/harpy/workflows/align/strobe/") +@click.command(no_args_is_help = True, context_settings={"allow_interspersed_args" : False}, epilog= "Documentation: https://pdimens.github.io/harpy/workflows/align/standard/") @click.option('-w', '--depth-window', panel = "Parameters", default = 50000, show_default = True, type = click.IntRange(min = 50), help = 'Interval size (in bp) for depth stats') @click.option('-x', '--extra-params', panel = "Parameters", type = StrobeAlignParams(), help = 'Additional strobealign parameters, in quotes') @click.option('-u', '--keep-unmapped', panel = "Parameters", is_flag = True, default = False, help = 'Include unmapped sequences in output') @click.option('-q', '--min-quality', panel = "Parameters", default = 30, show_default = True, type = click.IntRange(0, 40, clamp = True), help = 'Minimum mapping quality to output') -@click.option('-d', '--molecule-distance', panel = "Parameters", default = 0, show_default = True, type = click.IntRange(min = 0), help = 'Distance cutoff for molecule assignment (bp)') +@click.option('-d', '--molecule-distance', panel = "Parameters", default = 50000, show_default = True, type = click.IntRange(min = 0), help = 'Distance cutoff for molecule assignment (bp)') @click.option('-O', '--output', panel = "Workflow Options", type = click.Path(exists = False, resolve_path = True), default = "Align/strobealign", show_default=True, help = 'Output directory name') @click.option('-@', '--threads', panel = "Workflow Options", default = 4, show_default = True, type = click.IntRange(4,999, clamp = True), help = 'Number of threads to use') @click.option('-U','--unlinked', panel = "Parameters", is_flag = True, default = False, help = "Treat input data as not linked reads") @@ -117,7 +117,7 @@ def strobe(reference, inputs, output, unlinked, keep_unmapped, depth_window, thr but can be deliberately ignored using `-U`. Setting `--molecule-distance` to `>0` activates alignment-distance based barcode deconvolution. """ - workflow = Workflow("align_strobe", "align_strobe.smk", output, container, clean, quiet) + workflow = Workflow("align_strobe", "align.smk", output, container, clean, quiet) workflow.setup_snakemake(threads, hpc, snakemake, no_temp) workflow.notebook_files = ["align_stats.ipynb", "align_lrstats.ipynb", "samtools_stats.ipynb"] workflow.conda = ["align", "qc"] @@ -148,5 +148,72 @@ def strobe(reference, inputs, output, unlinked, keep_unmapped, depth_window, thr workflow.initialize(setup) + +@click.command(no_args_is_help = True, context_settings={"allow_interspersed_args" : False}, epilog= "Documentation: https://pdimens.github.io/harpy/workflows/align/standard/") +@click.option('-w', '--depth-window', panel = "Parameters", default = 50000, show_default = True, type = click.IntRange(min = 50), help = 'Interval size (in bp) for depth stats') +@click.option('-x', '--extra-params', panel = "Parameters", type = MinimapParams(), help = 'Additional minimap2 parameters, in quotes') +@click.option('-u', '--keep-unmapped', panel = "Parameters", is_flag = True, default = False, help = 'Include unmapped sequences in output') +@click.option('-q', '--min-quality', panel = "Parameters", default = 30, show_default = True, type = click.IntRange(0, 40, clamp = True), help = 'Minimum mapping quality to output') +@click.option('-d', '--molecule-distance', panel = "Parameters", default = 50000, show_default = True, type = click.IntRange(min = 0), help = 'Distance cutoff for molecule assignment (bp)') +@click.option('-O', '--output', panel = "Workflow Options", type = click.Path(exists = False, resolve_path = True), default = "Align/minimap", show_default=True, help = 'Output directory name') +@click.option('-t', '--technology', panel = "Workflow Options", type = click.Choice(["sr", "pb", "hifi", "ont", "iclr"], case_sensitive=False), default = "sr", show_default=True, help = 'Sequence type (for minimap presets)') +@click.option('-@', '--threads', panel = "Workflow Options", default = 4, show_default = True, type = click.IntRange(4,999, clamp = True), help = 'Number of threads to use') +@click.option('-U','--unlinked', panel = "Parameters", is_flag = True, default = False, help = "Treat input data as not linked reads") +@click.option('--clean', hidden = True, panel = "Workflow Options", type = str, help = 'Delete the log (`l`), .snakemake (`s`), and/or workflow (`w`) folders when done') +@click.option('-C', '--container', panel = "Workflow Options", is_flag = True, default = False, help = 'Use a container instead of conda', callback=container_ok) +@click.option('-H', '--hpc', panel = "Workflow Options", type = HPCProfile(), help = 'HPC submission YAML configuration file') +@click.option('-Q', '--quiet', panel = "Workflow Options", default = 0, type = click.IntRange(0,2,clamp=True), help = '`0` all output, `1` progress bar, `2` no output') +@click.option('-T', '--no-temp', hidden = True, panel = "Workflow Options", is_flag = True, default = False, help = 'Don\'t delete temporary files') +@click.option('-N', '--setup', panel = "Workflow Options", is_flag = True, hidden = True, default = False, help = 'Setup the workflow and exit') +@click.option('-R', '--skip-reports', panel = "Workflow Options", is_flag = True, show_default = True, default = False, help = 'Don\'t generate HTML reports') +@click.option('-S', '--snakemake', panel = "Workflow Options", type = SnakemakeParams(), help = 'Additional Snakemake parameters, in quotes') +@click.help_option('--help', hidden = True) +@click.argument('reference', type=FASTAfile(), required = True, nargs = 1) +@click.argument('inputs', required=True, type=FASTQfile(), nargs=-1) +def minimap(reference, inputs, output, depth_window, unlinked, threads, keep_unmapped, extra_params, min_quality, technology, molecule_distance, snakemake, skip_reports, quiet, hpc, clean, container, no_temp, setup): + """ + Align sequences to reference genome using minimap2 + + Provide the reference fasta followed by input fastq files and/or directories at the end of the command as individual + files/folders, using shell wildcards (e.g. `data/echidna*.fastq.gz`), or both. + + Minimap2 is an ultra-fast aligner tuned for long reads (e.g., pacbio, nanopore) that does not use barcodes when mapping. + Presence and type of linked-read data is auto-detected, but can be deliberately ignored using `-U`. + Setting `--molecule-distance` to `>0` activates alignment-distance based barcode deconvolution for reporting only (the barcodes remain unmodified). + """ + workflow = Workflow("align_minimap", "align.smk", output, container, clean, quiet) + workflow.setup_snakemake(threads, hpc, snakemake, no_temp) + workflow.notebook_files = ["align_stats.ipynb", "align_lrstats.ipynb", "samtools_stats.ipynb"] + workflow.conda = ["align", "qc"] + + ## checks and validations ## + fastq = FASTQ(inputs, detect_bc = not unlinked, quiet = quiet) + fasta = FASTA(reference, quiet = quiet) + + workflow.linkedreads["type"] = fastq.lr_type + workflow.linkedreads["standardized"] = {"BX" : fastq.bx_tag, "VX": fastq.vx_tag} + workflow.notebooks["skip"] = skip_reports + workflow.input(fasta.file, "reference") + workflow.input(fastq.files, "fastq") + workflow.param(technology.lower(), "aligner-technology") + workflow.param(fastq.illumina_old, "illumina-format-old") + workflow.param(molecule_distance, "distance-threshold") + workflow.param(min_quality, "min-map-quality") + workflow.param(keep_unmapped, "keep-unmapped") + workflow.param(depth_window, "depth-windowsize") + if extra_params: + workflow.param(extra_params, "extra") + + workflow.info = { + "Samples": fastq.count, + "Linked-Read Type": fastq.lr_type, + "Technology": technology, + "Reference": os.path.basename(reference), + "Output Folder" : os.path.relpath(output) + "/" + } + + workflow.initialize(setup) + align.add_command(bwa) align.add_command(strobe) +align.add_command(minimap) diff --git a/harpy/commands/diagnose.py b/harpy/commands/diagnose.py index 8bc83456d..2a5d14c61 100644 --- a/harpy/commands/diagnose.py +++ b/harpy/commands/diagnose.py @@ -146,7 +146,7 @@ def rule(directory): if failed_rule: hp.log(f"Failing rule: [yellow]{failed_rule.lstrip()}", newline=True) else: - hp.log(f"No errors found in {os.path.basename(latest_log)}", style = "green", markup=False, highlight=False) + hp.log(f"No errors found in {os.path.basename(latest_log)}") sys.exit(0) if infiles: if not os.path.exists(CONFIG_FILE): diff --git a/harpy/commands/impute.py b/harpy/commands/impute.py index d6b1c8472..2d3740602 100644 --- a/harpy/commands/impute.py +++ b/harpy/commands/impute.py @@ -62,7 +62,7 @@ def impute(parameters, vcf, inputs, output, strategy, buffer, grid_size, threads ## checks and validations ## params = ImputeParams(parameters, quiet) alignments = XAM(inputs, quiet = quiet) - vcffile = VCF(vcf, workflow.workflow_directory, quiet) + vcffile = VCF(vcf, workflow.workflow_directory, quiet, threads = threads) vcffile.find_biallelic_contigs() vcffile.match_samples(alignments.files, vcf_samples) if vcf_samples: diff --git a/harpy/commands/phase.py b/harpy/commands/phase.py index f577d13f0..7a50341a5 100644 --- a/harpy/commands/phase.py +++ b/harpy/commands/phase.py @@ -57,7 +57,7 @@ def bam(vcf, inputs, output, threads, unlinked, vcf_samples, molecule_distance, ## checks and validations ## alignments = XAM(inputs, detect_bc= not unlinked, quiet = quiet) - vcffile = VCF(vcf, workflow.workflow_directory, quiet = quiet) + vcffile = VCF(vcf, workflow.workflow_directory, quiet = quiet, threads = threads) vcffile.check_phase() vcffile.match_samples(alignments.files, vcf_samples) if vcf_samples: @@ -126,7 +126,7 @@ def snp(vcf, inputs, output, threads, unlinked, min_map_quality, min_base_qualit ## checks and validations ## alignments = XAM(inputs, detect_bc= not unlinked, quiet = quiet) - vcffile = VCF(vcf, workflow.workflow_directory, quiet = quiet) + vcffile = VCF(vcf, workflow.workflow_directory, quiet = quiet, threads = threads) vcffile.match_samples(alignments.files, vcf_samples) if contigs: vcffile.match_contigs(contigs) diff --git a/harpy/commands/preprocess.py b/harpy/commands/preprocess.py index 50986ad06..13f7afb17 100644 --- a/harpy/commands/preprocess.py +++ b/harpy/commands/preprocess.py @@ -11,8 +11,6 @@ from harpy.common.workflow import Workflow from harpy.validation.fastq import FASTQ - -#TODO update docs link to dedicated pages @click.group(options_metavar='') @click.help_option('--help', hidden = True) def preprocess(): diff --git a/harpy/commands/report.py b/harpy/commands/report.py index a925f1639..0448324a2 100644 --- a/harpy/commands/report.py +++ b/harpy/commands/report.py @@ -9,10 +9,10 @@ from rich.live import Live from rich.panel import Panel + from harpy.common.printing import HarpyPrint from harpy.report.render import ReportRender from harpy.report.static import ReportStatic -from harpy.report.utilities import check_tool from harpy.common.cli_filetypes import IPYNBfile @click.group(options_metavar='') @@ -130,28 +130,13 @@ def static(notebooks, debug, self_contained): created in the same directories as their source `.ipynb` files, but will lack the nicer features and formatting of a proper MyST-MD website. """ - check_tool( - "jupyter", - "It can be installed with using one of these methods:\npip: [green]pip install -U nbconvert[/]\n" - "conda: [green]conda install -c conda-forge nbconvert[/]\n" - "pixi: [green]pixi add nbconvert[/]" - ) - if self_contained: - check_tool( - "monolith", - "Monolith is required to flatten an HTML notebook and bundle the Javascript and CSS elements within it. " - "Harpy does not provide it, but it can be installed using " - "cargo ([green]cargo install monolith[/]) or by downloading and adding " - "a pre-built binary from [blue]https://github.com/Y2Z/monolith/releases[/] to your PATH." - ) all_notebooks = [nb for group in notebooks for nb in group] + rs = ReportStatic(quiet = not debug, static = self_contained) n = len(all_notebooks) if n > 1 : print(f"Converting {n} notebooks into HTML files.", file = sys.stderr) - rs = ReportStatic("", quiet = not debug, static = self_contained) for nb in all_notebooks: - rs.notebook = nb - rs.convert() + rs.convert(nb) report.add_command(live) report.add_command(static) \ No newline at end of file diff --git a/harpy/common/cli_params.py b/harpy/common/cli_params.py index 57856fafa..45f4f4ace 100644 --- a/harpy/common/cli_params.py +++ b/harpy/common/cli_params.py @@ -34,19 +34,19 @@ class BwaParams(click.ParamType): """A class for a click type that validates bwa extra-params.""" name = "bwa_params" def convert(self, value, param, ctx): - harpy_options = "-C -v -t -R -T -m".split() - valid_options = "-k -w -d -r -y -c -D -W -S -P -A -B -O -E -L -U -p -R -H -j -5 -q -K -h -a -V -Y -M -I".split() + harpy_options = "-y -x -R".split() + valid_options = "-l -b --hic --meth -k -c -g -w -W -m -p -N --chain-only -A -B -O -E -s -P --rescue -I --outn --outs --xa -Y -H -5 -K --mmap".split() opts = 0 - docs = "https://github.com/bwa-mem2/bwa-mem2" + docs = "https://github.com/lh3/minibwa" for i in shellsplit(value): if i.startswith("-"): opts += 1 if i in harpy_options: - self.fail(f"{i} is already used by Harpy when calling bwa-mem2 mem.", param, ctx) + self.fail(f"{i} is already used by Harpy when calling minibwa.", param, ctx) if i not in valid_options: - self.fail(f"{i} is not a valid bwa-mem2 option. See the bwa documentation for a list of available options: {docs}.", param, ctx) + self.fail(f"{i} is not a valid minibwa option. See the minibwa documentation for a list of available options: {docs}.", param, ctx) if opts < 1: - self.fail(f"No valid options recognized. Available bwa-mem2 options begin with one dash (e.g. -M). See the bwa-mem2 documentation: {docs}.", param, ctx) + self.fail(f"No valid options recognized. See the minibwa documentation: {docs}. Available minibwa options are: {' '.join(valid_options)}", param, ctx) return sanitize_shell(value) class StrobeAlignParams(click.ParamType): @@ -68,6 +68,25 @@ def convert(self, value, param, ctx): self.fail(f"No valid options recognized. Available strobealign options begin with one or two dashes (e.g. --eqx or -L). See the strobealign documentation for a list of available options: {docs}.", param, ctx) return sanitize_shell(value) +class MinimapParams(click.ParamType): + """A class for a click type that validates minimap extra-params.""" + name = "minimap_params" + def convert(self, value, param, ctx): + harpy_options = "-y --MD -y -a -x -ax".split() + valid_options = "--secondary -H -k -w -I -f -g -G -F -r -n -m -X -p -N -A -B -O -E -z -s -u -J -j -L --cs --ds --eqx -Y -K".split() + opts = 0 + docs = "https://github.com/lh3/minimap2" + for i in shellsplit(value): + if i.startswith("-"): + opts += 1 + if i in harpy_options: + self.fail(f"{i} is already used by Harpy when calling minimap2.", param, ctx) + if i not in valid_options: + self.fail(f"{i} is not a valid minimap2 option. See the minimap2 documentation for a list of available options: {docs}.", param, ctx) + if opts < 1: + self.fail(f"No valid options recognized. See the minimap2 documentation: {docs}. Available minimap2 options are: {' '.join(valid_options)}", param, ctx) + return sanitize_shell(value) + class SpadesParams(click.ParamType): """A class for a click type that validates spades extra-params.""" name = "spades_params" diff --git a/harpy/common/environments.py b/harpy/common/environments.py index 22b42dfde..92ef19802 100644 --- a/harpy/common/environments.py +++ b/harpy/common/environments.py @@ -20,8 +20,10 @@ class HarpyEnvs(): def __init__(self): self.__environments__: dict = { "align" : [ - "bioconda::bwa-mem2", + #"bioconda::arachne", "bioconda::bwa", + "bioconda::minibwa", + "bioconda::minimap2", "bioconda::samtools=1.23", "bioconda::seqtk", "bioconda::strobealign", diff --git a/harpy/common/summaries.py b/harpy/common/summaries.py index 721d372de..b1c5da4e0 100644 --- a/harpy/common/summaries.py +++ b/harpy/common/summaries.py @@ -29,18 +29,18 @@ def align_bwa(self): unmapped = "" if keep_unmapped else "-F 4" bx_mode = "--barcode-tag BX" if not ignore_bx else "" - bwa_static = "-C -v 2" if is_standardized else "-v 2" + bwa_static = "-y -x sr" if is_standardized else "-x sr" extra = extra - align = "Sequences were aligned with BWA using:\n" - align += f'\tbwa mem {bwa_static} {extra} -R "@RG\\tID:SAMPLE\\tSM:SAMPLE" genome forward_reads reverse_reads |\n' + align = "Sequences were aligned with minibwa using:\n" + align += f'\tminibwa map {bwa_static} {extra} -R "@RG\\tID:SAMPLE\\tSM:SAMPLE" genome forward_reads reverse_reads |\n' align += f"\tsamtools view -h {unmapped} -q {quality}" duplicates = "Duplicates in the alignments were marked following:\n" duplicates += "\tsamtools collate |\n" duplicates += "\tsamtools fixmate |\n" - duplicates += f"\tsamtools sort -T SAMPLE -m 2000M |\n" + duplicates += "\tsamtools sort -T SAMPLE -m 2000M |\n" duplicates += f"\tsamtools markdup -S {bx_mode} -d 100 (2500 for novaseq)" - standardization = "Barcodes were standardized to BX + VX format in the aligments using:\n" + standardization = "If linked reads, barcodes were standardized to BX + VX format in the aligments using:\n" standardization += "\tdjinn-standardize {input.bam} > {output.bam}" self.summary.append("The harpy align bwa workflow ran using these parameters:") self.summary.append(f"The provided genome: {genomefile}") @@ -74,7 +74,7 @@ def align_strobe(self): duplicates += "\tsamtools fixmate |\n" duplicates += f"\tsamtools sort -T SAMPLE --reference {genomefile} -m 2000M |\n" duplicates += f"\tsamtools markdup -S {bx_mode} -d 100 (2500 for novaseq)" - standardization = "Barcodes were standardized in the aligments using:\n" + standardization = "If linked reads, barcodes were standardized in the aligments using:\n" standardization += "\tstandardize-barcodes-sam > {output} < {input}" self.summary.append("The harpy align strobe workflow ran using these parameters:") self.summary.append(f"The provided genome: {genomefile}") @@ -83,6 +83,40 @@ def align_strobe(self): self.summary.append(standardization) self.summary.append(duplicates) + def align_minimap(self): + ignore_bx = self.WORKFLOW.get("linkedreads", {}).get("type", 'none') == "none" + bx_tag = self.WORKFLOW.get("linkedreads", {}).get("standardized", {}).get("BX", False) + vx_tag = self.WORKFLOW.get("linkedreads", {}).get("standardized", {}).get("VX", False) + tech = self.PARAMETERS.get("aligner-technology", "sr") + is_standardized = bx_tag and vx_tag + keep_unmapped = self.PARAMETERS.get("keep-unmapped", False) + extra = self.PARAMETERS.get("extra", "") + genomefile = self.INPUTS["reference"] + quality = self.PARAMETERS.get("min-map-quality", 30) + + unmapped = "" if keep_unmapped else "-F 4" + bx_mode = "--barcode-tag BX" if not ignore_bx else "" + tech = f"-ax map-{tech}" if tech != "sr" else "-ax map sr" + static = "-y --MD" if is_standardized else "--MD" + extra = self.PARAMETERS.get("extra", "") + + align = "Sequences were aligned with minimap2 using:\n" + align += f"\tminimap2 {tech} {static} -R \"@RG\\tID:SAMPLE\\tSM:SAMPLE\" {extra} genome reads.F.fq reads.R.fq |\n" + align += f"\t\tsamtools view -h {unmapped} -q {quality}" + duplicates = "Duplicates in the alignments were marked following:\n" + duplicates += "\tsamtools collate |\n" + duplicates += "\tsamtools fixmate |\n" + duplicates += f"\tsamtools sort -T SAMPLE --reference {genomefile} -m 2000M |\n" + duplicates += f"\tsamtools markdup -S {bx_mode} -d 100 (2500 for novaseq)" + standardization = "If linked reads, barcodes were standardized in the aligments using:\n" + standardization += "\tstandardize-barcodes-sam > {output} < {input}" + self.summary.append("The harpy align minimap workflow ran using these parameters:") + self.summary.append(f"The provided genome: {genomefile}") + self.summary.append(align) + if not ignore_bx: + self.summary.append(standardization) + self.summary.append(duplicates) + def assembly(self): # SPADES max_mem = self.PARAMETERS.get("spades", {}).get("max-memory", 'auto') diff --git a/harpy/common/workflow.py b/harpy/common/workflow.py index 82e83d1e4..b72cbf00b 100644 --- a/harpy/common/workflow.py +++ b/harpy/common/workflow.py @@ -188,6 +188,25 @@ def fetch_snakefile(self): f"The required snakefile [blue bold]{self.snakefile}[/] was not found in the Harpy installation.", "There may be an issue with your Harpy installation, which would require reinstalling Harpy. Alternatively, there may be an issue with your conda/mamba environment or configuration." ) + # get complimentary snakefiles for alignment workflows + if self.snakefile == "align.smk": + # self.name should be 'align_bwa' and such + dest_file = os.path.join(self.workflow_directory, f"{self.name}.smk") + source_file = resources.files("harpy.snakefiles") / f"{self.name}.smk" + try: + with ( + resources.as_file(source_file) as _source, + open(_source, 'r') as smk_in, + open(dest_file, 'w') as smk_out + ): + smk_out.write(f"## SOURCE: Harpy version {self.version}\n" + smk_in.read()) + #shutil.copy2(_source, dest_file) + except (FileNotFoundError, KeyError): + self.print.error( + "snakefile missing", + f"The required snakefile [blue bold]{self.snakefile}[/] was not found in the Harpy installation.", + "There may be an issue with your Harpy installation, which would require reinstalling Harpy. Alternatively, there may be an issue with your conda/mamba environment or configuration." + ) def fetch_scripts(self) -> None: """ diff --git a/harpy/notebooks/align_lrstats.ipynb b/harpy/notebooks/align_lrstats.ipynb index ca660450c..0f404353e 100644 --- a/harpy/notebooks/align_lrstats.ipynb +++ b/harpy/notebooks/align_lrstats.ipynb @@ -33,8 +33,7 @@ "from harpy.report.utilities import StopExecution\n", "import os\n", "from pathlib import Path\n", - "import polars as pl\n", - "_ = pl.Config.set_tbl_hide_column_data_types()\n" + "import polars as pl\n" ] }, { @@ -78,13 +77,13 @@ " if tot_reads == 0:\n", " return {\n", " \"Sample\": samplename,\n", - " \"Total Reads\": 0,\n", + " \"Fragments\": 0,\n", " \"Valid\": 0,\n", " \"Linked\": 0,\n", " \"Barcodes\": 0,\n", " \"Unique Molecules\": 0,\n", " \"Molecules/Barcode\": 0,\n", - " \"Reads/Molecule\": 0,\n", + " \"Fragments/Molecule\": 0,\n", " \"Insert Molecule Coverage\": 0,\n", " \"Exact Molecule Coverage\": 0,\n", " \"N50\": 0,\n", @@ -109,13 +108,13 @@ "\n", " return {\n", " \"Sample\": samplename,\n", - " \"Total Reads\": tot_reads,\n", + " \"Fragments\": tot_reads,\n", " \"Valid\": round(tot_valid_reads / tot_reads, 4),\n", " \"Linked\": round(tot_linked_reads / tot_reads, 4),\n", " \"Barcodes\": tot_uniq_bx,\n", " \"Unique Molecules\": tot_mol,\n", " \"Molecules/Barcode\": round(tot_mol / tot_uniq_bx, 2),\n", - " \"Reads/Molecule\": avg_reads_per_mol,\n", + " \"Fragments/Molecule\": avg_reads_per_mol,\n", " \"Insert Molecule Coverage\": avg_mol_cov,\n", " \"Exact Molecule Coverage\": avg_mol_cov_bp,\n", " \"N50\": n50,\n", @@ -139,13 +138,13 @@ "aggregate_df = (\n", " pl.DataFrame(records, schema={\n", " \"Sample\": pl.String,\n", - " \"Total Reads\": pl.Int64,\n", + " \"Fragments\": pl.Int64,\n", " \"Valid\": pl.Float64,\n", " \"Linked\": pl.Float64,\n", " \"Barcodes\": pl.Int64,\n", " \"Unique Molecules\": pl.Int64,\n", " \"Molecules/Barcode\": pl.Float64,\n", - " \"Reads/Molecule\": pl.Float64,\n", + " \"Fragments/Molecule\": pl.Float64,\n", " \"Insert Molecule Coverage\": pl.Float64,\n", " \"Exact Molecule Coverage\": pl.Float64,\n", " \"N50\": pl.Int64,\n", @@ -199,15 +198,15 @@ "_bpmsd = round(aggregate_df['Molecules/Barcode'].replace(0, None).std() or 0, 2)\n", "_linked = round(aggregate_df['Linked'].mean() or 0, 2)\n", "_linkedsd = round(aggregate_df['Linked'].std() or 0, 2)\n", - "_rpm = round(aggregate_df['Reads/Molecule'].replace(0, None).mean() or 0, 2)\n", - "_rpmsd = round(aggregate_df['Reads/Molecule'].replace(0, None).std() or 0, 2)\n", + "_rpm = round(aggregate_df['Fragments/Molecule'].replace(0, None).mean() or 0, 2)\n", + "_rpmsd = round(aggregate_df['Fragments/Molecule'].replace(0, None).std() or 0, 2)\n", "_valid = round(aggregate_df['Valid'].replace(0, None).mean() or 0, 2)\n", "\n", "(\n", " StatsBox()\n", " .add(len(aggregate_df), \"Samples\")\n", " .add(_bpm, \"Molecules/Barcode\", plus_minus=_bpmsd)\n", - " .conditional(_rpm, \"Reads/Mol\", 2, plus_minus=_rpmsd)\n", + " .conditional(_rpm, \"Fragments/Mol\", 2, plus_minus=_rpmsd)\n", " .conditional(_linked, \"Avg Linked\", 0.4, as_percent=True)\n", " .conditional(_valid, \"Valid BX\", 0.4, as_percent=True)\n", " .add(_n50, \"N50\", units = \"kb\")\n", @@ -285,24 +284,28 @@ "metadata": {}, "source": [ "This table [^1] is an aggregation of data for each sample based on their `*.lrstats.gz` file.\n", - "Every column after `Barcodes` ignores singletons in its calculations.\n", + "Every column after `Barcodes` ignores singletons in its calculations. Note the term **fragment**[^2] instead\n", + "of \"reads\". See [supporting info](#supporting-info) for a thorough explanation.\n", "\n", "[^1]:\n", - " | Column | Description |\n", - " |:-------|:------------|\n", - " | `Sample` | name of the sample |\n", - " | `Total Reads`| total number of alignments |\n", - " | `Valid` | proportion of valid barcoded alignments |\n", - " | `Unique Molecules` | the unique DNA molecules as inferred from linked-read barcodes |\n", - " | `Linked` | molecules composed of two or more single/paired-end sequences, in other words, molecules with linked-read information |\n", - " | `Barcodes` | number of unique barcodes, which may differ from unique molecules after deconvolution |\n", - " | `Molecules/Barcode` | molecule-to-barcode ratio, which helps benchmark deconvolution performance, if performed |\n", - " | `Reads/Molecule` | average number of reads per unique molecule |\n", - " | `Insert Molecule Coverage` | average percent of a molecule that is covered by a read, where coverage includes unsequenced gaps between linked reads |\n", - " | `Exact Molecule Coverage` | average percent molecule coverage, where coverage only includes sequences and not the gaps between linked reads |\n", - " | `N50` | N50 of inferred molecules |\n", - " | `N75` | N75 of inferred molecules |\n", - " | `N90` | N90 of inferred molecules |" + " | Column | Description |\n", + " | :------------------------- | :--------------------------------------------------------------------------------------------------------------------------- |\n", + " | `Sample` | name of the sample |\n", + " | `Fragments` | total number of \"fragments\", which are a DNA fragment from a molecule that appears as either a paired-end or single-end read |\n", + " | `Valid` | proportion of valid barcoded alignments |\n", + " | `Unique Molecules` | the unique DNA molecules as inferred from linked-read barcodes |\n", + " | `Linked` | molecules composed of two or more single/paired-end sequences, in other words, molecules with linked-read information |\n", + " | `Barcodes` | number of unique barcodes, which may differ from unique molecules after deconvolution |\n", + " | `Molecules/Barcode` | molecule-to-barcode ratio, which helps benchmark deconvolution performance, if performed |\n", + " | `Fragments/Molecule` | average number of reads per unique molecule |\n", + " | `Insert Molecule Coverage` | average percent of a molecule that is covered by a read, where coverage includes unsequenced gaps between linked reads |\n", + " | `Exact Molecule Coverage` | average percent molecule coverage, where coverage only includes sequences and not the gaps between linked reads |\n", + " | `N50` | N50 of inferred molecules |\n", + " | `N75` | N75 of inferred molecules |\n", + " | `N90` | N90 of inferred molecules |\n", + "\n", + "[^2]:\n", + " The term 'fragment' more accurately captures how many representative samples a source molecule has by not double-counting mates of read pairs" ] }, { @@ -390,24 +393,24 @@ }, "outputs": [], "source": [ - "_hist = binned_histogram(aggregate_df['Reads/Molecule'], 0.25, True)\n", + "_hist = binned_histogram(aggregate_df['Fragments/Molecule'], 0.25, True)\n", "with SafeRender(_hist):\n", " (\n", " alt.Chart(_hist)\n", " .mark_area(interpolate = \"monotone\")\n", " .encode(\n", " x=alt.X('interval:O')\n", - " .axis(title='Reads Per Molecule', tickMinStep=0.5, labelAngle=-40)\n", + " .axis(title='Fragments Per Molecule', tickMinStep=0.5, labelAngle=-40)\n", " .scale(domainMin=0,padding=0),\n", " y=alt.Y('proportion:Q', title='Percent of Samples').axis(format = '%'),\n", " tooltip = [\n", - " alt.Tooltip('interval', title = \"Reads Per Molecule\"),\n", + " alt.Tooltip('interval', title = \"Fragments Per Molecule\"),\n", " alt.Tooltip('proportion', title = \"Percent of Samples\").format('.1%')\n", " ]\n", " )\n", " .properties(\n", - " title=alt.Title('Reads Per Molecule', subtitle = 'Values derived using non-singleton molecules'),\n", - " usermeta={'embedOptions': {'downloadFileName': f'alignments.readspermol'}}\n", + " title=alt.Title('Fragments Per Molecule', subtitle = 'Values derived using non-singleton molecules'),\n", + " usermeta={'embedOptions': {'downloadFileName': 'alignments.fragspermol'}}\n", " )\n", " ).display()\n" ] @@ -480,7 +483,7 @@ }, "outputs": [], "source": [ - "_ycols = [\"Valid\", \"Unique Molecules\", \"Linked\", \"Barcodes\", \"Molecules/Barcode\", \"Reads/Molecule\"]\n", + "_ycols = [\"Valid\", \"Unique Molecules\", \"Linked\", \"Barcodes\", \"Molecules/Barcode\", \"Fragments/Molecule\"]\n", "dropdowny = alt.binding_select(\n", " options= _ycols,\n", " labels=[f\" {i} \" for i in _ycols],\n", @@ -490,18 +493,18 @@ "\n", "with SafeRender(aggregate_df):\n", " (\n", - " alt.Chart(aggregate_df[['Sample', 'Valid', 'Total Reads', \"Unique Molecules\", \"Linked\", \"Barcodes\", \"Molecules/Barcode\", \"Reads/Molecule\"]])\n", + " alt.Chart(aggregate_df[['Sample', 'Valid', 'Fragments', \"Unique Molecules\", \"Linked\", \"Barcodes\", \"Molecules/Barcode\", \"Fragments/Molecule\"]])\n", " .mark_point(stroke = \"white\", strokeWidth = 1)\n", " .add_params(ycol_param)\n", " .transform_calculate(Selected=f'datum[{ycol_param.name}]')\n", " .encode(\n", " y=alt.Y('Selected:Q', title = ''),\n", - " x='Total Reads:Q',\n", + " x='Fragments:Q',\n", " color=alt.Color('Selected:Q').scale(scheme='turbo', reverse = True),\n", - " tooltip = ['Sample', 'Selected:Q', 'Total Reads:Q']\n", + " tooltip = ['Sample', 'Selected:Q', 'Fragments:Q']\n", " )\n", " .properties(\n", - " title=alt.Title('Proportion against Total Reads'),\n", + " title=alt.Title('Proportion against number of fragments'),\n", " usermeta={\n", " 'embedOptions': {'downloadFileName': f'alignments.percent.linked'},\n", " 'padding': {'left': 80, 'right': 20, 'top': 10, 'bottom': 10}\n", @@ -524,6 +527,13 @@ "would be the first number (starting from the biggest) that sums up to at least 5 (50% of 10), which is `3`, because `4` + `3` = 7. The `N90` would be the first number that sums up to at least 9 (90% of 10), which is `2` because `4` + `3` + `2` = 9.\n", ":::\n", "\n", + ":::{dropdown} Fragments vs Reads\n", + "Rather than using 'reads', we report 'fragments', defined as either one paired-end read pair or one single-end\n", + "read (that doesn't have a mate). Alignment records that were designated as duplicates, secondary, supplementary,\n", + "or multimapped across different contigs are omitted in the fragment count. Taken together, this means there will\n", + "**always** be considerably fewer fragments than there are reads and this is typically not a cause for concern.\n", + ":::\n", + "\n", ":::{dropdown} Understanding inferred molecule lengths\n", "Harpy uses the highest and lowest mapping positions of a read cluster sharing the same barcode (incorporating distance-based deconvolution thresholds, if configured to do so) to infer the original molecule length. Given that it's mostly impossible to understand how much of a source molecule the mapped reads actually covered, it would be more appropriate to think of the inferred molecule lengths (and NX measures) to be more like _\"at least this long\"_, since the reads likely did not originate from the ends of the source DNA molecule. The absolute length of the source DNA molecule would be a more useful diagnostic in troubleshooting/modifying the actual linked-read chemistry, whereas the the inferred molecule length describes what length of it is actually represented in the sequences and thus more useful in downstream/analytical contexts.\n", ":::" diff --git a/harpy/notebooks/align_stats.ipynb b/harpy/notebooks/align_stats.ipynb index 137f8c9a9..7d2c83bb8 100644 --- a/harpy/notebooks/align_stats.ipynb +++ b/harpy/notebooks/align_stats.ipynb @@ -603,7 +603,7 @@ " density_plot(binned_df, 'Inferred Molecule Length', samplename, 'inferredlen')\n", " .transform_calculate(bin_kb='datum.bin / 1000')\n", " .encode(\n", - " x=alt.X('bin_kb:Q', title=\"Kilobases (kbp)\").scale(type='log', padding=0).axis(labelAngle=-40),\n", + " x=alt.X('bin_kb:Q', title=\"Kilobases (kbp)\").axis(labelAngle=-40),\n", " y=alt.Y('value:Q', title='Cumulative % Molecules')\n", " .axis(format='%')\n", " .stack(False)\n", diff --git a/harpy/notebooks/hapcut.ipynb b/harpy/notebooks/hapcut.ipynb index 6a9a864b3..30f3d0ddc 100644 --- a/harpy/notebooks/hapcut.ipynb +++ b/harpy/notebooks/hapcut.ipynb @@ -94,7 +94,7 @@ " .add(df['n_snp'].sum(), 'Total SNPs')\n", " .add(round(df['n_snp'].mean(), 0), 'Mean SNPs')\n", " .add(df['n_snp'].median(), 'Median SNPs')\n", - " .add(df['block_length'].max(), 'Longest Haplotype')\n", + " .add(df['block_length'].max(), 'Longest Haplotype', units = \"bp\")\n", " .add(nxx(df['block_length'], 50) / 1000, 'N50', units = \"kb\")\n", " .add(nxx(df['block_length'], 75) / 1000, 'N75', units = \"kb\")\n", " .add(nxx(df['block_length'], 90) / 1000, 'N90', units = \"kb\")\n", @@ -170,18 +170,22 @@ "\n", "with SafeRender(_hist):\n", " _chart = (\n", - " alt.Chart(_hist) \n", + " alt.Chart(_hist)\n", " .transform_calculate(\n", - " bin_kb = 'datum.bin / 1000',\n", - " interval_kbp ='datum.interval + \" bp\"'\n", + " bin_kb='datum.bin / 1000',\n", + " interval_kbp='datum.interval + \" bp\"'\n", " )\n", - " .mark_area()\n", + " .transform_window(\n", + " cumulative_proportion='sum(proportion)',\n", + " sort=[alt.SortField('bin_kb', order='ascending')]\n", + " )\n", + " .mark_line(interpolate='step-after')\n", " .encode(\n", - " x=alt.X('bin_kb:Q', title=\"Haplotype Length (kbp)\").axis(tickMinStep = 5,labelAngle=-40),\n", - " y=alt.Y('proportion:Q', title='Percent of Haplotypes').scale(domain = [0,1]).axis(format='%'),\n", - " tooltip = [\n", + " x=alt.X('bin_kb:Q', title=\"Haplotype Length (kbp)\").axis(tickMinStep=5, labelAngle=-40),\n", + " y=alt.Y('cumulative_proportion:Q', title='Cumulative % of Haplotypes').scale(domain=[0, 1]).axis(format='%'),\n", + " tooltip=[\n", " alt.Tooltip(\"interval_kbp:N\", title=\"Length\"),\n", - " alt.Tooltip(\"proportion:Q\", format = '.2%', title = \"% Haplotypes\")\n", + " alt.Tooltip(\"cumulative_proportion:Q\", format='.2%', title=\"Cumulative % Haplotypes\")\n", " ]\n", " )\n", " ).interactive()\n", @@ -214,7 +218,7 @@ "source": [ "### NX Information\n", "\n", - "An **NX** metric (e.g. **N50**) is the length of the shortest molecule in the group of longest molecules that together\n", + "An **NX** metric (e.g. **N50**) is the length (in **kilobases**) of the shortest molecule in the group of longest molecules that together\n", "represent at least **X%** of the total molecules by length. For example, `N50` would be the shortest molecule in the \n", "group of longest molecules that together represent **50%** of the total molecules by length (sort of like a cumulative median)." ] @@ -364,19 +368,24 @@ "with SafeRender(_hist):\n", " (\n", " alt.Chart(_hist)\n", - " .mark_area()\n", + " .mark_line()\n", " .transform_calculate(y=f'datum[{ycol_param.name}]')\n", + " .transform_window(\n", + " cum_y='sum(y)',\n", + " sort=[alt.SortField('bin')],\n", + " frame=[None, 0]\n", + " )\n", " .encode(\n", - " x=alt.X('bin:Q', bin = 'binned', title = \"Haplotype Length (bp)\").axis(tickMinStep = 5000,labelAngle=-40),\n", - " y=alt.Y('y:Q').title('% of Haplotypes'),\n", - " tooltip = [\n", - " alt.Tooltip('interval:N', title = \"Haplotype Length (bp)\"),\n", - " alt.Tooltip('y:Q', title = '% Haplotypes', format = '.2%')\n", + " x=alt.X('bin:Q', bin='binned', title=\"Haplotype Length (bp)\").axis(tickMinStep=5000, labelAngle=-40),\n", + " y=alt.Y('cum_y:Q').title('Cumulative % of Haplotypes').scale(domain=[0, 1]),\n", + " tooltip=[\n", + " alt.Tooltip('interval:N', title=\"Haplotype Length (bp)\"),\n", + " alt.Tooltip('cum_y:Q', title='Cumulative % Haplotypes', format='.2%')\n", " ]\n", " )\n", " .add_params(ycol_param)\n", " .properties(\n", - " title= \"Haplotype Lengths\",\n", + " title=\"Haplotype Lengths\",\n", " usermeta={'embedOptions': {'downloadFileName': 'haplotypes.contig'}}\n", " )\n", " ).display()\n" diff --git a/harpy/notebooks/samtools_stats.ipynb b/harpy/notebooks/samtools_stats.ipynb index 296e959f5..873b19d84 100644 --- a/harpy/notebooks/samtools_stats.ipynb +++ b/harpy/notebooks/samtools_stats.ipynb @@ -14,7 +14,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "metadata": { "tags": [ "remove-cell" @@ -29,8 +29,7 @@ "from harpy.report.utilities import StopExecution\n", "import os\n", "import polars as pl\n", - "from pathlib import Path\n", - "_ = pl.Config.set_tbl_hide_column_data_types()\n" + "from pathlib import Path\n" ] }, { @@ -318,7 +317,27 @@ "if len(dat_filt) == 0:\n", " print_html(\"No statsfiles were found, or were unable to be properly parsed. Skipping report.\")\n", " raise StopExecution\n", - "dat_filt = dat_filt.select(pl.col(\"*\").exclude(\"% Duplicates (PCR)\", \"% Duplicates (Opt)\"))\n" + "dat_filt = dat_filt.select(pl.col(\"*\").exclude(\"% Duplicates (PCR)\", \"% Duplicates (Opt)\"))\n", + "\n", + "depths = []\n", + "for samplename in dat_filt['Sample']:\n", + " df = pl.read_csv(\n", + " os.path.join(indir, \"coverage\", f\"{samplename}.regions.bed.gz\"),\n", + " separator='\\t',\n", + " has_header=False,\n", + " schema_overrides = [pl.String],\n", + " new_columns=[\"Contig\", \"Position\", \"Position End\", \"Read Depth\"]\n", + " )\n", + " depth = df.select(\n", + " (\n", + " (pl.col(\"Read Depth\") * (pl.col(\"Position End\") - pl.col(\"Position\"))).sum()\n", + " / (pl.col(\"Position End\") - pl.col(\"Position\")).sum()\n", + " ).alias(\"depth\")\n", + " ).item()\n", + " depths.append(round(depth, 2))\n", + " #depths.append(round(df['Read Depth'].mean(), 2))\n", + "\n", + "dat_filt = dat_filt.insert_column(2, pl.Series(\"Depth\", depths))\n" ] }, { @@ -337,6 +356,9 @@ "_mapped = dat_filt[\"% Mapped\"].replace(0, None).mean() or 0\n", "_mappedstd = dat_filt[\"% Mapped\"].replace(0, None).std() or 0\n", "\n", + "_depth = round(dat_filt[\"Depth\"].replace(0, None).mean() or 0, 2)\n", + "_depthstd = round(dat_filt[\"Depth\"].replace(0, None).std() or 0, 2)\n", + "\n", "_isize = round(dat_filt['Average Insert (bp)'].replace(0, None).mean() or 0)\n", "_isizestd = round(dat_filt['Average Insert (bp)'].replace(0, None).std() or 0)\n", "\n", @@ -348,6 +370,7 @@ " .add(len(dat_filt), \"Samples\")\n", " .conditional(_pp, \"Properly Paired\", 80, plus_minus=_ppstd, add_percent=True, digits = 1)\n", " .conditional(_mapped, \"Mapped\", 70, plus_minus=_mappedstd, add_percent=True, digits = 1)\n", + " .add(_depth, \"Coverage Depth\", units = \"X\", plus_minus=_depthstd)\n", " .conditional(_err, \"Error Rate\", 10, plus_minus=_errstd, lower_bad = False, add_percent=True, digits = 3)\n", " .add(_isize, \"Insert Size\", plus_minus=_isizestd, units = \"bp\")\n", ").render()" diff --git a/harpy/notebooks/validate_bam.ipynb b/harpy/notebooks/validate_bam.ipynb index 7fbad410a..5475d9e2c 100644 --- a/harpy/notebooks/validate_bam.ipynb +++ b/harpy/notebooks/validate_bam.ipynb @@ -69,9 +69,11 @@ " .height\n", ")\n", "\n", - "noMItag = (data['noMI'] == data['records']).sum()\n", "noBXtag = (data['noBX'] == data['records']).sum()\n", "noVXtag = (data['noVX'] == data['records']).sum()\n", + "bx_records = data['records'] - data['noBX']\n", + "noMItag = ((bx_records > 0) & (data['noMI'] == bx_records)).sum()\n", + "noVXtag = ((bx_records > 0) & (data['noVX'] == bx_records)).sum()\n", "bxnotlast = (data['bxNotLast'] > 0).sum()\n", "format_issues = (data['badBX'] > 0).sum()\n" ] @@ -106,7 +108,7 @@ "## Metrics\n", "The `harpy validate bam` command created a `validate.bam.tsv` file\n", "that summarizes the results included in this report. This file contains\n", - "a tab-delimited table with the columns: `file`,\t`nameMismatch`, `alignments`, `format`, `noBX`, `noVX`, `noMI`, and `badBX`.\n", + "a tab-delimited table with the columns: `file`, `records`, `nameMismatch`, `noMI`, `noBX`, `noVX`, `bxNotLast`, and `badBX`.\n", "These columns are defined in **Supporting info** below.\n", "\n", "Look over the metric descriptions to better understand what these validation assessments mean and their severities. At the bottom of this report\n", @@ -114,12 +116,12 @@ "\n", ":::{note} Validation logic\n", ":open:\n", - "Severity: 🔶 moderate | 🛑 serious\n", + "Severity: ⬜ mild | 🔶 moderate | 🛑 serious\n", "| column | severity | pass condition | fail condition |\n", "| :--------------- | :------- | :-------------------------------------------------------------------------------------- | :------------------------------------------------------------------------------ |\n", "| **nameMismatch** | 🛑 | the file name matches the `@RG ID:` tag in the header | file name does not match `@RG ID:` in the header |\n", "| **noBX** | 🛑 | **any** `BX:Z` tags present | **all** alignments lack `BX:Z` tag |\n", - "| **noMI** | 🔶 | **any** alignments with `BX:Z` tag also have `MI` | **all** reads that have `BX:Z` tag present don't have `MI` tag |\n", + "| **noMI** | ⬜/🔶 | **any** alignments with `BX:Z` tag also have `MI` | **all** reads that have `BX:Z` tag present don't have `MI` tag |\n", "| **noVX** | 🔶/🛑 | **any** alignments with `BX:Z` tag also have `VX:i` | **all** reads that have `BX:Z` tag present don't have `VX:i` tag |\n", "| **bxNotLast** | 🔶 | **all** reads have `BX:Z` as final tag in alignment records | **at least 1 read** doesn't have `BX:Z` tag as final tag |\n", "| **badBX** | 🛑 | **all** alignments with `BX:Z` tag have properly formatted barcodes for their chemistry | **any** `BX:Z` barcodes have incorrectly formatted barcodes for their chemistry |\n", @@ -167,7 +169,7 @@ "If you expect all or some of your reads should have `BX:Z` tags, then further investigation is necessary\n", "\n", "**noMI**: Alignment records lack `MI` tag (`MI:i` or `MI:Z`)\n", - "- `MI` tags are optional, they represent a confident molecule assignment, usually following some kind of deconvolution. Some software recognizes it, but most linked-read software aims for the `BX:Z` tag\n", + "- `MI` tags are optional, they represent a confident molecule assignment, usually following some kind of deconvolution. Some software recognizes it, but most linked-read software aims for the `BX:Z` tag. Presence of an `MI` tag is only checked if the record already has a `BX` tag.\n", "\n", "**noVX**: Alignment records that have a `BX:Z` tag lack a `VX:i` validation tag compliant with the standard linked-read data format\n", "- `VX:i` tags are not mandatory _per se_, however we strongly encourage and promote their use to unify the differences between linked-read barcode styles\n", diff --git a/harpy/notebooks/validate_fastq.ipynb b/harpy/notebooks/validate_fastq.ipynb index d6d2f7a74..ce120dd56 100644 --- a/harpy/notebooks/validate_fastq.ipynb +++ b/harpy/notebooks/validate_fastq.ipynb @@ -65,7 +65,7 @@ "attention = (\n", " data\n", " .with_columns((pl.col('reads') == pl.col('noBX')).alias('allmissing'))\n", - " .select(pl.sum_horizontal(['badBX', 'badSamSpec', 'allmissing']).alias('attention_sum'))\n", + " .select(pl.sum_horizontal(['badBX', 'badSamSpec', 'allmissing', 'noVX']).alias('attention_sum'))\n", " .filter(pl.col('attention_sum') > 0)\n", " .height\n", ")\n" @@ -96,7 +96,7 @@ "## Metrics\n", "The `harpy validate fastq` command created a `validate.fastq.tsv` file in the specified\n", "output directory that summarizes the results that are included in this report. This file\n", - "contains a tab-delimited table with the columns: `file`, `reads`, `noBX`, `badBX`, and `badSamSpec`.\n", + "contains a tab-delimited table with the columns: `file`, `reads`, `noBX`, `noVX`, `bxNotLast`, `badBX`, and `badSamSpec`.\n", "\n", ":::{note} Validation logic\n", ":open: false\n", @@ -107,11 +107,11 @@ "| **badSamSpec** | 🛑 | **all** reads have proper `TAG:TYPE:VALUE` comments | **any** reads have incorrectly formatted comments |\n", "| **bxNotLast** | ⬜/🔶 | **all** reads have `BX:Z` as final comment [^2] | **at least 1 read** with a `BX:Z` tag doesn’t have it as the terminal tag [^2] |\n", "| **noBX** | ⬜/🔶 | **any** barcodes present | **all** reads lack barcodes |\n", - "| **noVX** | ⬜/🔶 | reads with `BX:Z` tag have also have `VX:Z` tag | reads with `BX:Z` tag don't have a `VX:Z` tag |\n", + "| **noVX** | ⬜/🔶 | reads with `BX:Z` tag have also have `VX:i` tag | reads with `BX:Z` tag don't have a `VX:i` tag |\n", "\n", ":::\n", "[^1]:\n", - " Standard format refers to a barcode in a `BX:Z` SAM tag and its validation in a `VX:i` SAM tag. See [the LASTQ standard](https://pdimens.github.io/lastq/standard/).\n", + " Standard format refers to a barcode in a `BX:Z` SAM tag and its validation in a `VX:i` SAM tag. See [the LASTQ standard](https://blinkseq.github.io/lastq/standard/).\n", "\n", "[^2]:\n", " Things regarding `BX:Z` tags are ignored in non-haplotagging chemistries" diff --git a/harpy/report/static.py b/harpy/report/static.py index 9637d0e76..1b629290e 100644 --- a/harpy/report/static.py +++ b/harpy/report/static.py @@ -5,6 +5,8 @@ import tempfile from pathlib import Path +from rich.console import Group + from harpy.common.printing import HarpyPrint try: @@ -15,15 +17,52 @@ FRONTMATTER_RE = re.compile(r"\A---\s*\n(.*?\n)---\s*\n?", re.DOTALL) +def has_nbconvert(): + hp = HarpyPrint() + try: + has = subprocess.run( + ["jupyter", "nbconvert", "--version"], + capture_output=True + ).returncode == 0 + except FileNotFoundError: + has = False + if not has: + _table = hp.table() + _table.add_column("tool") + _table.add_column("installation command", style = "green") + _table.add_row('pip', 'pip install -U nbconvert') + _table.add_row('conda', 'conda install -c conda-forge nbconvert') + _table.add_row('pixi', 'pixi add nbconvert') + hp.error( + "Missing dependency", + "jupyter nbconvert is not found on the PATH and is required to proceed.", + Group("It can be installed using one of these methods:", _table) + ) + +def has_monolith(): + if not shutil.which('monolith'): + hp = HarpyPrint() + _table = hp.table() + _table.add_column("tool") + _table.add_column("installation", style = "green") + _table.add_row('cargo', 'cargo install monolith') + _table.add_row('prebuilt binary', 'add binary to your PATH from https://github.com/Y2Z/monolith/releases') + hp.error( + "Missing dependency", + "Monolith is required to flatten an HTML notebook but was not found on the PATH.", + Group("Harpy does not provide it, but it can be installed using:", _table) + ) class ReportStatic(): - def __init__(self, notebook: str, quiet: bool, static: bool): + def __init__(self, quiet: bool, static: bool): self.quiet: bool = quiet self.static: bool = static - self.notebook: str = notebook self.hp = HarpyPrint() self.hp.console.soft_wrap = True self.nbc_log = "ERROR" if quiet else 30 + has_nbconvert() + if static: + has_monolith() def render_frontmatter_cell(self, nb: dict) -> None: """ @@ -80,8 +119,8 @@ def run(self, cmd: list[str], **kwargs) -> None: subprocess.run(cmd, check=True, **kwargs) - def convert(self): - nb_path: Path = Path(self.notebook).resolve() + def convert(self, notebook: str): + nb_path: Path = Path(notebook).resolve() nb_name = nb_path.stem out_path: Path = nb_path.with_name(f"{nb_name}.html") diff --git a/harpy/report/utilities.py b/harpy/report/utilities.py index 869e562d6..ddde3af7d 100644 --- a/harpy/report/utilities.py +++ b/harpy/report/utilities.py @@ -259,14 +259,4 @@ def __init__(self, report_dir): self.gc_curves_after = pl.DataFrame(gc_curves["after"]) self.gc_curves_r2 = pl.DataFrame(gc_curves_r2["before"]) self.gc_curves_r2_after = pl.DataFrame(gc_curves_r2["after"]) - - -def check_tool(name: str, hint: str) -> None: - if shutil.which(name) is None: - hp = HarpyPrint() - hp.error( - "Missing dependency", - f"{name} is not found on the PATH and required to proceed.", - hint - ) \ No newline at end of file diff --git a/harpy/snakefiles/align.smk b/harpy/snakefiles/align.smk new file mode 100644 index 000000000..3c714afd8 --- /dev/null +++ b/harpy/snakefiles/align.smk @@ -0,0 +1,242 @@ +import os + +localrules: all +wildcard_constraints: + sample = r"[a-zA-Z0-9._-]+" + +WORKFLOW = config.get('Workflow') or {} +PARAMETERS = config.get('Parameters') or {} +REPORTS = WORKFLOW.get("reports") or {} +INPUTS = config['Inputs'] +VERSION = WORKFLOW.get('harpy-version', 'latest') + + +lr_type = WORKFLOW.get("linkedreads", {}).get("type", 'none') +bx_tag = WORKFLOW.get("linkedreads", {}).get("standardized", {}).get("BX", False) +vx_tag = WORKFLOW.get("linkedreads", {}).get("standardized", {}).get("VX", False) +skip_reports = REPORTS.get("skip", False) +molecule_distance = PARAMETERS.get("distance-threshold", 0) +keep_unmapped = PARAMETERS.get("keep-unmapped", False) +extra = PARAMETERS.get("extra", "") +windowsize = PARAMETERS.get("depth-windowsize", 50000) +genomefile = INPUTS["reference"] + +ignore_bx = lr_type == "none" +bn = os.path.basename(genomefile) +workflow_geno = f"workflow/reference/{bn}" + +aligner = WORKFLOW.get("name", "align_bwa").split("_")[-1] +include: f"align_{aligner}.smk" + +rule sort: + retries: 3 + input: + ref = workflow_geno, + bam = f"{aligner}/{{sample}}.{aligner}.bam" + output: + bam = temp("sort/{sample}.sort.bam"), + stats = "reports/data/samtools_stats/{sample}.raw.stats", + tmp = temp(directory("sort/{sample}_tmp")) + log: + "logs/sort/{sample}.sort.log" + params: + sortthreads = lambda wc, threads: threads - 1 + threads: + 4 + resources: + tmpdir = lambda wc: f"sort/{wc.sample}_tmp", + mem_mb_per_thread = lambda wc, attempt: 3000 // attempt + shell: + """ + mkdir -p {resources.tmpdir} + {{ + samtools fixmate -z on -m -u {input.bam} - | + samtools sort -@ {params.sortthreads} -M -T {resources.tmpdir} -o {output.bam} -u -l 0 -m {resources.mem_mb_per_thread}M - + samtools stats -@ {params.sortthreads} -d -x -r {input.ref} {output.bam} > {output.stats} + }} 2> {log} + """ + +rule mark_duplicates: + priority: 1 + input: + fq = get_fq, + bam = "sort/{sample}.sort.bam" + output: + bam = "{sample}.bam" if lr_type == "none" or (bx_tag and vx_tag) else temp("markdup/{sample}.bam"), + stats = "reports/data/markdup/{sample}.markdup", + tmp = temp(directory("markdup/{sample}_tmp")) + log: + "logs/markdup/{sample}.markdup.log" + params: + bx_mode = "-S --barcode-tag BX" if not ignore_bx else "-S", + quality = PARAMETERS.get('min-map-quality', 30), + unmapped = "-F 4" if not keep_unmapped else "", + mdthreads = lambda wc, threads: threads - 1 + resources: + tmpdir = lambda wc: f"markdup/{wc.sample}_tmp" + threads: + 4 + shell: + """ + mkdir -p {resources.tmpdir} + OPT=$(harpy-utils optical-dist-fq {input.fq}) + {{ + samtools view -h -u -q {params.quality} {params.unmapped} {input.bam} | + samtools markdup -@ {params.mdthreads} -T {resources.tmpdir} {params.bx_mode} -d $OPT -f {output.stats} - {output.bam} + }} 2> {log} + """ + +if lr_type != "none" and not (bx_tag and vx_tag): + rule standardize: + input: + "markdup/{sample}.bam" + output: + "{sample}.bam" + log: + "logs/{sample}.std.log" + threads: + 2 + shell: + "djinn-standardize --threads {threads} {input} > {output} 2> {log}" + +rule depth_stats: + input: + "{sample}.bam.bai", + bam = "{sample}.bam" + output: + "reports/data/coverage/{sample}.regions.bed.gz" + params: + f"-b {windowsize}", + "-n --fast-mode" + log: + "logs/depthstats/{sample}.mosdepth.log" + threads: + 2 + conda: + "envs/qc.yaml" + container: + f"docker://pdimens/harpy:qc_{VERSION}" + shell: + """ + mosdepth {params} -t 1 reports/data/coverage/{wildcards.sample} {input.bam} 2> {log} + rm -f reports/data/coverage/{wildcards.sample}.mosdepth* reports/data/coverage/{wildcards.sample}*.csi + """ + +rule sample_stats: + input: + "{sample}.bam" + output: + temp("{sample}.bam.bai"), + stats = "reports/data/samtools_stats/{sample}.filtered.stats" + log: + "logs/stats/{sample}.stats.log" + threads: + 2 + shell: + """ + {{ + samtools index {input} + samtools stats -@ 1 -x -d {input} > {output.stats} + }} 2> {log} + """ + +rule molecule_coverage: + input: + fai = f"{workflow_geno}.fai", + stats = "reports/data/lrstats/{sample}.lrstats.gz" + output: + "reports/data/coverage/{sample}.molcov.gz" + log: + "logs/stats/{sample}.molstats.log" + params: + windowsize + shell: + "harpy-utils molecule-coverage -w {params} {input} 2> {log} | gzip > {output}" + +rule molecule_stats: + input: + "{sample}.bam" + output: + "reports/data/lrstats/{sample}.lrstats.gz" + log: + "logs/molcov/{sample}.molcov.log" + params: + molecule_distance + shell: + "harpy-utils bx-stats-sam -d {params} {input} 2> {log} | gzip > {output}" + +rule alignment_report: + input: + collect("reports/data/markdup/{sample}.markdup", sample = samplenames), + collect("reports/data/samtools_stats/{sample}.{data}.stats", sample = samplenames, data = ["raw", "filtered"]), + collect("reports/data/coverage/{sample}.regions.bed.gz", sample = samplenames), + ipynb = f"workflow/samtools_stats.ipynb" + output: + tmp = temp(f"reports/{aligner}.summary.tmp.ipynb"), + ipynb = f"reports/{aligner}.summary.ipynb" + params: + lr_type = lr_type, + indir = "-p indir " + os.path.abspath("reports/data") + log: + f"logs/reports/{aligner}.report.log" + shell: + """ + export IPYTHONDIR=/tmp/ipython-align-stats + {{ + papermill -k ipython-harpy --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params.indir} + harpy-utils process-notebook {output.tmp} {params.lr_type} > {output.ipynb} + }} 2> {log} + """ + +rule sample_reports: + input: + lrstats = "reports/data/lrstats/{sample}.lrstats.gz", + coverage = "reports/data/coverage/{sample}.regions.bed.gz", + molcov = "reports/data/coverage/{sample}.molcov.gz", + ipynb = f"workflow/align_stats.ipynb" + output: + tmp = temp("reports/{sample}.tmp.ipynb"), + ipynb = "reports/{sample}.ipynb" + params: + placeholders = f'{aligner} {lr_type}', + papermill = f'-p platform {lr_type} -p basedir {os.path.abspath("reports/data")} -p mol_dist {molecule_distance} -p windowsize {windowsize}', + samplename = lambda wc: "-p samplename " + wc.get("sample") + log: + "logs/reports/{sample}.report.log" + shell: + """ + export IPYTHONDIR=/tmp/ipython-{wildcards.sample}.rpt + {{ + papermill -k ipython-harpy --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params.papermill} {params.samplename} + harpy-utils process-notebook {output.tmp} {wildcards.sample} {params.placeholders} > {output.ipynb} + }} 2> {log} + """ + +rule linked_read_report: + input: + collect("reports/data/lrstats/{sample}.lrstats.gz", sample = samplenames), + ipynb = f"workflow/align_lrstats.ipynb" + output: + tmp = temp("reports/linkedreads.summary.tmp.ipynb"), + ipynb = "reports/linkedreads.summary.ipynb" + params: + lr_type = lr_type, + indir = "-p indir " + os.path.abspath("reports/data/lrstats") + log: + f"logs/reports/lrstats.report.log" + shell: + """ + {{ + export IPYTHONDIR=/tmp/ipython-lr-stats + papermill -k ipython-harpy --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params.indir} + harpy-utils process-notebook {output.tmp} {params.lr_type} > {output.ipynb} + }} 2> {log} + """ + +rule all: + default_target: True + input: + bams = collect("{sample}.bam", sample = samplenames), + reports = collect("reports/{sample}.ipynb", sample = samplenames) if not skip_reports and not ignore_bx else [], + align_report = f"reports/{aligner}.summary.ipynb" if (not skip_reports and len(samplenames) > 1) else [], + bx_report = "reports/linkedreads.summary.ipynb" if (not skip_reports and not ignore_bx and len(samplenames) > 1) else [] diff --git a/harpy/snakefiles/align_arachne.smk b/harpy/snakefiles/align_arachne.smk new file mode 100644 index 000000000..f09a43720 --- /dev/null +++ b/harpy/snakefiles/align_arachne.smk @@ -0,0 +1,363 @@ +import os +import re + +localrules: all +wildcard_constraints: + sample = r"[a-zA-Z0-9._-]+" + +WORKFLOW = config.get('Workflow') or {} +PARAMETERS = config.get('Parameters') or {} +REPORTS = WORKFLOW.get("reports") or {} +INPUTS = config['Inputs'] +VERSION = WORKFLOW.get('harpy-version', 'latest') + +skip_reports = REPORTS.get("skip", False) +molecule_distance = PARAMETERS.get("distance-threshold", 0) +keep_unmapped = PARAMETERS.get("keep-unmapped", False) +extra = PARAMETERS.get("extra", "") +windowsize = PARAMETERS.get("depth-windowsize", 50000) +fqlist = INPUTS["fastq"] +genomefile = INPUTS["reference"] +centromeres = INPUTS.get("centromeres", None) +lr_type = WORKFLOW.get("linkedreads", {}).get("type", 'none') + +bn = os.path.basename(genomefile) +workflow_geno = f"workflow/reference/{bn}" +genome_zip = True if bn.lower().endswith(".gz") else False +geno_idx = f"{workflow_geno}.gzi" if genome_zip else f"{workflow_geno}.fai" +bn_r = r"([_\.][12]|[_\.][FR]|[_\.]R[12](?:\_00[0-9])*)?\.((fastq|fq)(\.gz)?)$" +samplenames = {re.sub(bn_r, "", os.path.basename(i), flags = re.IGNORECASE) for i in fqlist} +d = dict(zip(samplenames, samplenames)) + +def get_fq(wildcards): + # returns a list of fastq files for read 1 based on *wildcards.sample* e.g. + r = re.compile(fr".*/({re.escape(wildcards.sample)}){bn_r}", flags = re.IGNORECASE) + return sorted(list(filter(r.match, fqlist))[:2]) + +rule process_reference: + input: + genomefile + output: + geno = workflow_geno, + #TODO THREADS SUPPORT? + bwa_idx = multiext(workflow_geno, '.amb', '.ann', '.bwt', '.pac', '.sa', '.l2b', '.mbw'), + fai = f"{workflow_geno}.fai", + gzi = f"{workflow_geno}.gzi" if genome_zip else [] + log: + f"{workflow_geno}.preprocess.log" + params: + genome_zip + threads: + workflow.cores + conda: + "envs/align.yaml" + container: + f"docker://pdimens/harpy:align_{VERSION}" + shell: + """ + {{ + if (file {input} | grep -q compressed ) ;then + # is regular gzipped, needs to be BGzipped + seqtk seq {input} | bgzip -c > {output.geno} + else + cp -f {input} {output.geno} + fi + + if [ "{params}" = "True" ]; then + samtools faidx --gzi-idx {output.gzi} --fai-idx {output.fai} {output.geno} + else + samtools faidx --fai-idx {output.fai} {output.geno} + fi + + minibwa index -t {threads} {output.geno} + arachne index -t {threads} {output.geno} + }} 2> {log} + """ + +rule process_fastq: + input: + get_fq + output: + collect("arachne/prep/{{sample}}.{prefix}.{FR}.fq.gz", prefix = ['arachne', 'invalid'], FR = ['R1', 'R2']) + threads: + 4 + conda: + "envs/align.yaml" + container: + f"docker://pdimens/harpy:align_{VERSION}" + shell: + """ + mkdir -p arachne/prep + arachne prep -t {threads} arachne/prep/{wildcards.sample} {input} + """ + +rule arachne_align: + input: + multiext(workflow_geno, '.amb', '.ann', '.bwt', '.pac', '.sa'), + ref = workflow_geno, + R1 = "arachne/prep/{sample}.arachne.R1.fq.gz", + R2 = "arachne/prep/{sample}.arachne.R2.fq.gz", + centromeres = centromeres if centromeres else [] + output: + temp("arachne/align/{sample}.arachne.bam") + log: + "logs/arachne/{sample}.arachne.log" + params: + RG_tag = lambda wc: "-s " + wc.get("sample"), + dist = f"-d {molecule_distance}", + extra = extra + threads: + 12 + conda: + "envs/align.yaml" + container: + f"docker://pdimens/harpy:align_{VERSION}" + shell: + """ + mkdir -p {resources.tmpdir} + arachne align -t {threads} {params} {input.ref} {input.R1} {input.R2} | + samtools view -l 0 -O BAM - | samtools sort -u > {output} 2> {log} + """ + +rule bwa_align: + input: + multiext(workflow_geno, ".l2b", ".mbw"), + ref = workflow_geno, + R1 = "arachne/prep/{sample}.invalid.R1.fq.gz", + R2 = "arachne/prep/{sample}.invalid.R2.fq.gz" + output: + bam = temp("bwa/{sample}.bwa.bam"), + tmp = temp(directory("bwa/{sample}_tmp")) + log: + "logs/bwa/{sample}.bwa.log" + params: + RG_tag = lambda wc: "-R \"@RG\\tID:" + wc.get("sample") + "\\tSM:" + wc.get("sample") + "\"", + extra = extra + threads: + 12 + resources: + tmpdir = lambda wc: f"bwa/{wc.sample}_tmp" + conda: + "envs/align.yaml" + container: + f"docker://pdimens/harpy:align_{VERSION}" + shell: + """ + mkdir -p {resources.tmpdir} + {{ + minibwa map -t {threads} {params} {input.ref} {input.R1} {input.R2} | + samtools collate -T {resources.tmpdir} -O -u - + }} 2> {log} > {output.bam} + """ + +rule sort: + retries: 3 + input: + ref = workflow_geno, + bam = "bwa/{sample}.bwa.bam" + output: + bam = temp("sort/{sample}.sort.bam"), + stats = "reports/data/samtools_stats/{sample}.raw.stats", + tmp = temp(directory("sort/{sample}_tmp")) + log: + "logs/sort/{sample}.sort.log" + params: + sortthreads = lambda wc, threads: threads - 1 + threads: + 4 + resources: + tmpdir = lambda wc: f"sort/{wc.sample}_tmp", + mem_mb_per_thread = lambda wc, attempt: 3000 // attempt + shell: + """ + mkdir -p {resources.tmpdir} + {{ + samtools fixmate -z on -m -u {input.bam} - | + samtools sort -@ {params.sortthreads} -M -T {resources.tmpdir} -o {output.bam} -u -l 0 -m {resources.mem_mb_per_thread}M - + samtools stats -@ {params.sortthreads} -d -x -r {input.ref} {output.bam} > {output.stats} + }} 2> {log} + """ + +rule mark_duplicates: + priority: 1 + input: + fq = get_fq, + bam = "sort/{sample}.sort.bam" + output: + bam = "markdup/{sample}.bam", + stats = "reports/data/markdup/{sample}.markdup", + tmp = temp(directory("markdup/{sample}_tmp")) + log: + "logs/markdup/{sample}.markdup.log" + params: + bx_mode = "-S", + quality = PARAMETERS.get('min-map-quality', 30), + unmapped = "-F 4" if not keep_unmapped else "", + mdthreads = lambda wc, threads: threads - 1 + resources: + tmpdir = lambda wc: f"markdup/{wc.sample}_tmp" + threads: + 4 + shell: + """ + mkdir -p {resources.tmpdir} + OPT=$(harpy-utils optical-dist-fq {input.fq}) + {{ + samtools view -h -u -q {params.quality} {params.unmapped} {input.bam} | + samtools markdup -@ {params.mdthreads} -T {resources.tmpdir} {params.bx_mode} -d $OPT -f {output.stats} - {output.bam} + }} 2> {log} + """ + +rule combine_alignments: + input: + "arachne/align/{sample}.arachne.bam", + "markdup/{sample}.bam" + output: + "{sample}.bam" + log: + "logs/concat/{sample}.concat.log" + shell: + "samtools cat -o {output} {input} 2> {log}" + +rule depth_stats: + input: + "{sample}.bam.bai", + bam = "{sample}.bam" + output: + "reports/data/coverage/{sample}.regions.bed.gz" + params: + f"-b {windowsize}", + "-n --fast-mode" + log: + "logs/depthstats/{sample}.mosdepth.log" + threads: + 2 + conda: + "envs/qc.yaml" + container: + f"docker://pdimens/harpy:qc_{VERSION}" + shell: + """ + mosdepth {params} -t 1 reports/data/coverage/{wildcards.sample} {input.bam} 2> {log} + rm -f reports/data/coverage/{wildcards.sample}.mosdepth* reports/data/coverage/{wildcards.sample}*.csi + """ + +rule sample_stats: + input: + "{sample}.bam" + output: + temp("{sample}.bam.bai"), + stats = "reports/data/samtools_stats/{sample}.filtered.stats" + log: + "logs/stats/{sample}.stats.log" + threads: + 2 + shell: + """ + {{ + samtools index {input} + samtools stats -@ 1 -x -d {input} > {output.stats} + }} 2> {log} + """ + +rule molecule_coverage: + input: + fai = f"{workflow_geno}.fai", + stats = "reports/data/lrstats/{sample}.lrstats.gz" + output: + "reports/data/coverage/{sample}.molcov.gz" + log: + "logs/stats/{sample}.molstats.log" + params: + windowsize + shell: + "harpy-utils molecule-coverage -w {params} {input} 2> {log} | gzip > {output}" + +rule molecule_stats: + input: + "{sample}.bam" + output: + "reports/data/lrstats/{sample}.lrstats.gz" + log: + "logs/molcov/{sample}.molcov.log" + params: + molecule_distance + shell: + "harpy-utils bx-stats-sam -d {params} {input} 2> {log} | gzip > {output}" + +rule alignment_report: + input: + collect("reports/data/markdup/{sample}.markdup", sample = samplenames), + collect("reports/data/samtools_stats/{sample}.{data}.stats", sample = samplenames, data = ["raw", "filtered"]), + collect("reports/data/coverage/{sample}.regions.bed.gz", sample = samplenames), + ipynb = f"workflow/samtools_stats.ipynb" + output: + tmp = temp("reports/bwa.summary.tmp.ipynb"), + ipynb = "reports/bwa.summary.ipynb" + params: + indir = "-p indir " + os.path.abspath("reports/data") + log: + f"logs/reports/bwa.report.log" + shell: + """ + export IPYTHONDIR=/tmp/ipython-bwa-stats + {{ + papermill -k ipython-harpy --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params.indir} + harpy-utils process-notebook {output.tmp} > {output.ipynb} + }} 2> {log} + """ + +rule sample_reports: + input: + lrstats = "reports/data/lrstats/{sample}.lrstats.gz", + coverage = "reports/data/coverage/{sample}.regions.bed.gz", + molcov = "reports/data/coverage/{sample}.molcov.gz", + ipynb = f"workflow/align_stats.ipynb" + output: + tmp = temp("reports/{sample}.tmp.ipynb"), + ipynb = "reports/{sample}.ipynb" + params: + lr_type = lr_type, + basedir = "-p basedir " + os.path.abspath("reports/data"), + mol_dist = f"-p mol_dist {molecule_distance}", + window_size = f"-p windowsize {windowsize}", + samplename = lambda wc: "-p samplename " + wc.get("sample"), + log: + "logs/reports/{sample}.report.log" + shell: + """ + export IPYTHONDIR=/tmp/ipython-{wildcards.sample}.rpt + {{ + papermill -k ipython-harpy --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} -p platform {params} + harpy-utils process-notebook {output.tmp} {wildcards.sample} minibwa {params.lr_type} > {output.ipynb} + }} 2> {log} + """ + +rule linked_read_report: + input: + collect("reports/data/lrstats/{sample}.lrstats.gz", sample = samplenames), + ipynb = f"workflow/align_lrstats.ipynb" + output: + tmp = temp("reports/linkedreads.summary.tmp.ipynb"), + ipynb = "reports/linkedreads.summary.ipynb" + params: + lr_type = lr_type, + indir = "-p indir " + os.path.abspath("reports/data/lrstats") + log: + f"logs/reports/lrstats.report.log" + shell: + """ + {{ + export IPYTHONDIR=/tmp/ipython-bwa.lr + papermill -k ipython-harpy --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params.indir} + harpy-utils process-notebook {output.tmp} {params.lr_type} > {output.ipynb} + }} 2> {log} + """ + +rule all: + default_target: True + input: + bams = collect("{sample}.bam", sample = samplenames), + reports = collect("reports/{sample}.ipynb", sample = samplenames) if not skip_reports else [], + align_report = "reports/bwa.summary.ipynb" if (not skip_reports and len(samplenames) > 1) else [], + bx_report = "reports/linkedreads.summary.ipynb" if (not skip_reports and len(samplenames) > 1) else [] diff --git a/harpy/snakefiles/align_bwa.smk b/harpy/snakefiles/align_bwa.smk index d8952d7ca..6af48d8fd 100644 --- a/harpy/snakefiles/align_bwa.smk +++ b/harpy/snakefiles/align_bwa.smk @@ -1,7 +1,6 @@ import os import re -localrules: all wildcard_constraints: sample = r"[a-zA-Z0-9._-]+" @@ -11,19 +10,12 @@ REPORTS = WORKFLOW.get("reports") or {} INPUTS = config['Inputs'] VERSION = WORKFLOW.get('harpy-version', 'latest') -lr_type = WORKFLOW.get("linkedreads", {}).get("type", 'none') -bx_tag = WORKFLOW.get("linkedreads", {}).get("standardized", {}).get("BX", False) -vx_tag = WORKFLOW.get("linkedreads", {}).get("standardized", {}).get("VX", False) -skip_reports = REPORTS.get("skip", False) illumina_old = PARAMETERS.get("illumina-format-old", False) -molecule_distance = PARAMETERS.get("distance-threshold", 0) -keep_unmapped = PARAMETERS.get("keep-unmapped", False) extra = PARAMETERS.get("extra", "") windowsize = PARAMETERS.get("depth-windowsize", 50000) fqlist = INPUTS["fastq"] genomefile = INPUTS["reference"] -ignore_bx = lr_type == "none" bn = os.path.basename(genomefile) workflow_geno = f"workflow/reference/{bn}" genome_zip = True if bn.lower().endswith(".gz") else False @@ -42,13 +34,15 @@ rule process_reference: genomefile output: geno = workflow_geno, - bwa_idx = multiext(workflow_geno, ".0123", ".amb", ".ann", ".bwt.2bit.64", ".pac"), + bwa_idx = multiext(workflow_geno, ".l2b", ".mbw"), fai = f"{workflow_geno}.fai", gzi = f"{workflow_geno}.gzi" if genome_zip else [] log: f"{workflow_geno}.preprocess.log" params: genome_zip + threads: + workflow.cores conda: "envs/align.yaml" container: @@ -69,13 +63,13 @@ rule process_reference: samtools faidx --fai-idx {output.fai} {output.geno} fi - bwa-mem2 index {output.geno} + minibwa index -t {threads} {output.geno} }} 2> {log} """ rule align: input: - multiext(workflow_geno, ".0123", ".amb", ".ann", ".bwt.2bit.64", ".pac"), + multiext(workflow_geno, ".l2b", ".mbw"), ref = workflow_geno, fastq = get_fq output: @@ -85,7 +79,7 @@ rule align: "logs/bwa/{sample}.bwa.log" params: RG_tag = lambda wc: "-R \"@RG\\tID:" + wc.get("sample") + "\\tSM:" + wc.get("sample") + "\"", - static = "-m 10 -C -v 2 -T 10" if illumina_old else "-v 2 -T 10 -m 10", + static = "-x sr -y" if illumina_old else "-x sr", extra = extra threads: 12 @@ -99,220 +93,7 @@ rule align: """ mkdir -p {resources.tmpdir} {{ - bwa-mem2 mem -t {threads} {params} {input.ref} {input.fastq} | + minibwa map -t {threads} {params} {input.ref} {input.fastq} | samtools collate -T {resources.tmpdir} -O -u - }} 2> {log} > {output.bam} """ - -rule sort: - retries: 3 - input: - ref = workflow_geno, - bam = "bwa/{sample}.bwa.bam" - output: - bam = temp("sort/{sample}.sort.bam"), - stats = "reports/data/samtools_stats/{sample}.raw.stats", - tmp = temp(directory("sort/{sample}_tmp")) - log: - "logs/sort/{sample}.sort.log" - params: - sortthreads = lambda wc, threads: threads - 1 - threads: - 4 - resources: - tmpdir = lambda wc: f"sort/{wc.sample}_tmp", - mem_mb_per_thread = lambda wc, attempt: 3000 // attempt - shell: - """ - mkdir -p {resources.tmpdir} - {{ - samtools fixmate -z on -m -u {input.bam} - | - samtools sort -@ {params.sortthreads} -M -T {resources.tmpdir} -o {output.bam} -u -l 0 -m {resources.mem_mb_per_thread}M - - samtools stats -@ {params.sortthreads} -d -x -r {input.ref} {output.bam} > {output.stats} - }} 2> {log} - """ - -rule mark_duplicates: - priority: 1 - input: - fq = get_fq, - bam = "sort/{sample}.sort.bam" - output: - bam = "{sample}.bam" if lr_type == "none" or (bx_tag and vx_tag) else temp("markdup/{sample}.bam"), - stats = "reports/data/markdup/{sample}.markdup", - tmp = temp(directory("markdup/{sample}_tmp")) - log: - "logs/markdup/{sample}.markdup.log" - params: - bx_mode = "-S --barcode-tag BX" if not ignore_bx else "-S", - quality = PARAMETERS.get('min-map-quality', 30), - unmapped = "-F 4" if not keep_unmapped else "", - mdthreads = lambda wc, threads: threads - 1 - resources: - tmpdir = lambda wc: f"markdup/{wc.sample}_tmp" - threads: - 4 - shell: - """ - mkdir -p {resources.tmpdir} - OPT=$(harpy-utils optical-dist-fq {input.fq}) - {{ - samtools view -h -u -q {params.quality} {params.unmapped} {input.bam} | - samtools markdup -@ {params.mdthreads} -T {resources.tmpdir} {params.bx_mode} -d $OPT -f {output.stats} - {output.bam} - }} 2> {log} - """ - -if lr_type != "none" and not (bx_tag and vx_tag): - rule standardize: - input: - "markdup/{sample}.bam" - output: - "{sample}.bam" - log: - "logs/{sample}.std.log" - threads: - 2 - shell: - "djinn-standardize --threads {threads} {input} > {output} 2> {log}" - -rule depth_stats: - input: - "{sample}.bam.bai", - bam = "{sample}.bam" - output: - "reports/data/coverage/{sample}.regions.bed.gz" - params: - f"-b {windowsize}", - "-n --fast-mode" - log: - "logs/depthstats/{sample}.mosdepth.log" - threads: - 2 - conda: - "envs/qc.yaml" - container: - f"docker://pdimens/harpy:qc_{VERSION}" - shell: - """ - mosdepth {params} -t 1 reports/data/coverage/{wildcards.sample} {input.bam} 2> {log} - rm -f reports/data/coverage/{wildcards.sample}.mosdepth* reports/data/coverage/{wildcards.sample}*.csi - """ - -rule sample_stats: - input: - "{sample}.bam" - output: - temp("{sample}.bam.bai"), - stats = "reports/data/samtools_stats/{sample}.filtered.stats" - log: - "logs/stats/{sample}.stats.log" - threads: - 2 - shell: - """ - {{ - samtools index {input} - samtools stats -@ 1 -x -d {input} > {output.stats} - }} 2> {log} - """ - -rule molecule_coverage: - input: - fai = f"{workflow_geno}.fai", - stats = "reports/data/lrstats/{sample}.lrstats.gz" - output: - "reports/data/coverage/{sample}.molcov.gz" - log: - "logs/stats/{sample}.molstats.log" - params: - windowsize - shell: - "harpy-utils molecule-coverage -w {params} {input} 2> {log} | gzip > {output}" - -rule molecule_stats: - input: - "{sample}.bam" - output: - "reports/data/lrstats/{sample}.lrstats.gz" - log: - "logs/molcov/{sample}.molcov.log" - params: - molecule_distance - shell: - "harpy-utils bx-stats-sam -d {params} {input} 2> {log} | gzip > {output}" - -rule alignment_report: - input: - collect("reports/data/markdup/{sample}.markdup", sample = samplenames), - collect("reports/data/samtools_stats/{sample}.{data}.stats", sample = samplenames, data = ["raw", "filtered"]), - ipynb = f"workflow/samtools_stats.ipynb" - output: - tmp = temp("reports/bwa.summary.tmp.ipynb"), - ipynb = "reports/bwa.summary.ipynb" - params: - indir = "-p indir " + os.path.abspath("reports/data") - log: - f"logs/reports/bwa.report.log" - shell: - """ - export IPYTHONDIR=/tmp/ipython-bwa-stats - {{ - papermill -k xpython --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params.indir} - harpy-utils process-notebook {output.tmp} > {output.ipynb} - }} 2> {log} - """ - -rule sample_reports: - input: - lrstats = "reports/data/lrstats/{sample}.lrstats.gz", - coverage = "reports/data/coverage/{sample}.regions.bed.gz", - molcov = "reports/data/coverage/{sample}.molcov.gz", - ipynb = f"workflow/align_stats.ipynb" - output: - tmp = temp("reports/{sample}.tmp.ipynb"), - ipynb = "reports/{sample}.ipynb" - params: - lr_type = lr_type, - basedir = "-p basedir " + os.path.abspath("reports/data"), - mol_dist = f"-p mol_dist {molecule_distance}", - window_size = f"-p windowsize {windowsize}", - samplename = lambda wc: "-p samplename " + wc.get("sample"), - log: - "logs/reports/{sample}.report.log" - shell: - """ - export IPYTHONDIR=/tmp/ipython-{wildcards.sample}.rpt - {{ - papermill -k xpython --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} -p platform {params} - harpy-utils process-notebook {output.tmp} {wildcards.sample} BWA-MEM2 {params.lr_type} > {output.ipynb} - }} 2> {log} - """ - -rule linked_read_report: - input: - collect("reports/data/lrstats/{sample}.lrstats.gz", sample = samplenames), - ipynb = f"workflow/align_lrstats.ipynb" - output: - tmp = temp("reports/linkedreads.summary.tmp.ipynb"), - ipynb = "reports/linkedreads.summary.ipynb" - params: - lr_type = lr_type, - indir = "-p indir " + os.path.abspath("reports/data/lrstats") - log: - f"logs/reports/lrstats.report.log" - shell: - """ - {{ - export IPYTHONDIR=/tmp/ipython-bwa.lr - papermill -k xpython --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params.indir} - harpy-utils process-notebook {output.tmp} {params.lr_type} > {output.ipynb} - }} 2> {log} - """ - -rule all: - default_target: True - input: - bams = collect("{sample}.bam", sample = samplenames), - reports = collect("reports/{sample}.ipynb", sample = samplenames) if not skip_reports and not ignore_bx else [], - align_report = "reports/bwa.summary.ipynb" if (not skip_reports and len(samplenames) > 1) else [], - bx_report = "reports/linkedreads.summary.ipynb" if (not skip_reports and not ignore_bx and len(samplenames) > 1) else [] diff --git a/harpy/snakefiles/align_minimap.smk b/harpy/snakefiles/align_minimap.smk new file mode 100644 index 000000000..659534f60 --- /dev/null +++ b/harpy/snakefiles/align_minimap.smk @@ -0,0 +1,97 @@ +import os +import re + +wildcard_constraints: + sample = r"[a-zA-Z0-9._-]+" + +WORKFLOW = config.get('Workflow') or {} +PARAMETERS = config.get('Parameters') or {} +REPORTS = WORKFLOW.get("reports") or {} +INPUTS = config['Inputs'] +VERSION = WORKFLOW.get('harpy-version', 'latest') + +technology = PARAMETERS.get("aligner-technology", "sr") +illumina_old = PARAMETERS.get("illumina-format-old", False) +extra = PARAMETERS.get("extra", "") +fqlist = INPUTS["fastq"] +genomefile = INPUTS["reference"] +tech_opt = f"-ax map-{technology}" if technology != "sr" else "-ax map sr" + +bn = os.path.basename(genomefile) +workflow_geno = f"workflow/reference/{bn}" +genome_zip = True if bn.lower().endswith(".gz") else False +geno_idx = f"{workflow_geno}.gzi" if genome_zip else f"{workflow_geno}.fai" +bn_r = r"([_\.][12]|[_\.][FR]|[_\.]R[12](?:\_00[0-9])*)?\.((fastq|fq)(\.gz)?)$" +samplenames = {re.sub(bn_r, "", os.path.basename(i), flags = re.IGNORECASE) for i in fqlist} +d = dict(zip(samplenames, samplenames)) + +def get_fq(wildcards): + # returns a list of fastq files for read 1 based on *wildcards.sample* e.g. + r = re.compile(fr".*/({re.escape(wildcards.sample)}){bn_r}", flags = re.IGNORECASE) + return sorted(list(filter(r.match, fqlist))[:2]) + +rule process_reference: + input: + genomefile + output: + geno = workflow_geno, + idx = multiext(workflow_geno, ".mmi"), + fai = f"{workflow_geno}.fai", + gzi = f"{workflow_geno}.gzi" if genome_zip else [] + log: + f"{workflow_geno}.preprocess.log" + params: + genome_zip + conda: + "envs/align.yaml" + container: + f"docker://pdimens/harpy:align_{VERSION}" + shell: + """ + {{ + if (file {input} | grep -q compressed ) ;then + # is regular gzipped, needs to be BGzipped + seqtk seq {input} | bgzip -c > {output.geno} + else + cp -f {input} {output.geno} + fi + + if [ "{params}" = "True" ]; then + samtools faidx --gzi-idx {output.gzi} --fai-idx {output.fai} {output.geno} + else + samtools faidx --fai-idx {output.fai} {output.geno} + fi + minimap2 -d {output.idx} {output.geno} + }} 2> {log} + """ + +rule align: + input: + ref = workflow_geno + ".mmi", + fastq = get_fq + output: + bam = temp("minimap/{sample}.minimap.bam"), + tmp = temp(directory("minimap/{sample}_tmp")) + log: + "logs/minimap/{sample}.minimap.log" + params: + RG_tag = lambda wc: "-R \"@RG\\tID:" + wc.get("sample") + "\\tSM:" + wc.get("sample") + "\"", + tech = tech_opt, + static = "--MD -y" if illumina_old else "--MD", + extra = extra + threads: + 12 + resources: + tmpdir = lambda wc: f"minimap/{wc.sample}_tmp" + conda: + "envs/align.yaml" + container: + f"docker://pdimens/harpy:align_{VERSION}" + shell: + """ + mkdir -p {resources.tmpdir} + {{ + minimap2 -t {threads} {params} {input} | + samtools collate -T {resources.tmpdir} -O -u - + }} 2> {log} > {output.bam} + """ diff --git a/harpy/snakefiles/align_strobe.smk b/harpy/snakefiles/align_strobe.smk index 1c86c0217..210d89677 100644 --- a/harpy/snakefiles/align_strobe.smk +++ b/harpy/snakefiles/align_strobe.smk @@ -1,7 +1,6 @@ import os import re -localrules: all wildcard_constraints: sample = r"[a-zA-Z0-9._-]+" @@ -11,22 +10,14 @@ REPORTS = WORKFLOW.get("reports") or {} INPUTS = config['Inputs'] VERSION = WORKFLOW.get('harpy-version', 'latest') -lr_type = WORKFLOW.get("linkedreads", {}).get("type", 'none') -bx_tag = WORKFLOW.get("linkedreads", {}).get("standardized", {}).get("BX", False) -vx_tag = WORKFLOW.get("linkedreads", {}).get("standardized", {}).get("VX", False) -skip_reports = REPORTS.get("skip", False) illumina_old = PARAMETERS.get("illumina-format-old", False) -windowsize = PARAMETERS.get("depth-windowsize", 50000) -molecule_distance = PARAMETERS.get("distance-threshold", 0) -keep_unmapped = PARAMETERS.get("keep-unmapped", False) extra = PARAMETERS.get("extra", "") fqlist = INPUTS["fastq"] genomefile = INPUTS["reference"] bn = os.path.basename(genomefile) -bn_r = r"([_\.][12]|[_\.][FR]|[_\.]R[12](?:\_00[0-9])*)?\.((fastq|fq)(\.gz)?)$" -ignore_bx = lr_type == "none" bn = bn[:-3] if bn.lower().endswith(".gz") else bn +bn_r = r"([_\.][12]|[_\.][FR]|[_\.]R[12](?:\_00[0-9])*)?\.((fastq|fq)(\.gz)?)$" workflow_geno = f"workflow/reference/{bn}" samplenames = {re.sub(bn_r, "", os.path.basename(i), flags = re.IGNORECASE) for i in fqlist} d = dict(zip(samplenames, samplenames)) @@ -57,10 +48,10 @@ rule align: workflow_geno, get_fq, output: - bam = temp("strobealign/{sample}.strobe.bam"), - tmp = temp(directory("strobealign/{sample}_tmp")) + bam = temp("strobe/{sample}.strobe.bam"), + tmp = temp(directory("strobe/{sample}_tmp")) log: - "logs/strobealign/{sample}.strobealign.log" + "logs/strobealign/{sample}.strobe.log" params: static = "-N 2 -C" if illumina_old else "-N 2", RGid = lambda wc: f"--rg-id={wc.get('sample')}", @@ -69,7 +60,7 @@ rule align: threads: 12 resources: - tmpdir = lambda wc: f"strobealign/{wc.sample}_tmp" + tmpdir = lambda wc: f"strobe/{wc.sample}_tmp" conda: "envs/align.yaml" container: @@ -82,217 +73,3 @@ rule align: samtools collate -T {resources.tmpdir} -O -u -l 0 - }} 2> {log} > {output.bam} """ - -rule sort: - retries: 3 - input: - ref = workflow_geno, - bam = "strobealign/{sample}.strobe.bam" - output: - bam = temp("sort/{sample}.sort.bam"), - stats = "reports/data/samtools_stats/{sample}.raw.stats", - tmp = temp(directory("sort/{sample}_tmp")) - log: - "logs/sort/{sample}.sort.log" - params: - sortthreads = lambda wc, threads: threads - 1 - threads: - 4 - resources: - tmpdir = lambda wc: f"sort/{wc.sample}_tmp", - mem_mb_per_thread = lambda wc, attempt: 3000 // attempt - shell: - """ - mkdir -p {resources.tmpdir} - {{ - samtools fixmate -z on -m -u {input.bam} - | - samtools sort -@ {params.sortthreads} -M -T {resources.tmpdir} -o {output.bam} -u -l 0 -m {resources.mem_mb_per_thread}M - - samtools stats -@ {params.sortthreads} -d -x -r {input.ref} {output.bam} > {output.stats} - }} 2> {log} - """ - -rule mark_duplicates: - priority: 1 - input: - fq = get_fq, - bam = "sort/{sample}.sort.bam" - output: - bam = "{sample}.bam" if lr_type == "none" or (bx_tag and vx_tag) else temp("markdup/{sample}.bam"), - stats = "reports/data/markdup/{sample}.markdup", - tmp = temp(directory("markdup/{sample}_tmp")) - log: - "logs/markdup/{sample}.markdup.log" - params: - bx_mode = "-S --barcode-tag BX" if not ignore_bx else "-S", - quality = PARAMETERS.get('min-map-quality', 30), - unmapped = "-F 4" if not keep_unmapped else "", - mdthreads = lambda wc, threads: threads - 1 - resources: - tmpdir = lambda wc: f"markdup/{wc.sample}_tmp" - threads: - 4 - shell: - """ - mkdir -p {resources.tmpdir} - OPT=$(harpy-utils optical-dist-fq {input.fq}) - {{ - samtools view -h -u -q {params.quality} {params.unmapped} {input.bam} | - samtools markdup -@ {params.mdthreads} -T {resources.tmpdir} {params.bx_mode} -d $OPT -f {output.stats} - {output.bam} - }} 2> {log} - """ - -if lr_type != "none" and not (bx_tag and vx_tag): - rule standardize: - input: - "markdup/{sample}.bam" - output: - "{sample}.bam" - log: - "logs/{sample}.std.log" - threads: - 2 - shell: - "djinn-standardize --threads {threads} {input} > {output} 2> {log}" - -rule depth_stats: - input: - "{sample}.bam.bai", - bam = "{sample}.bam" - output: - "reports/data/coverage/{sample}.regions.bed.gz" - params: - f"-b {windowsize}", - "-n --fast-mode" - log: - "logs/depthstats/{sample}.mosdepth.log" - threads: - 2 - conda: - "envs/qc.yaml" - container: - f"docker://pdimens/harpy:qc_{VERSION}" - shell: - """ - mosdepth {params} -t 1 reports/data/coverage/{wildcards.sample} {input.bam} 2> {log} - rm -f reports/data/coverage/{wildcards.sample}.mosdepth* reports/data/coverage/{wildcards.sample}*.csi - """ - -rule sample_stats: - input: - "{sample}.bam" - output: - temp("{sample}.bam.bai"), - stats = "reports/data/samtools_stats/{sample}.filtered.stats" - log: - "logs/stats/{sample}.stats.log" - threads: - 2 - shell: - """ - {{ - samtools index {input} - samtools stats -@ 1 -x -d {input} > {output.stats} - }} 2> {log} - """ - -rule molecule_coverage: - input: - fai = f"{workflow_geno}.fai", - stats = "reports/data/lrstats/{sample}.lrstats.gz" - output: - "reports/data/coverage/{sample}.molcov.gz" - log: - "logs/stats/{sample}.molstats.log" - params: - windowsize - shell: - "harpy-utils molecule-coverage -w {params} {input} 2> {log} | gzip > {output}" - -rule molecule_stats: - input: - "{sample}.bam" - output: - "reports/data/lrstats/{sample}.lrstats.gz" - log: - "logs/molcov/{sample}.molcov.log" - params: - molecule_distance - shell: - "harpy-utils bx-stats-sam -d {params} {input} 2> {log} | gzip > {output}" - -rule alignment_report: - input: - collect("reports/data/samtools_stats/{sample}.{data}.stats", sample = samplenames, data = ["raw","filtered"]), - collect("reports/data/markdup/{sample}.markdup", sample = samplenames), - ipynb = f"workflow/samtools_stats.ipynb" - output: - tmp = temp("reports/strobealign.summary.tmp.ipynb"), - ipynb = "reports/strobealign.summary.ipynb" - params: - lr_type = lr_type, - indir = "-p indir " + os.path.abspath("reports/data") - log: - f"logs/reports/strobealign.report.log" - shell: - """ - export IPYTHONDIR=/tmp/ipython-stobe-stats - {{ - papermill -k xpython --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params.indir} - harpy-utils process-notebook {output.tmp} {params.lr_type} > {output.ipynb} - }} 2> {log} - """ - -rule sample_reports: - input: - lrstats = "reports/data/lrstats/{sample}.lrstats.gz", - coverage = "reports/data/coverage/{sample}.regions.bed.gz", - molecule_coverage = "reports/data/coverage/{sample}.molcov.gz", - ipynb = f"workflow/align_stats.ipynb" - output: - tmp = temp("reports/{sample}.tmp.ipynb"), - ipynb = "reports/{sample}.ipynb" - params: - lr_type = lr_type, - basedir = "-p basedir " + os.path.abspath("reports/data"), - mol_dist = f"-p mol_dist {molecule_distance}", - window_size = f"-p windowsize {windowsize}", - samplename = lambda wc: "-p samplename " + wc.get("sample"), - log: - "logs/reports/{sample}.report.log" - shell: - """ - export IPYTHONDIR=/tmp/ipython-{wildcards.sample}.rpt - {{ - papermill -k xpython --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} -p platform {params} - harpy-utils process-notebook {output.tmp} {wildcards.sample} strobealign {params.lr_type} > {output.ipynb} - }} 2> {log} - """ - -rule linked_read_report: - input: - collect("reports/data/lrstats/{sample}.lrstats.gz", sample = samplenames), - ipynb = f"workflow/align_lrstats.ipynb" - output: - tmp = temp("reports/linkedreads.summary.tmp.ipynb"), - ipynb = "reports/linkedreads.summary.ipynb" - params: - lr_type = lr_type, - indir = "-p indir " + os.path.abspath("reports/data/lrstats") - log: - f"logs/reports/lrstats.report.log" - shell: - """ - export IPYTHONDIR=/tmp/ipython-strobe.lr - {{ - papermill -k xpython --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params.indir} - harpy-utils process-notebook {output.tmp} {params.lr_type} > {output.ipynb} - }} 2> {log} - """ - -rule all: - default_target: True - input: - bams = collect("{sample}.bam", sample = samplenames), - reports = collect("reports/{sample}.ipynb", sample = samplenames) if not skip_reports and not ignore_bx else [], - align_report = "reports/strobealign.summary.ipynb" if (not skip_reports and len(samplenames) > 1) else [], - bx_report = "reports/linkedreads.summary.ipynb" if (not skip_reports and not ignore_bx and len(samplenames) > 1) else [] diff --git a/harpy/snakefiles/impute.smk b/harpy/snakefiles/impute.smk index b0ae4b0d8..fabfaebce 100644 --- a/harpy/snakefiles/impute.smk +++ b/harpy/snakefiles/impute.smk @@ -262,7 +262,7 @@ rule contig_report: export IPYTHONDIR=/tmp/ipython-{wildcards.paramset}.{wildcards.contig} {{ bcftools stats -s "-" {input.vcf} > {output.stats} - papermill -k xpython --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params} + papermill -k ipython-harpy --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params} harpy-utils process-notebook {output.tmp} {wildcards.contig} {wildcards.paramset} > {output.ipynb} }} 2> {log} """ @@ -296,7 +296,7 @@ rule impute_reports: {{ bcftools stats -s "-" {input.orig} {input.impute} | grep \"GCTs\" > {output.comparison} bcftools query -f '%CHROM\\t%POS\\t%INFO/INFO_SCORE\\n' {input.impute} > {output.infoscore} - papermill -k xpython --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params} + papermill -k ipython-harpy --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params} harpy-utils process-notebook {output.tmp} {wildcards.paramset} > {output.ipynb} }} 2> {log} """ diff --git a/harpy/snakefiles/phase_snp.smk b/harpy/snakefiles/phase_snp.smk index 912e9a806..5b71551d8 100644 --- a/harpy/snakefiles/phase_snp.smk +++ b/harpy/snakefiles/phase_snp.smk @@ -257,7 +257,7 @@ rule phase_report: """ export IPYTHONDIR=/tmp/ipython-phase {{ - papermill -k xpython --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params} + papermill -k ipython-harpy --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params} harpy-utils process-notebook {output.tmp} > {output.ipynb} }} 2> {log} """ diff --git a/harpy/snakefiles/preprocess_gih.smk b/harpy/snakefiles/preprocess_gih.smk index e22973688..8e6f0a6bc 100644 --- a/harpy/snakefiles/preprocess_gih.smk +++ b/harpy/snakefiles/preprocess_gih.smk @@ -179,7 +179,7 @@ rule barcode_report: """ export IPYTHONDIR=/tmp/ipython-pre-gih {{ - papermill -k xpython --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params.indir} + papermill -k ipython-harpy --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params.indir} harpy-utils process-notebook {output.tmp} > {output.ipynb} }} 2> {log} """ diff --git a/harpy/snakefiles/qc.smk b/harpy/snakefiles/qc.smk index bfc6f169b..35c641b58 100644 --- a/harpy/snakefiles/qc.smk +++ b/harpy/snakefiles/qc.smk @@ -98,7 +98,7 @@ rule barcode_report: """ export IPYTHONDIR=/tmp/ipython-lrstats {{ - papermill -k xpython --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params.indir} + papermill -k ipython-harpy --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params.indir} harpy-utils process-notebook {output.tmp} {params.lr} > {output.ipynb} }} 2> {log} """ @@ -118,7 +118,7 @@ rule qc_report: """ export IPYTHONDIR=/tmp/ipython-fastp {{ - papermill -k xpython --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params} + papermill -k ipython-harpy --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params} harpy-utils process-notebook {output.tmp} > {output.ipynb} }} 2> {log} """ diff --git a/harpy/snakefiles/snp_freebayes.smk b/harpy/snakefiles/snp_freebayes.smk index a6859c538..329ba7eed 100644 --- a/harpy/snakefiles/snp_freebayes.smk +++ b/harpy/snakefiles/snp_freebayes.smk @@ -169,8 +169,8 @@ rule variant_report: export IPYTHONDIR=/tmp/ipython-snp-{wildcards.type} {{ bcftools stats -s "-" --fasta-ref {input.genome} {input.bcf} > {output.data} - papermill -k xpython --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params} - harpy-utils process-notebook {output.tmp} variants.{wildcards.type} > {output.ipynb} + papermill -k ipython-harpy --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params} + harpy-utils process-notebook {output.tmp} "Variants ({wildcards.type})" > {output.ipynb} }} 2> {log} """ diff --git a/harpy/snakefiles/snp_mpileup.smk b/harpy/snakefiles/snp_mpileup.smk index 3ddaa23ba..fe2f4972b 100644 --- a/harpy/snakefiles/snp_mpileup.smk +++ b/harpy/snakefiles/snp_mpileup.smk @@ -196,8 +196,8 @@ rule variant_report: export IPYTHONDIR=/tmp/ipython.snp.{wildcards.type} {{ bcftools stats -s "-" --fasta-ref {input.genome} {input.bcf} > {output.data} - papermill -k xpython --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params} - harpy-utils process-notebook {output.tmp} variants.{wildcards.type} > {output.ipynb} + papermill -k ipython-harpy --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params} + harpy-utils process-notebook {output.tmp} "Variants ({wildcards.type})" > {output.ipynb} }} 2> {log} """ diff --git a/harpy/snakefiles/sv_leviathan.smk b/harpy/snakefiles/sv_leviathan.smk index bcfbb0b57..b291c437a 100644 --- a/harpy/snakefiles/sv_leviathan.smk +++ b/harpy/snakefiles/sv_leviathan.smk @@ -232,7 +232,7 @@ rule report: """ export IPYTHONDIR=/tmp/ipython-leviathan {{ - papermill -k xpython --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params} + papermill -k ipython-harpy --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params} harpy-utils process-notebook {output.tmp} LEVIATHAN > {output.ipynb} }} 2> {log} """ diff --git a/harpy/snakefiles/sv_naibr.smk b/harpy/snakefiles/sv_naibr.smk index bb6143b76..31ab7c7ff 100644 --- a/harpy/snakefiles/sv_naibr.smk +++ b/harpy/snakefiles/sv_naibr.smk @@ -205,7 +205,7 @@ rule report: """ export IPYTHONDIR=/tmp/ipython-sv.naibr {{ - papermill -k xpython --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params} + papermill -k ipython-harpy --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params} harpy-utils process-notebook {output.tmp} NAIBR > {output.ipynb} }} 2> {log} """ diff --git a/harpy/snakefiles/validate_bam.smk b/harpy/snakefiles/validate_bam.smk index dd381a906..4ec7e3856 100644 --- a/harpy/snakefiles/validate_bam.smk +++ b/harpy/snakefiles/validate_bam.smk @@ -57,7 +57,7 @@ rule create_report: """ export IPYTHONDIR=/tmp/ipython-validate-xam {{ - papermill -k xpython --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params.infile} + papermill -k ipython-harpy --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params.infile} harpy-utils process-notebook {output.tmp} {params.lr_platform} > {output.ipynb} }} 2> {log} """ diff --git a/harpy/snakefiles/validate_fastq.smk b/harpy/snakefiles/validate_fastq.smk index 623c7593f..a04e5574b 100644 --- a/harpy/snakefiles/validate_fastq.smk +++ b/harpy/snakefiles/validate_fastq.smk @@ -72,7 +72,7 @@ rule create_report: """ export IPYTHONDIR=/tmp/ipython-validate-fastq {{ - papermill -k xpython --cwd . --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params.infile} + papermill -k ipython-harpy --cwd . --no-progress-bar --log-level ERROR {input.ipynb} {output.tmp} {params.infile} harpy-utils process-notebook {output.tmp} {params.lr_platform} > {output.ipynb} }} 2> {log} """ diff --git a/harpy/utils/check_fastq.py b/harpy/utils/check_fastq.py index 3ea709a9e..371cd5b1a 100755 --- a/harpy/utils/check_fastq.py +++ b/harpy/utils/check_fastq.py @@ -8,7 +8,7 @@ class FormatChecker: """Base class: shared counters + SAM-spec/BX-position validation.""" - SAMSPEC = re.compile(r'[A-Z][A-Z]:[AifZHB]:') + SAMSPEC = re.compile(r'[A-Za-z][A-Za-z0-9]:[AifZHB]:') def __init__(self): self.N_READS = 0 self.NO_BX = 0 diff --git a/harpy/utils/standardize/standardize.go b/harpy/utils/standardize/standardize.go index 2b3d62ed4..bb3cdb376 100644 --- a/harpy/utils/standardize/standardize.go +++ b/harpy/utils/standardize/standardize.go @@ -9,7 +9,6 @@ import ( "os" "regexp" "runtime" - "strconv" "github.com/biogo/hts/bam" "github.com/biogo/hts/sam" @@ -65,44 +64,6 @@ func GetStringTag(r *sam.Record, tag string) (string, bool) { return "", false } -func copyHeaderWithPG(src *sam.Header, infile string) (*sam.Header, error) { - // Clone the header via marshal/unmarshal - b, err := src.MarshalText() - if err != nil { - return nil, err - } - dst := &sam.Header{} - if err := dst.UnmarshalText(b); err != nil { - return nil, err - } - - // Build the @PG record - pg := sam.NewProgram( - "djinn", // ID - "djinn", // name (PN) - "djinn standardize "+infile, // command line (CL) - lastPGID(src), // previous PG ID (PP), or "" if none - "1.0", // version (VN) — set as appropriate - ) - - if err := dst.AddProgram(pg); err != nil { - return nil, err - } - - return dst, nil -} - -// lastPGID returns the ID of the last @PG record in the header, -// which becomes the PP (previous program) of the new entry. -// Returns "" if there are no existing @PG records. -func lastPGID(h *sam.Header) string { - progs := h.Progs() - if len(progs) == 0 { - return "" - } - return strconv.Itoa(progs[len(progs)-1].ID()) -} - func main() { nThreads := flag.Int("threads", 1, "Number of threads for BAM io. No real value beyond ~4-5.") flag.Usage = func() { @@ -139,14 +100,28 @@ func main() { defer br.Close() // ── BAM writer ────────────────────────────────────────────────────────────── - newHeader, err := copyHeaderWithPG(br.Header(), infile) - if err != nil { + hdr := br.Header() + progs := hdr.Progs() + + var prev string + if len(progs) > 0 { + prev = progs[len(progs)-1].UID() + } + + pg := sam.NewProgram( + "djinn", // ID + "djinn", // name (PN) + "djinn standardize "+infile, // command line (CL) + prev, // previous PG ID (PP), or "" if none + "1.0", // version (VN) — set as appropriate + ) + + if err := hdr.AddProgram(pg); err != nil { log.Fatal(err) } // WRITES SAM OUTPUT - //w, err := sam.NewWriter(os.Stdout, newHeader, sam.FlagDecimal) - w, err := bam.NewWriterLevel(os.Stdout, newHeader, 4, writeThread) + w, err := bam.NewWriterLevel(os.Stdout, hdr, 4, writeThread) if err != nil { log.Fatal(err) } diff --git a/harpy/validation/fasta.py b/harpy/validation/fasta.py index 3d84d9679..99a85b853 100644 --- a/harpy/validation/fasta.py +++ b/harpy/validation/fasta.py @@ -1,11 +1,13 @@ import os +import re import pysam -from harpy.common.file_ops import safe_read from harpy.common.printing import HarpyPrint +alphanum = re.compile(r'^[a-zA-Z0-9_\-\|]+$') +nuc = re.compile(r'[ACGTURYKMSWBDHVN\-]+$', flags = re.IGNORECASE) class FASTA(): ''' @@ -13,68 +15,15 @@ class FASTA(): ''' def __init__(self, fasta, quiet:int = 0): self.file = fasta + self.file_base = os.path.basename(self.file) self.print = HarpyPrint(quiet) - self.print.log("Correct FASTA file format", newline = False) - # validate fasta file contents - line_num = 0 - seq_id = 0 - seq = 0 - last_header = False - with safe_read(self.file) as fasta: - for line in fasta: - line_num += 1 - if line.startswith(">"): - seq_id += 1 - if last_header: - self.print.validation(False) - self.print.error( - "consecutive contig names", - f"All contig names must be followed by at least one line of nucleotide sequences, but two consecutive lines of contig names were detected. This issue was identified at line [bold]{line_num}[/] in [blue]{self.file}[/], but there may be others further in the file.", - "See the FASTA file spec and try again after making the appropriate changes: https://www.ncbi.nlm.nih.gov/genbank/fastaformat/" - ) - else: - last_header = True - if len(line.rstrip()) == 1: - self.print.validation(False) - self.print.error( - "unnamed contigs", - f"All contigs must have an alphanumeric name, but a contig was detected without a name. This issue was identified at line [bold]{line_num}[/] in [blue]{self.file}[/], but there may be others further in the file.", - "See the FASTA file spec and try again after making the appropriate changes: https://www.ncbi.nlm.nih.gov/genbank/fastaformat/" - ) - if line.startswith("> "): - self.print.validation(False) - self.print.error( - "invalid contig names", - f"All contig names must be named [green bold]>contig_name[/], without a space, but a contig was detected with a space between the [green bold]>[/] and contig_name. This issue was identified at line [bold]{line_num}[/] in [blue]{self.file}[/], but there may be others further in the file.", - "See the FASTA file spec and try again after making the appropriate changes: https://www.ncbi.nlm.nih.gov/genbank/fastaformat/" - ) - elif line == "\n": - self.print.validation(False) - self.print.error( - "empty lines", - f"Empty lines are not permitted in FASTA files, but one was detected at line [bold]{line_num}[/] in [blue]{self.file}[/]. The scan ended at this error, but there may be others further in the file.", - "See the FASTA file spec and try again after making the appropriate changes: https://www.ncbi.nlm.nih.gov/genbank/fastaformat/" - ) - else: - seq += 1 - last_header = False - solutiontext = "FASTA files must have at least one contig name followed by sequence data on the next line. Example:\n" - solutiontext += "[green] >contig_name\n ATACAGGAGATTAGGCA[/]\n" - # make sure there is at least one of each - if seq_id == 0: - self.print.validation(False) - self.print.error( - "contig names absent", - f"No contig names detected in [blue]{self.file}[/].", - f"{solutiontext}\nSee the FASTA file spec and try again after making the appropriate changes: https://www.ncbi.nlm.nih.gov/genbank/fastaformat/" - ) - if seq == 0: + self.print.log("Validating input FASTA file", newline=False) + try: + with pysam.FastxFile(self.file, persist=False) as fa: + pass + except Exception as e: self.print.validation(False) - self.print.error( - "sequences absent", - f"No sequences detected in [blue]{self.file}[/].", - f"{solutiontext}\nSee the FASTA file spec and try again after making the appropriate changes: https://www.ncbi.nlm.nih.gov/genbank/fastaformat/" - ) + self.print.error("bad FASTA file", e) self.print.validation(True) @@ -90,13 +39,12 @@ def match_contigs(self, contigs: str): if i not in valid_contigs: bad_names.append(i) if bad_names: - shortname = os.path.basename(self.file) self.print.validation(False) self.print.error( "contigs absent", - f"Some of the provided contigs were not found in [blue]{shortname}[/]. This will definitely cause plotting errors in the workflow.", + f"Some of the provided contigs were not found in [blue]{self.file_base}[/]. This will definitely cause plotting errors in the workflow.", "Check that your contig names are correct, including uppercase and lowercase.", - f"Contigs absent in {shortname}", + f"Contigs absent in {self.file_base}", ",".join([i for i in bad_names]) ) self.print.validation(True) @@ -150,8 +98,8 @@ def validate_region(self, regioninput) -> None: if row[0] not in contigs: self.print.error( "missing contig", - f"The contig listed at row {idx} ([bold yellow]{row[0]}[/]) is not present in ([blue]{os.path.basename(self.file)}[/]). This is the first row triggering this error, but it may not be the only one.", - f"Check that all the contigs listed in [blue]{os.path.basename(regioninput)}[/] are also present in [blue]{os.path.basename(self.file)}[/]", + f"The contig listed at row {idx} ([bold yellow]{row[0]}[/]) is not present in ([blue]{self.file_base}[/]). This is the first row triggering this error, but it may not be the only one.", + f"Check that all the contigs listed in [blue]{os.path.basename(regioninput)}[/] are also present in [blue]{self.file_base}[/]", "Row triggering this error", line ) diff --git a/harpy/validation/vcf.py b/harpy/validation/vcf.py index fe341144f..ae74f74ab 100644 --- a/harpy/validation/vcf.py +++ b/harpy/validation/vcf.py @@ -13,19 +13,25 @@ class VCF(): ''' A class to contain and validate a VCF input file. ''' - def __init__(self, filename:str, workdir:str, quiet:int = 0): + def __init__(self, filename:str, workdir:str, quiet:int = 0, threads: int = 2): os.makedirs(workdir, exist_ok = True) self.file: str = filename + self.file_base: str = os.path.basename(filename) self.workdir: str = workdir self.biallelic_file: str = "" self.contigs: dict[str,int] = {} self.print = HarpyPrint(quiet) + self.threads = threads - 1 if self.file.lower().endswith("bcf") and not os.path.exists(f"{self.file}.csi"): - pysam.bcftools.index(self.file) + self.print.log("Indexing BCF" , newline=False) + pysam.bcftools.index("--threads", str(self.threads), self.file) + self.print.validation(True) if self.file.lower().endswith("vcf.gz") and not os.path.exists(f"{self.file}.tbi"): - pysam.bcftools.index("--tbi", self.file) + self.print.log("Indexing VCF" , newline=False) + pysam.bcftools.index("--threads", str(self.threads),"--tbi", self.file) + self.print.validation(True) def get_contigs(self): """reads the header of a vcf/bcf file and populate `self.contigs` with the contigs (keys) and their lengths (values)""" @@ -52,7 +58,7 @@ def find_biallelic_contigs(self): keep = False # Use bcftools to count the number of biallelic SNPs in the contig viewcmd = subprocess.Popen( - [bcftools, 'view', '-H', '-r', str(contig), '-v', 'snps', '-m2', '-M2', '-c', '2', self.file], + [bcftools, 'view', '-H', '--threads', str(self.threads), '-r', str(contig), '-v', 'snps', '-m2', '-M2', '-c', '2', self.file], stdout=subprocess.PIPE, text = True ) @@ -84,7 +90,7 @@ def find_biallelic_contigs(self): def check_phase(self): """Check to see if the input VCf file is phased or not, determined by the presence of ID=PS or ID=HP tags""" self.print.log("VCF file is phased ([green]PS[/] or [green]HP[/] tags)", newline=False) - with pysam.VariantFile(self.file) as _vcf: + with pysam.VariantFile(self.file, threads = self.threads) as _vcf: formats = list(_vcf.header.formats) if 'PS' not in formats and 'HP' not in formats: bn = os.path.basename(self.file) @@ -115,7 +121,7 @@ def match_samples(self, bamlist: list[str], prioritize_vcf: bool) -> None: self.print.error( "mismatched inputs", f"There are [bold]{len(missing_samples)}[/] samples found in [blue]{fromthis}[/] that are not in [blue]{inthis}[/]. Terminating Harpy to avoid downstream errors.", - f"[blue]{fromthis}[/] cannot contain samples that are absent in [blue]{inthis}[/]. Check the spelling or remove those samples from [blue]{fromthis}[/] or remake the vcf file to include/omit these samples. Alternatively, toggle [green]--vcf-samples[/] to aggregate the sample list from the input files or [blue]{self.file}[/].", + f"[blue]{fromthis}[/] cannot contain samples that are absent in [blue]{inthis}[/]. Check the spelling or remove those samples from [blue]{fromthis}[/] or remake the vcf file to include/omit these samples. Alternatively, toggle [green]--vcf-samples[/] to aggregate the sample list from the input files or [blue]{self.file_base}[/].", "The samples causing this error are", ", ".join(sorted(missing_samples)) + "\n" ) @@ -132,13 +138,12 @@ def match_contigs(self, contigs: list[str]): if i not in self.contigs: bad_names.append(i) if bad_names: - shortname = os.path.basename(self.file) self.print.validation(False) self.print.error( "contigs absent", - f"Some of the provided contigs were not found in [blue]{shortname}[/]. This will definitely cause plotting errors in the workflow.", + f"Some of the provided contigs were not found in [blue]{self.file_base}[/]. This will definitely cause plotting errors in the workflow.", "Check that your contig names are correct, including uppercase and lowercase. It's possible that you listed a contig in the genome that isn't in the variant call file due to filtering.", - f"Contigs absent in {shortname}", + f"Contigs absent in {self.file_base}", ",".join([i for i in bad_names]) ) self.print.validation(True) diff --git a/harpy/validation/xam.py b/harpy/validation/xam.py index e4db6643c..9c3ca19cf 100644 --- a/harpy/validation/xam.py +++ b/harpy/validation/xam.py @@ -43,7 +43,7 @@ def __init__(self, filenames, detect_bc:bool = False, nonlinked_ok:bool = True, uniqs.add(bn) self.count += 1 try: - with pysam.AlignmentFile(i, 'r', require_index=False): + with pysam.AlignmentFile(i, 'r', require_index=False, threads = 2): pass except (ValueError, OSError): badfiles.append(i) @@ -126,7 +126,7 @@ def __init__(self, filenames, detect_bc:bool = False, nonlinked_ok:bool = True, def is_phased(self, file_path: str) -> bool: ''' Scan the `file_path` to determine if the file has `PS` or `HP` tags''' - with pysam.AlignmentFile(file_path, require_index=False) as alnfile: + with pysam.AlignmentFile(file_path, require_index=False, threads = 2) as alnfile: for i, record in enumerate(alnfile.fetch(until_eof = True), 1): if i > 100: break @@ -139,7 +139,7 @@ def which_linkedread(self, file_path: str) -> str: Scans the first `self.max_records` records of a SAM/BAM file and tries to determine the barcode technology Returns one of: "haplotagging", "stlfr", "tellseq", or "none" """ - with pysam.AlignmentFile(file_path, require_index=False) as alnfile: + with pysam.AlignmentFile(file_path, require_index=False, threads = 2) as alnfile: for i, record in enumerate(alnfile.fetch(until_eof = True), 1): if i > 100: break diff --git a/pyproject.toml b/pyproject.toml index bd815eb28..2a8c84812 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -66,13 +66,14 @@ samtools = ">=1.23" seqtk = "*" snakemake-minimal = ">=9.19.0,<10" vl-convert-python = ">=2.0.0rc1,<3" -xeus-python = ">=0.18.1,<0.20" +ipykernel = ">=7.3.0,<8" [tool.pixi.pypi-dependencies] harpy = { path = ".", editable = true } [tool.pixi.activation] env = { JUPYTER_NOTARY_DB = ":memory:" } +scripts = ["resources/register-kernel.sh"] [tool.pixi.tasks] build-stagger = {cmd = "go build -C harpy/utils/stagger -o ../gih-stagger -ldflags='-s -w' stagger.go", inputs = ["harpy/utils/stagger/stagger.go"], outputs = ["harpy/utils/gih-stagger"]} diff --git a/resources/build.sh b/resources/build.sh index 743654523..a8c46e079 100644 --- a/resources/build.sh +++ b/resources/build.sh @@ -10,15 +10,16 @@ mv gih-stagger gih-convert djinn-standardize ${PREFIX}/bin/ } -## CLI completions +## activate/deactive processes mkdir -p $PREFIX/etc/conda/activate.d mkdir -p $PREFIX/etc/conda/deactivate.d -#echo "export JUPYTER_NOTARY_DB=':memory:'" > ${PREFIX}/etc/conda/activate.d/my-package-activate.sh -#echo "unset JUPYTER_NOTARY_DB" > ${PREFIX}/etc/conda/deactivate.d/my-package-deactivate.sh cat > ${PREFIX}/etc/conda/activate.d/harpy-activate.sh <<'EOF' export _HARPY_OLD_JUPYTER_NOTARY_DB="${JUPYTER_NOTARY_DB-__UNSET__}" -export JUPYTER_NOTARY_DB=':memory:' +export JUPYTER_NOTARY_DB=':memory:' + +python -m ipykernel install --prefix "$CONDA_PREFIX" --name ipython-harpy \ + --display-name "Python (harpy)" EOF cat > ${PREFIX}/etc/conda/deactivate.d/harpy-deactivate.sh <<'EOF' diff --git a/resources/meta.yaml b/resources/meta.yaml index 362e0b8e2..2190ad819 100644 --- a/resources/meta.yaml +++ b/resources/meta.yaml @@ -36,6 +36,7 @@ requirements: - conda >=24.8 - djinn >=2.3 - htslib >=1.23 + - ipykernel >=7.3 - jupyter-book >=2.1.0 - papermill >=2.6 - pandas >=2.3.3 @@ -48,7 +49,6 @@ requirements: - seqtk - snakemake-minimal >=9.19.0 - vl-convert-python >=2.0 - - xeus-python >=0.18.1 test: commands: @@ -60,7 +60,7 @@ about: license: "GPL-3.0-or-later" license_family: GPL3 license_file: LICENSE - summary: "Process raw haplotagging data, from raw sequences to phased haplotypes." + summary: "Process linked-read or WGS data, from raw sequences to phased haplotypes." description: | Harpy is a command-line tool to easily process platform-agnostic linked-read or WGS data. It uses Snakemake under the hood to execute different workflows (quality control, trimming, diff --git a/resources/register-kernel.sh b/resources/register-kernel.sh new file mode 100644 index 000000000..23f7ba086 --- /dev/null +++ b/resources/register-kernel.sh @@ -0,0 +1,4 @@ +#! /usr/bin/env bash + +python -m ipykernel install --prefix "${CONDA_PREFIX:?CONDA_PREFIX is required}" --name ipython-harpy \ + --display-name "Python (harpy)" \ No newline at end of file