Picard: directly support Sentieon

[vladsavelyev]

Sentieon internally uses Picard tools to generate metrics. However, it adds a slightly different header into the outputs. Compare Picard's CollectGcBiasMetrics output:

## htsjdk.samtools.metrics.StringHeader
# picard.analysis.CollectGcBiasMetrics INPUT=SRR4238351_subsamp.sorted.bam <...>
## htsjdk.samtools.metrics.StringHeader
# Started on: Wed Apr 05 11:14:59 UTC 2017

## METRICS CLASS        picard.analysis.GcBiasDetailMetrics
ACCUMULATION_LEVEL      GC      WINDOWS READ_STARTS     MEAN_BASE_QUALITY       NORMALIZED_COVERAGE     ERROR_BAR_WIDTH SAMPLE  LIBRARY READ_GROUP
...

And Sentieon's --algo GCBias output:

#SentieonCommandLine: /usr/local/sentieon-genomics-201911/libexec/driver -t 36 -i sorted_0.bam -i sorted_1.bam -r genome/hs37d5.fa --algo GCBias --summary gc_summary.txt gc_metric.txt
ACCUMULATION_LEVEL     GC      WINDOWS READ_STARTS     MEAN_BASE_QUALITY       NORMALIZED_COVERAGE     ERROR_BAR_WIDTH SAMPLE  LIBRARY READ_GROUP
...

The module for Sentieon was implemented completely independent from Picard, however, the actual module code was copied from Picard almost identically (apart from treating the headers and extracting the sample name). It led to some overtime discrepancy, when Picard was updated, whereas Sentieon wasn't (e.g.).

Decoupling module dependencies is a great idea in general, but in my opinion that makes sense when modules exist for different tools, even if very similar. But because Sentieon literally uses the very same Picard tools, I think the code should always stay in sync, so I think it makes sense for the Sentieon module to just import the Picard module code. We already do a similar thing for Bracken that just imports Kraken.

Adds support for all Picard-compatible Sentieon QC outputs, fixes:

• Sentieon HsMetrics Module to match Picard HsMetrics #1661
• Request for support of Sentieon-Dedup #1936
• Picard outputs from Sentieon not parsed #1180 (merged, but not for all provided test examples)
• New module: Sentieon Dedup Metrics #1943 (Sentieon Dedup metrics)

UPD:

It makes sense to remove the Sentieon MultiQC module completely, and adjust Picard to support both header formats. This way, MultiQC won't break up the general stats table into two columns, and all the plots sections as well:

TODO:

• Check linting. We probably shouldn't demand to add doi and version in every submodule, but rather only in the main picard.py.

[ewels]

x-ref after discussion: #2152

[ewels]

Needs a special changelog entry warning Senteion users about changes:

• -m sentieon will no longer work
• Exported file names will have picard instead of senteion
• Report naming will change

[vladsavelyev]

@ewels would be great if you can take another look. I'd love to merge it before completing the ruff linting refactoring!

[ewels]

Feedback on the Slack thread in nf-core has been positive, so I think we're ok on that front 👍🏻

I'll have a quick look / play now.

[ewels]

In testing, I was looking for the Senteion sample NA12878 in the report generated from this PR. I didn't find it, but more worryingly I got a JS error (probably from another change?)

Uncaught TypeError: s.replace is not a function
    at Object.format (multiqc_report.html:1562:8475)
    at buildCache (multiqc_report.html:1561:2774)
    at HTMLTableElement.<anonymous> (multiqc_report.html:1562:5250)
    at Function.each (multiqc_report.html:1075:2815)
    at r.fn.init.each (multiqc_report.html:1075:1003)
    at r.fn.init.construct [as tablesorter] (multiqc_report.html:1562:4923)
    at HTMLDocument.<anonymous> (multiqc_report.html:1758:21)
    at j (multiqc_report.html:1075:29948)
    at k (multiqc_report.html:1075:30262)

Need to stop now, but are you able to replicate @vladsavelyev ? This is in a report generated from the picard directory in test data, with this PR checked out and installed. Using the highlighting toolbox functionality for NA12878.

[vladsavelyev]

NA12878 should be reported as sorted_0 - per the command line in the file header, e.g.:

#SentieonCommandLine: /usr/local/sentieon-genomics-201911/libexec/driver -t 36 -i sorted_0.bam -i sorted_1.bam -r genome/hs37d5.fa --algo InsertSizeMetricAlgo is_metric.txt

(in ./InsertSizeMetrics/sentieon_pr_1239/NA12878.InsertSize_metrics.txt)

[vladsavelyev]

Regarding the java script error - that's coming from a change in the previous PR.

Fixed it now: b88eb1b

[ewels]

So flicking between old / new reports before / after this PR, there is some discrepancies:

• ✅ 63 General stats rows goes to 57
  • Looking at a diff between the sample names column, I guess that this is all due to us now taking the input filename, so several Senteion samples are probably now overwriting one another as sorted_0 for example. Probably fine.
• ✅ 32 General stats columns to 29
  • Explained by the loss of the Senteion specific columns
• 🐛 In "Base Distribution Metrics" there are samples called None_R1 and None_R2 - sounds like we need a fallback there for not finding sample name from inputs?
• 🐛 Some unexpected fields in multiqc_data seem to be going missing
  • multiqc_picard_IlluminaLaneMetrics.yaml - loses fields PHASING_APPLIED, PREPHASING_APPLIED, TYPE_NAME?
  • multiqc_picard_ExtractIlluminaBarcodes.yaml - loses LANE

[vladsavelyev]

🐛 In "Base Distribution Metrics" there are samples called None_R1 and None_R2 - sounds like we need a fallback there for not finding sample name from inputs?

Fixed that, thanks!

🐛 Some unexpected fields in multiqc_data seem to be going missing
multiqc_picard_IlluminaLaneMetrics.yaml - loses fields PHASING_APPLIED, PREPHASING_APPLIED, TYPE_NAME?

There are two types of IlluminaLaneMetrics outputs, reporting different metrics:

## METRICS CLASS	picard.illumina.IlluminaPhasingMetrics
LANE	TYPE_NAME	PHASING_APPLIED	PREPHASING_APPLIED
...

and

## METRICS CLASS	picard.illumina.IlluminaLaneMetrics
CLUSTER_DENSITY	LANE
...

After refactoring, one started to override another one (like we normally do for duplicated samples). I modified back to updating dicts instead of overwriting for this sort of reports.

multiqc_picard_ExtractIlluminaBarcodes.yaml - loses LANE

LANE was added there to act as a suffix for the sample name, so I assumed it redundant. Should have kept it for compatibility, and it might be useful for the users downstream. Fixed.

[ewels]

Ok, happy with the two minor bug fixes 👍🏻

Let's merge, then see if anyone complains 😬

Thanks for the major refactoring work here!