Short read QC

[jmtsuji]

This code adds basic short read QC functionality to rotary to address #38. Summary of changes:

• A new master checkpoint for QC was created in rule qc. This selects (on a per-run basis, for the MVP version) whether short + long read QC will be performed or just long read QC
• Short read QC starts with tossing out junk reads (e.g., with ambiguous bases) and checking file integrity with rule short_read_reformat
• Adapter removal (from the right side of reads) and quality trimming are then performed in the same step, rule short_read_adapter_and_quality_trimming. Quality trimming is performed after adapter removal so that adapters near the ends of reads can be more accurately identified (e.g., before some bases are clipped off the ends). I think ATLAS takes the same approach.
  • TODO: the adapters file from ATLAS is directly downloaded from Zenodo - do we want to make and host our own file?
  • Note that adapter trimming and quality trimming are currently non-optional. We can add the ability to toggle these on/off if you would like.
  • I set a relatively aggressive q-trim threshold of Q20
• Contaminant filtration is the last step, via rule short_read_contamination_filter. This step is optional based on a variable in the config.
  • I built in PhiX genome and human genome contaminant filtration into rotary. These genomes will be auto-downloaded into the DB folder if the user requests to screen with them.
  • The user can optionally add a single additional FastA nucleotide file as a screening reference.

I did not add good logging for QC stats. The user needs to manually look at log files to understand what happened during QC. In future, it could be nice to add a rule that parses the log files to summarize QC results.

Lastly, the current code passes --dry-run on my machine but has only been tested in smaller chunks with real data. A proper end-to-end test will be needed at some point (maybe once everything is merged into #71 ?).

@LeeBergstrand Mind adding your review? It looks like there will be no merge conflicts with improve_multi_sample_pathing (the branch that this branch was derived from), which is good news! Thanks.

[LeeBergstrand]

@jmtsuji This looks good so far. I left a few comments.

[jmtsuji]

@LeeBergstrand Addressed your comments -- mind taking another look?

Also, what are your thoughts on the following (pasted from the PR description above)?

• TODO: the adapters file from ATLAS is directly downloaded from Zenodo - do we want to make and host our own file?
• Note that adapter trimming and quality trimming are currently non-optional. We can add the ability to toggle these on/off if you would like.
• I set a relatively aggressive q-trim threshold of Q20 [as default]

[LeeBergstrand]

@jmtsuji

• TODO: the adapters file from ATLAS is directly downloaded from Zenodo - do we want to make and host our own file?

We could put this down as a long-term TODO. I wonder if the archive will go down or if ATLAS will stop using this URL.

• Note that adapter trimming and quality trimming are currently non-optional. We can add the ability to toggle these on/off if you would like.

This should be implemented. However, I would like to merge the QC code into develop beforehand.

• I set a relatively aggressive q-trim threshold of Q20 [as default]

This is fine for Illumina, which should have reasonable error rates. We should have more stringent cutoffs since we are polishing with the short reads.

[jmtsuji]

@LeeBergstrand

We could put this down as a long-term TODO. I wonder if the archive will go down or if ATLAS will stop using this URL.

Addressed - the code now pulls from a rotary-specific data repo that contains the adapters file: https://zenodo.org/records/10087395 Attribution is given to the ATLAS repo.

[jmtsuji]

This should be implemented. However, I would like to merge the QC code into develop beforehand.

OK, sounds good.

This is fine for Illumina, which should have reasonable error rates. We should have more stringent cutoffs since we are polishing with the short reads.

Great, let's proceed with Q20 for now, then.

[jmtsuji]

Note that adapter trimming and quality trimming are currently non-optional. We can add the ability to toggle these on/off if you would like.

This should be implemented. However, I would like to merge the QC code into develop beforehand.

@LeeBergstrand FYI, I had a bit of extra time, so I went ahead and added the ability to toggle adapter and quality trimming on/off into the code. Basically:

• I had to split adapter trimming and quality trimming into 2 separate rules
• The adapter trimming rule does not run if the config flag perform_adapter_trimming is set to False
• The quality trimming rule always runs, but if quality_trim_direction is set to "f" in the config (an actual setting in bbduk), then the reads are just passed through the rule without actual trimming. This simplified the logic in the snakefile quite a bit.

[LeeBergstrand]

Improve by putting contaminant NCBI IDs in the config and also a list of paths to custom contaminant FASTA. As discussed below:

[jmtsuji]

Improve by putting contaminant NCBI IDs in the config and also a list of paths to custom contaminant FASTA.

Done! Let me know your thoughts on the updated code.

Aside from addressing the core changes you requested, some other things were also added:

• warning was added in the config file about ~60 GB memory needed for human genome-based decontamination
• pigz / unpigz is also forced on now

[LeeBergstrand]

@jmtsuji I'm working on an elegant solution to the short-read contaminant set up. It should address most of my review comments. See: https://github.com/jmtsuji/rotary/pull/82

[jmtsuji]

@LeeBergstrand Thanks for the proposed revision! I just added my comments at #82

[jmtsuji]

@LeeBergstrand I think this PR is only waiting on #82, correct? If so, then feel free to merge this PR after finishing the #82 merge. Otherwise, feel free to keep this PR open if you'd like to do another code review.