Alignment deduplication#
import oxbow as ox
import polars as pl
Walkthrough#
To illustrate step by step, let’s grab a small sample SAM file and materialize it in memory as a Polars DataFrame.
url = "https://oxbow-ngs.s3.us-east-2.amazonaws.com/Col0_C1.100k.sam"
# Let's use oxbow to read in the SAM file as a polars dataframe
df = ox.from_sam(url).to_polars()
df.head()
shape: (5, 12)
| qname | flag | rname | pos | mapq | cigar | rnext | pnext | tlen | seq | qual | end |
|---|---|---|---|---|---|---|---|---|---|---|---|
| str | u16 | cat | i32 | u8 | str | cat | i32 | i32 | str | str | i32 |
| "HWI-ST486:305:C0RH5ACXX:1:2104… | 99 | "1" | 1150 | 40 | "43S13M6872N45M" | "1" | 8048 | 170 | "CAGACAGGAACTAGCAATGCTTGAAATCAA… | "CCCFFFFFHHHHHJIJJJJIJJJJJJJJJJ… | 8079 |
| "HWI-ST486:305:C0RH5ACXX:1:2301… | 163 | "1" | 3654 | 40 | "96M" | "1" | 3658 | 105 | "GAGAAAACAAATACATAATCGGAGAAATAC… | "@@BDFFFFHGHHHJJJJGIIIJGEIGFEHB… | 3749 |
| "HWI-ST486:305:C0RH5ACXX:1:2301… | 83 | "1" | 3658 | 40 | "101M" | "1" | 3654 | -105 | "AAACAAATACATAATCGGAGAAATACAGAT… | "C?CDDACDDDDDB?5DDDEEEDCDCDCDCA… | 3758 |
| "HWI-ST486:305:C0RH5ACXX:1:2105… | 99 | "1" | 3666 | 40 | "98M" | "1" | 3743 | 178 | "ACATAATCGGAGAAATACAGATTACAGAGA… | "@C@FFFFFHHDHHJJDGGIJJIFIJJGIIG… | 3763 |
| "HWI-ST486:305:C0RH5ACXX:1:2304… | 99 | "1" | 3675 | 40 | "91M" | "1" | 3737 | 163 | "GAGAAATACAGATTACAGAGAGCGAGAGAG… | "@@@FBDEDHHFHFGBGBECGGIGEHEEEDE… | 3765 |
Helper functions#
# If the bit 0x10 is set, the read is on the reverse strand
STRAND_BIT = 0x10
def parse_cigar(cigar_str: str) -> list[tuple[str, int]]:
"""Parse the CIGAR string into a list of tuples (operation, length)
Example:
>>> parse_cigar("76M")
[('M', 76)]
>>> parse_cigar("10M1I65M")
[('M', 10), ('I', 1), ('M', 65)]
"""
result = []
current_number = ""
for char in cigar_str:
if char.isdigit():
current_number += char
else:
if current_number:
result.append((char, int(current_number)))
current_number = ""
return result
def get_unclipped_5p_start(row) -> int:
"""
Get the unclipped 5′ start position from the CIGAR string.
Accounts for both soft clips (S) and hard clips (H), matching the
reference implementation in htsjdk.
Args:
row: Row from the SAM file
Returns:
int: Unclipped 5′ start position
"""
pos = row["pos"]
cigar = row["cigar"]
flag = row["flag"]
# Parse CIGAR string
cigar_ops = parse_cigar(cigar)
# Check if reverse strand (bit 0x10 set)
is_reverse = flag & STRAND_BIT
if not is_reverse:
# Forward strand: 5′ end = POS - (number of leading soft/hard-clipped bases)
leading_clips = 0
for op, length in cigar_ops:
if op in ("S", "H"):
leading_clips += length
else:
break # Stop at first non-clip operation
return pos - leading_clips
else:
# Reverse strand: 5′ end is at POS + (aligned length) + (trailing soft/hard-clipped bases) - 1
aligned_length = 0
trailing_clips = 0
# Calculate aligned length (M, =, X, D, N operations)
for op, length in cigar_ops:
if op in ["M", "=", "X", "D", "N"]:
aligned_length += length
# Find trailing soft/hard clips (S or H operations at the end)
for op, length in reversed(cigar_ops):
if op in ("S", "H"):
trailing_clips += length
else:
break # Stop at first non-clip operation from the end
return pos + aligned_length + trailing_clips - 1
def get_quality_score_sum(qual_str):
"""Calculate the sum of quality scores from a string of quality scores"""
return sum(ord(c) - 33 for c in qual_str if c != " ")
def build_dedup_key(rnames, positions, strands):
"""Makes a dedup key for a read pair"""
items = sorted(zip(rnames, positions, strands))
if len(items) < 2:
print(f"WARNING: read is missing pair: {items}")
return None
return f"{items[0][0]}:{items[0][1]}:{items[0][2]}__{items[1][0]}:{items[1][1]}:{items[1][2]}"
Compute derived fields#
df = df.with_columns(
pl.struct(["pos", "cigar", "flag"])
.map_elements(get_unclipped_5p_start, 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")
)
df.head()
shape: (5, 15)
| qname | flag | rname | pos | mapq | cigar | rnext | pnext | tlen | seq | qual | end | 5p_start | strand | total_quality |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| str | u16 | cat | i32 | u8 | str | cat | i32 | i32 | str | str | i32 | i64 | str | i64 |
| "HWI-ST486:305:C0RH5ACXX:1:2104… | 99 | "1" | 1150 | 40 | "43S13M6872N45M" | "1" | 8048 | 170 | "CAGACAGGAACTAGCAATGCTTGAAATCAA… | "CCCFFFFFHHHHHJIJJJJIJJJJJJJJJJ… | 8079 | 1107 | "+" | 3909 |
| "HWI-ST486:305:C0RH5ACXX:1:2301… | 163 | "1" | 3654 | 40 | "96M" | "1" | 3658 | 105 | "GAGAAAACAAATACATAATCGGAGAAATAC… | "@@BDFFFFHGHHHJJJJGIIIJGEIGFEHB… | 3749 | 3654 | "+" | 3385 |
| "HWI-ST486:305:C0RH5ACXX:1:2301… | 83 | "1" | 3658 | 40 | "101M" | "1" | 3654 | -105 | "AAACAAATACATAATCGGAGAAATACAGAT… | "C?CDDACDDDDDB?5DDDEEEDCDCDCDCA… | 3758 | 3758 | "-" | 3690 |
| "HWI-ST486:305:C0RH5ACXX:1:2105… | 99 | "1" | 3666 | 40 | "98M" | "1" | 3743 | 178 | "ACATAATCGGAGAAATACAGATTACAGAGA… | "@C@FFFFFHHDHHJJDGGIJJIFIJJGIIG… | 3763 | 3666 | "+" | 3371 |
| "HWI-ST486:305:C0RH5ACXX:1:2304… | 99 | "1" | 3675 | 40 | "91M" | "1" | 3737 | 163 | "GAGAAATACAGATTACAGAGAGCGAGAGAG… | "@@@FBDEDHHFHFGBGBECGGIGEHEEEDE… | 3765 | 3675 | "+" | 2766 |
Group the reads into pairs#
We assume that the qname corresponds to read pairs. We group by qname and carry the original alignment records through along with the fields needed for deduplication.
pairs_df = 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"),
]
)
pairs_df.head()
shape: (5, 6)
| qname | rnames | 5p_starts | strands | total_qualities | alignments |
|---|---|---|---|---|---|
| str | list[cat] | list[i64] | list[str] | list[i64] | list[struct[14]] |
| "HWI-ST486:305:C0RH5ACXX:1:2307… | ["1", "1"] | [48354, 48495] | ["+", "-"] | [3398, 3709] | [{163,"1",48356,40,"2S99M","1",48395,142,"GGGCAGCCATGTTGAGACACTCGAAGAACTGTCCGTTGTTACAAGTTCCATCTCCGAAAAAGGCGACAGTGACATCATCACAATCCTGTTTCAAGACTTCC","?=8=A?;DHBBFDGFA=CEBH@BE;CFDCDHG>C8CFE:8BG<FE0BFDBBD=@GGAGGGIHEEFDD?>CACCDDDDDDDDCDDDDDDCCDDEACCCCDCC",48454,48354,"+",3398}, {83,"1",48395,40,"101M","1",48356,-142,"CAAGTTCCATCTCCGAAAAAGGCGACAGTGACATCATCACAATCCTGTTTCAAGACTTCCCTCCTGTACTTGGAGCTAAAGGCAGCACCAGTGGCGACAGG","A>C?>ADCC=DDDDDDDDBDDCDD@CFFFHEEEFHEHE@CIHCGHBGDIIIIIHFHIGIIGEIGCIEGJJJIIIHHBJIJIIJIGFEJGHHHHDBFFFCC@",48495,48495,"-",3709}] |
| "HWI-ST486:305:C0RH5ACXX:1:2107… | ["1", "1"] | [47548, 47662] | ["+", "-"] | [3622, 3712] | [{163,"1",47550,40,"2S99M","1",47562,115,"CGACGAACTACAAGAAGAACCTTGTGATAGTTTAGTAACACCAAACAAAAAGCAAAAAGTGCAAATGATTCAAACTGATTCTCTGTAGCATTTAACTTAAT","@BBFD:DADHBHH@AHE:CAAFEHIIJIFIFFIIGCFIIICHIJJJGIIJGIFGJGIJB@GHHHFHGHEEFDFFEE@@ACC@C>CA@CACA;CCDDCAAC@",47648,47548,"+",3622}, {83,"1",47562,40,"100M1S","1",47550,-115,"AAGAACCTTGTGATAGTTTAGTAACACCAAACAAAAAGCAAAAAGTGCAAATGATTCAAACTGATTCTCTGTAGCATTTAACTTAATAAGATATATAGAAG","AAC;5?A;AC@FECFFDFHEEFEIIIIIIGGIHGIIIGIIIIIIGIIIIIIIIGG<IHGIGFH@GIIIIHHHFEIIIHEIIGIIIIGIHBHDHDFDDD:1+",47661,47662,"-",3712}] |
| "HWI-ST486:305:C0RH5ACXX:1:1202… | ["1", "1"] | [48968, 49085] | ["+", "-"] | [3265, 3532] | [{163,"1",48968,40,"98M","1",48985,118,"TGACAACTTCCTGGACAGAGACGACGGGAGATCGACGGGTGGCGTTGGAGTGATTGAGTCTGCGAAGGGAGAGGGAACGGGTGGATCCGAGGAAAGAA","CCCFFFFFHHHHHJJJJJJJ@HIJJJIIBGIIIJIIIJJ-=DBD<@BD3?,9>@CDCD>CCCAB;>BDD<B;BBB<BDDDD0<B?:@AA8>B-8C<9?",49065,48968,"+",3265}, {83,"1",48985,40,"101M","1",48968,-118,"GAGACGACGGGAGATCGACGGGTGGCGTTGGAGTGATTGAGTCTGCGAAGGGAGAGGGAACGGGTGGATCCGAGGAAAGAAGAAGGAGGAGCCAATCGGAT","<<?9988BDDCBA<&5<0<DDDB@B@DC>DEEEDCCBDFEDHBBHEGGGJJIJHEJIHF?GJIGGHIHGGEG@IIJIFIIGGIGGGIJHHGGGFFFFFCC@",49085,49085,"-",3532}] |
| "HWI-ST486:305:C0RH5ACXX:1:1303… | ["1", "1"] | [50508, 50625] | ["+", "-"] | [3870, 3728] | [{163,"1",50508,40,"101M","1",50525,118,"CAGCACCAACGTTCATGATGAGATCAGTCACGTTACGTTTCTCAGCCATCTTGGCGAATAGCATTGGCCAGTATGACTCAATACTAACACCAGCAGCTTTC","CCCFFFFFHHHHGJJJJJJJJJIJJJJIIJJJFHJJJHIIJIJIIJJGIJJJJJJJJJJJJJHHHHHHFFFDFEEEEEEEDDEDDDDDDDDDDDCDDADDC",50608,50508,"+",3870}, {83,"1",50525,40,"101M","1",50508,-118,"ATGAGATCAGTCACGTTACGTTTCTCAGCCATCTTGGCGAATAGCATTGGCCAGTATGACTCAATACTAACACCAGCAGCTTTCACCAAGGTCGCGATTTT",":@DEDDDCCDDB<@BB?<BCDCDDBDDDDDEECDFFFGHHHHJJIIIIIIJIHGJJIHFJIGIGHDIGHDD:IHCJJJJJJJJIGJIJHHHHHFFFFFCCC",50625,50625,"-",3728}] |
| "HWI-ST486:305:C0RH5ACXX:1:1203… | ["1", "1"] | [50617, 51174] | ["+", "-"] | [3631, 3642] | [{163,"1",50617,40,"15M251N83M3S","1",50912,145,"CGCGATTTTGTCAGCCGTGATAGCGATACCCTCGTCCTCGAGGATCATAACAGCGTAGCTGCAAGCAAGCTCTCCAACTGTCGACATTTTTTTATCCTAAG","BC@FFFFFHDFFH9EGI@FHIIJFHIJIIJIIJJHJJJIE<FHBHIIIIHIGEEEDEFBEEECEDBB?5?CDCCDACDADCCDBB@BCCDDDDACDDEDDC",50965,50617,"+",3631}, {83,"1",50912,40,"52M162N49M","1",50617,-145,"TCATAACAGCGTAGCTGCAAGCAAGCTCTCCAACTGTCGACATTTTTTTATCCTAAGAGTAAGATCGAAAATGCTTTTCGCAAGGCGGCGAGAGAGACAAA","CCDDC9DBDDDDCDDCDDDCDCCC???DDDDCCB?@ADC>DA=/GHGHHIGHEHJIIIIIJIIHFBE@JJJJJJJIIGGJIJIGGBGIGHHHHFFFFF@CC",51174,51174,"-",3642}] |
Build deduplication keys#
Build a deduplication key for each read pair and filter out unpaired reads.
pairs_df = pairs_df.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())
pairs_df[("dedup_key",)].head()
WARNING: read is missing pair: [('1', 99981, '+')]
WARNING: read is missing pair: [('1', 99987, '+')]
shape: (5, 1)
| dedup_key |
|---|
| str |
| "1:48354:+__1:48495:-" |
| "1:47548:+__1:47662:-" |
| "1:48968:+__1:49085:-" |
| "1:50508:+__1:50625:-" |
| "1:50617:+__1:51174:-" |
Resolve duplicates#
We choose the best read pair across duplicates by the highest total quality score. Sorting by dedup_key first minimizes the shuffle when the data is already coordinate-sorted.
best_pairs_df = pairs_df.sort(
["dedup_key", "total_qualities"], descending=[False, True]
).unique(
subset=["dedup_key"]
)
# Get the total number of duplicates
total_pair_dups = pairs_df.height - best_pairs_df.height
print("Total pair duplicates:", total_pair_dups)
best_pairs_df.head()
Total pair duplicates: 3123
shape: (5, 7)
| qname | rnames | 5p_starts | strands | total_qualities | alignments | dedup_key |
|---|---|---|---|---|---|---|
| str | list[cat] | list[i64] | list[str] | list[i64] | list[struct[14]] | str |
| "HWI-ST486:305:C0RH5ACXX:1:1105… | ["1", "1"] | [10761, 10937] | ["+", "-"] | [3934, 4016] | [{163,"1",10761,40,"101M","1",10837,177,"GAAAATTATGATCCGTAGAGACAGCATTTAAAAGTTCCTTACGTCCACGTAAAATAATATATCAATTTATACATATACATGTGTAAACTGTGTATATATAG","CCCFFFFFHHGHHJJGIJIJJIHHGGIJJJJJJIHHIJJJJJIGIJJJIIGHFHIJJJJJJJJJIJJJIJJJJJJHHHHHHFHFFFFFFFCE>CFEFEEFC",10861,10761,"+",3934}, {83,"1",10837,40,"101M","1",10761,-177,"ACATGTGTAAACTGTGTATATATAGGGTAGGTATATGTGTATATATATAGTAATTGACAAATGATTTAGGTTCTAACATATATTCTAAAAGTACTCATGAG","ADFFFFFFHHEHHJIIJJJIJJJJJJIIJIJGJJJJJJIJJJIJJJIJJIJJJIJIJJJJJJJJJJJJJJJJIJJJJJJJJJJJJIJJHHHHFFFFFFCC@",10937,10937,"-",4016}] | "1:10761:+__1:10937:-" |
| "HWI-ST486:305:C0RH5ACXX:1:2104… | ["1", "1"] | [1107, 8148] | ["+", "-"] | [3909, 3977] | [{99,"1",1150,40,"43S13M6872N45M","1",8048,170,"CAGACAGGAACTAGCAATGCTTGAAATCAAGAACTTGAATTGAAATAGTTTTTTACTGGATCAGAGACTACTCAATATCCCCAAACTTGGAAATTAGTTTG","CCCFFFFFHHHHHJIJJJJIJJJJJJJJJJJIJJJJJJIJJJJJJJJJBGGJJJJJJJJIIJJJJJJJIHHHHHFFFFFFFDEDEDDDDDCDCCDDDCCED",8079,1107,"+",3909}, {147,"1",8048,40,"101M","1",1150,-170,"ACTCAATATCCCCAAACTTGGAAATTAGTTTGTTGCTTGAGGTCTAAGATACTTCTATATATGGAAAAAGATTTTCAAAGCCAGATATTTCCACAAGTTTG","ADDEFDDDDDDFFFFEHEHHHJIJJJJJJJJIJJJJJJJJJIJJJJJJIHHJJJIJJIJJJJJJJJJJJJJJJJJJJJJJJJIJJJJJHHHHFFFFFFCCC",8148,8148,"-",3977}] | "1:1107:+__1:8148:-" |
| "HWI-ST486:305:C0RH5ACXX:1:1206… | ["1", "1"] | [11368, 11668] | ["+", "-"] | [3522, 3900] | [{163,"1",11369,40,"1S100M","1",11567,301,"TGGCACAGAGAGTACAATTCATGAAATTTATAAGCTTTTTTCCCACTCATCAATTATAATCTCAAGTTATAAATATCAAAACTGAAAAAAGAAGAAGATGA","?B@ADDFFFDDF+A:EGGHEDF@HEHEHIGHGGCD?CFGGDAD<FFH@FGCFHEDHG>BGBGBHCHFHGGIIGIECE>;AEH;@BCEFEDBDD?>CCCD@:",11468,11368,"+",3522}, {83,"1",11567,40,"90M1D11M","1",11369,-301,"AGTTATAGCGAATATTATGGATATAATTAGCTAACATTTGAGGCATGTGAACCTGTTATTTTATGTAAATTATATATATAGTTTATATACAAATTGAAAAG","DDDFFFFEEHHGIIJIHIGEGIJJJIIIJHGIHIEGHEGIGGEJJIHEIJIHEGGGIHIIHIJIIHIJIJJGGIJJHIJGJGIJIIHGAHHHFFFFFF@@C",11668,11668,"-",3900}] | "1:11368:+__1:11668:-" |
| "HWI-ST486:305:C0RH5ACXX:1:2302… | ["1", "1"] | [11483, 11608] | ["+", "-"] | [3953, 3986] | [{99,"1",11483,40,"101M","1",11508,126,"ATGAAATGCTGTAGGCACAGAGGTTCACCTAGTGTCAAGTATTAAGATTAACATATAACTTATGAAGATGATGAGTAACCAACGAGTTATAGCGAATATTA","CCCFFFFFHHHHHJJJJJJJIJJGGIJJJIJJFHHIJJJIIJIJJJJIJJJJJJJJJJJJJJJJJJJJJIIJJJIGHIJIJJHGFF@DEEECEDBDDDEEC",11583,11483,"+",3953}, {147,"1",11508,40,"101M","1",11483,-126,"CACCTAGTGTCAAGTATTAAGATTAACATATAACTTATGAAGATGATGAGTAACCAACGAGTTATAGCGAATATTATGGATATAATTAGCTAACATTTGAG","AADEEEEFFFFFFHHHHHHHJJJJIHJJJJJIGJIIIIJJJJJJIJJJIJIHFJIJIJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJHHHHHFFFFFCCC",11608,11608,"-",3986}] | "1:11483:+__1:11608:-" |
| "HWI-ST486:305:C0RH5ACXX:1:2304… | ["1", "1"] | [11582, 11716] | ["+", "-"] | [4000, 4000] | [{163,"1",11583,40,"1S18M1D82M","1",11616,135,"CATGGATATAATTAGCTAAATTTGAGGCATGTGAACCTGTTATTTTATGTAAATTATATATATAGTTTATATACAAAATTGAAAAGATGCGAGTTTCAACA","CCCFFFFFHHHHHJJJJJJJJJJJJJJIJIIHIJJJJJJIJJJJJJJJJIIIJJJJJJJJJJJJJIIJJJJJJJJJIJIJJGJJJHHHHHFFDCCDEEDDC",11683,11582,"+",4000}, {83,"1",11616,40,"101M","1",11583,-135,"AACCTGTTATTTTATGTAAATTATATATATAGTTTATATACAAAATTGAAAAGATGCGAGTTTCAACATGGTGACAAAAGCCTAATGATGATGAACATCAA","DCCEEEFFFFFFHHGHHHHHJJIJJJIJJJIJJJJJJJJIJJJIJJJIJJIJJJJJJJJJJIIIJJJJJJJJJHJIJJJJJJJJJJJJHHHGHFFFFFCCC",11716,11716,"-",4000}] | "1:11582:+__1:11716:-" |
Recover deduplicated alignments#
Explode and unnest the carried alignment records to recover the original fields.
deduped_df = best_pairs_df.select(
"qname", "alignments"
).explode("alignments").select(
pl.col("alignments").struct.unnest()
)
deduped_df.head()
shape: (5, 14)
| flag | rname | pos | mapq | cigar | rnext | pnext | tlen | seq | qual | end | 5p_start | strand | total_quality |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| u16 | cat | i32 | u8 | str | cat | i32 | i32 | str | str | i32 | i64 | str | i64 |
| 163 | "1" | 10761 | 40 | "101M" | "1" | 10837 | 177 | "GAAAATTATGATCCGTAGAGACAGCATTTA… | "CCCFFFFFHHGHHJJGIJIJJIHHGGIJJJ… | 10861 | 10761 | "+" | 3934 |
| 83 | "1" | 10837 | 40 | "101M" | "1" | 10761 | -177 | "ACATGTGTAAACTGTGTATATATAGGGTAG… | "ADFFFFFFHHEHHJIIJJJIJJJJJJIIJI… | 10937 | 10937 | "-" | 4016 |
| 99 | "1" | 1150 | 40 | "43S13M6872N45M" | "1" | 8048 | 170 | "CAGACAGGAACTAGCAATGCTTGAAATCAA… | "CCCFFFFFHHHHHJIJJJJIJJJJJJJJJJ… | 8079 | 1107 | "+" | 3909 |
| 147 | "1" | 8048 | 40 | "101M" | "1" | 1150 | -170 | "ACTCAATATCCCCAAACTTGGAAATTAGTT… | "ADDEFDDDDDDFFFFEHEHHHJIJJJJJJJ… | 8148 | 8148 | "-" | 3977 |
| 163 | "1" | 11369 | 40 | "1S100M" | "1" | 11567 | 301 | "TGGCACAGAGAGTACAATTCATGAAATTTA… | "?B@ADDFFFDDF+A:EGGHEDF@HEHEHIG… | 11468 | 11368 | "+" | 3522 |
Full streaming pipeline#
Here’s the entire deduplication pipeline chained together on a Polars LazyFrame:
ds = ox.from_sam(url)
ldf = ds.to_polars(lazy=True).with_columns(
pl.struct(["pos", "cigar", "flag"])
.map_elements(get_unclipped_5p_start, 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=[False, True]
).unique(
subset=["dedup_key"]
).select(
"qname", "alignments"
).explode(
"alignments"
).select(
pl.col("alignments").struct.unnest()
)
ldf.show_graph()
Let’s execute the query plan in streaming mode, writing the results to a Parquet file:
ldf.sink_parquet("data/Col0_C1.100k.dedup.pq")
WARNING: read is missing pair: [('1', 99987, '+')]
WARNING: read is missing pair: [('1', 99981, '+')]