Mark duplicates example on SAM file

[ddehueck]

• Added an example of finding duplicate read pairs in a SAM file.
• Added an example section to the docs and an example page for this example
• Validated the result with picard (https://broadinstitute.github.io/picard/command-line-overview.html#MarkDuplicates)

[ddehueck]

Not sure this is the best way to bring this sample data in - it's about 11MB and found it here: https://figshare.com/articles/dataset/Example_SAM_file/1460716

Happy to use something different if there are any recommendations.

[nvictus]

This looks good! I might just add it to the oxbow bucket.

[nvictus]

https://oxbow-ngs.s3.us-east-2.amazonaws.com/Col0_C1.100k.sam

[ddehueck]

Used jupytext to build the .md file from this. Not sure if this should be committed or not.

[nvictus]

No need to commit this if the .md can be roundtripped to ipynb correctly with jupytext, which I was able to do with enough metadata in the YAML frontmatter.

[ddehueck]

This could use a correctness check as this has all been pretty new territory for me :)

[nvictus]

Procedurally, this looks correct to me. The only difference I see from the reference implementation is that they include hard clips ("H"), whereas this seems to only account for soft clipping ("S").

[ddehueck]

Left the output as just the total count - should I add any other summary stats? Maybe mark duplicates in the original df's flag col and show the distribution of duplicates?

[nvictus]

Yeah, I think we could "mark" or even remove duplicates by carrying the full alignments through the pipeline and then exploding / unnesting the original fields at the end.

[nvictus]

e.g., assuming we passed the alignments records through, you could re-flatten them like this:

best_pairs_df.select(
    "qname", "alignments"
).explode("alignments").select(
    pl.col("alignments").struct.unnest()
)

[nvictus]

Sorry that this took so long to get to. It's fantastic! I left some comments and suggestions. Thanks again for making the contribution!

[nvictus]

I might throw in all the required fields upfront. Strand could be represented using the conventional +/- notation people are used to.

df = df.with_columns(
    pl.struct(["pos", "cigar", "flag"])
    .map_elements(get_unclipped_5_prime_start_position, return_dtype=pl.Int64)
    .alias("5p_start"),

    pl.when((pl.col("flag") & STRAND_BIT) == 0)
    .then(pl.lit("+"))
    .otherwise(pl.lit("-"))
    .alias("strand"),

    pl.col("qual").map_elements(get_quality_score_sum, return_dtype=pl.Int64)
    .alias("total_quality")
)

[nvictus]

If we want to carry the original alignment records through a single pass:

df.group_by("qname").agg(
    [
        pl.col("rname").alias("rnames"),
        pl.col("5p_start").alias("5p_starts"),
        pl.col("strand").alias("strands"),
        pl.col("total_quality").alias("total_qualities"),
        pl.struct("*").alias("alignments")
    ]
)

[nvictus]

This could be dropped if we calculate total_qualities as previously suggested.

[nvictus]

Thinking about this on a large dataset, sorting everything by total quality could be very costly. Often bams are coordinate sorted already, so maybe we sort first by dedup key, then by total_quality to minimize the shuffle?

[nvictus]

e.g., assuming we passed the alignments records through, you could re-flatten them like this:

best_pairs_df.select(
    "qname", "alignments"
).explode("alignments").select(
    pl.col("alignments").struct.unnest()
)

[nvictus]

This looks good! I might just add it to the oxbow bucket.

[nvictus]

When I try to convert this back into ipynb with jupytext, the code cells end up just being text.

Somehow, it works correctly by just adding this extra line of metadata:

Suggested change

	format_name: myst
	extension: .md
	format_name: myst

[nvictus]

https://oxbow-ngs.s3.us-east-2.amazonaws.com/Col0_C1.100k.sam

[nvictus]

It might be nice to show this chained into a full pipeline on a lazyframe at the end.

ds = ox.from_sam("data/Col0_C1.100k.sam")

ldf = ds.to_polars(lazy=True).with_columns(
    pl.struct(["pos", "cigar", "flag"])
    .map_elements(get_unclipped_5_prime_start_position, return_dtype=pl.Int64)
    .alias("5p_start"),

    pl.when((pl.col("flag") & STRAND_BIT) == 0)
    .then(pl.lit("+"))
    .otherwise(pl.lit("-"))
    .alias("strand"),

    pl.col("qual").map_elements(get_quality_score_sum, return_dtype=pl.Int64)
    .alias("total_quality")
).group_by("qname").agg(
    [
        pl.col("rname").alias("rnames"),
        pl.col("5p_start").alias("5p_starts"),
        pl.col("strand").alias("strands"),
        pl.col("total_quality").alias("total_qualities"),
        pl.struct(ds.schema.names).alias("alignments")
    ]
).with_columns(
    pl.struct(["rnames", "5p_starts", "strands"])
    .map_elements(
        lambda s: build_dedup_key(s["rnames"], s["5p_starts"], s["strands"]),
        return_dtype=pl.String
    )
    .alias("dedup_key"),
).filter(
    pl.col("dedup_key").is_not_null()
).sort(
    ["dedup_key", "total_qualities"], descending=True
).unique(
    subset=["dedup_key"]
).select(
    "qname", "alignments"
).explode(
    "alignments"
).select(
    pl.col("alignments").struct.unnest()
)

ldf

[nvictus]

Procedurally, this looks correct to me. The only difference I see from the reference implementation is that they include hard clips ("H"), whereas this seems to only account for soft clipping ("S").