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, 13)
| qname | flag | rname | pos | mapq | cigar | rnext | pnext | tlen | seq | qual | end | tags |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| str | u16 | cat | i32 | u8 | str | cat | i32 | i32 | str | str | i32 | struct[8] |
| "HWI-ST486:305:C0RH5ACXX:1:2104… | 99 | "1" | 1150 | 40 | "43S13M6872N45M" | "1" | 8048 | 170 | "CAGACAGGAACTAGCAATGCTTGAAATCAA… | "CCCFFFFFHHHHHJIJJJJIJJJJJJJJJJ… | 8079 | {1,"2G55",1,1,40,0,40,"45"} |
| "HWI-ST486:305:C0RH5ACXX:1:2301… | 163 | "1" | 3654 | 40 | "96M" | "1" | 3658 | 105 | "GAGAAAACAAATACATAATCGGAGAAATAC… | "@@BDFFFFHGHHHJJJJGIIIJGEIGFEHB… | 3749 | {1,"96",1,0,40,0,40,null} |
| "HWI-ST486:305:C0RH5ACXX:1:2301… | 83 | "1" | 3658 | 40 | "101M" | "1" | 3654 | -105 | "AAACAAATACATAATCGGAGAAATACAGAT… | "C?CDDACDDDDDB?5DDDEEEDCDCDCDCA… | 3758 | {1,"101",1,0,40,0,40,null} |
| "HWI-ST486:305:C0RH5ACXX:1:2105… | 99 | "1" | 3666 | 40 | "98M" | "1" | 3743 | 178 | "ACATAATCGGAGAAATACAGATTACAGAGA… | "@C@FFFFFHHDHHJJDGGIJJIFIJJGIIG… | 3763 | {1,"98",1,0,40,0,40,null} |
| "HWI-ST486:305:C0RH5ACXX:1:2304… | 99 | "1" | 3675 | 40 | "91M" | "1" | 3737 | 163 | "GAGAAATACAGATTACAGAGAGCGAGAGAG… | "@@@FBDEDHHFHFGBGBECGGIGEHEEEDE… | 3765 | {1,"91",1,0,40,0,40,null} |
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, 16)
| qname | flag | rname | pos | mapq | cigar | rnext | pnext | tlen | seq | qual | end | tags | 5p_start | strand | total_quality |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| str | u16 | cat | i32 | u8 | str | cat | i32 | i32 | str | str | i32 | struct[8] | i64 | str | i64 |
| "HWI-ST486:305:C0RH5ACXX:1:2104… | 99 | "1" | 1150 | 40 | "43S13M6872N45M" | "1" | 8048 | 170 | "CAGACAGGAACTAGCAATGCTTGAAATCAA… | "CCCFFFFFHHHHHJIJJJJIJJJJJJJJJJ… | 8079 | {1,"2G55",1,1,40,0,40,"45"} | 1107 | "+" | 3909 |
| "HWI-ST486:305:C0RH5ACXX:1:2301… | 163 | "1" | 3654 | 40 | "96M" | "1" | 3658 | 105 | "GAGAAAACAAATACATAATCGGAGAAATAC… | "@@BDFFFFHGHHHJJJJGIIIJGEIGFEHB… | 3749 | {1,"96",1,0,40,0,40,null} | 3654 | "+" | 3385 |
| "HWI-ST486:305:C0RH5ACXX:1:2301… | 83 | "1" | 3658 | 40 | "101M" | "1" | 3654 | -105 | "AAACAAATACATAATCGGAGAAATACAGAT… | "C?CDDACDDDDDB?5DDDEEEDCDCDCDCA… | 3758 | {1,"101",1,0,40,0,40,null} | 3758 | "-" | 3690 |
| "HWI-ST486:305:C0RH5ACXX:1:2105… | 99 | "1" | 3666 | 40 | "98M" | "1" | 3743 | 178 | "ACATAATCGGAGAAATACAGATTACAGAGA… | "@C@FFFFFHHDHHJJDGGIJJIFIJJGIIG… | 3763 | {1,"98",1,0,40,0,40,null} | 3666 | "+" | 3371 |
| "HWI-ST486:305:C0RH5ACXX:1:2304… | 99 | "1" | 3675 | 40 | "91M" | "1" | 3737 | 163 | "GAGAAATACAGATTACAGAGAGCGAGAGAG… | "@@@FBDEDHHFHFGBGBECGGIGEHEEEDE… | 3765 | {1,"91",1,0,40,0,40,null} | 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[15]] |
| "HWI-ST486:305:C0RH5ACXX:1:1202… | ["1", "1"] | [50174, 50297] | ["+", "-"] | [3883, 3810] | [{99,"1",50174,40,"101M","1",50197,124,"CATAAACCAACTAATAGTCCTTTAGCCAGATCAAACAAGTAATAATAACTCATACGACAAAATATAGTCTCCAATATCCAACTACAAAGAAGAAAGACAAG","CCCFFFFFHHHFHJJJJHJIIJJHIJJ>HHIIIIIIJJJJJJJGHIJJJIHIIJJIJJJIIJJIJIEGGHIIECHHHHFFFFFEECEEECDDDDDDDDDDD",50274,{1,"101",1,0,40,0,40,null},50174,"+",3883}, {147,"1",50197,40,"101M","1",50174,-124,"AGCCAGATCAAACAAGTAATAATAACTCATACGACAAAATATAGTCTCCAATATCCAACTACAAAGAAGAAAGACAAGTGACTGCGTTTAGTCAAACAAAC","DCCCEDDDCEEEECCCA;DDFFEHEHEHHHHJJJJJIHHBEJJIGJJJIFDHHDHFBDFDDIJGJJJIIJIIIHIJJGH@JJJJJJIGHHHHHFFFFFCBC",50297,{1,"101",1,0,40,0,40,null},50297,"-",3810}] |
| "HWI-ST486:305:C0RH5ACXX:1:1203… | ["1", "1"] | [50547, 50969] | ["+", "-"] | [3464, 3554] | [{99,"1",50547,40,"86M11S","1",50618,172,"TCTCAGCCATCTTGGCGAATAGCATTGGCCAGTATGACTCAATACTAACACCAGCAGCTTTCACCAAGGTCGCGATTTTGTCAGCCGTGATAGCGAT","@@CFFEFFGHDFHEHGIIEII@HH@GEEGHGGGHEIHHHGCDEGIGECGIHAEHIIJB;FECCHGGECCCHFFA>;<BCDDCDD@CB<<@CDACBBD",50632,{1,"86",1,0,40,0,40,null},50547,"+",3464}, {147,"1",50618,40,"14M251N83M4S","1",50547,-172,"GCGATTTTGTCAGCCGTGATAGCGATACCCTCGTCCTCGAGGATCATAACAGCGTAGCTGCAAGCAAGCTCTCCAACTGTCGACATTTTTTTATCCTAAGA","<?:CBCCCCA@BDB@CCC@??DDDC?8@DDBDB@CAADEFHGHE=CED@GIJJIIHGB@FIHGIGGGD?DCIIIHGFFIGGEC<IGIGFHHFFFFEFF@CB",50965,{1,"97",1,0,40,0,40,"45"},50969,"-",3554}] |
| "HWI-ST486:305:C0RH5ACXX:1:2104… | ["1", "1"] | [65849, 66322] | ["+", "-"] | [2577, 2779] | [{163,"1",65849,40,"16M242N54M101N28M","1",66132,131,"CTTTGTCCGACAAAACCTCAGGGGCAACGTAGTTTGGCGTTCCACAAGCTGTATGAAGCAAACCATCTTCCCGAACTTGTCGTGAGAAGGCGCTTAAC","7;?DDD8++@?<FH=?E;BFDF@1???0:)09BB>4(8<AAHIJGCD927AA=>=?@EDC>;>(6-(5;>;C/;;2,::>><@990(:?B155055AC",66289,{1,"98",1,0,40,0,40,"45"},65849,"+",2577}, {83,"1",66132,40,"29M101N61M","1",65849,-131,"CCACAAGCTGTATGAAGCAAACCATCTTCCCGAACTTGTCGTGAGAAGGCGCTTAACCCAAAATCAGAGACTTTCAAAACCCCATTTGCG","?:29BBBBB@;;:@@@A?<:AA>@5;8<@<?==.CA;EF;DCCBCA@FFEFGFD?6EGEGF>DGGABFC<FA+A;FFC0AF??D@?A@@<",66322,{1,"90",1,0,40,0,40,"45"},66322,"-",2779}] |
| "HWI-ST486:305:C0RH5ACXX:1:1206… | ["1", "1"] | [48341, 48495] | ["+", "-"] | [3562, 3578] | [{163,"1",48341,40,"101M","1",48395,155,"GGCAGTTTATAGAGAGCAGCCATGTTGAGACACTCGAAGAACTGTCCGTTGTTACAAGTTCCATCTCCGAAAAAGGCGACAGTGACATCATCACAATCCTG","@?<DFFADHHAHFIEEBHDGHG>@CHGHGIIHGIIJIGFDFIGCDB?FB?6DHIIEHIFHGGGIJIIGDEHFFF>=?>B6?B:>CCDDDDCDDDACAC9@C",48441,{1,"101",1,0,40,0,40,null},48341,"+",3562}, {83,"1",48395,40,"101M","1",48341,-155,"CAAGTTCCATCTCCGAAAAAGGCGACAGTGACATCATCACAATCCTGTTTCAAGACTTCCCTCCTGTACTTGGAGCTAAAGGCAGCACCAGTGGCGACAGG","CDDCCCA>52B@B=CADDBCDCCFFFCB@AAE@CE?4D@JIIIIIGGEGGFAGGFGHFEBGIHG<IGGFIIJHEEFJJIGHEHGGDAIHDHHFFFDDF@@@",48495,{1,"101",1,0,40,0,40,null},48495,"-",3578}] |
| "HWI-ST486:305:C0RH5ACXX:1:1307… | ["1", "1"] | [31514, 31726] | ["+", "-"] | [3803, 3702] | [{99,"1",31520,40,"6S83M12S","1",31536,123,"CCGGTCCATTGAGTACTGGATAGCTTCAACCGCAGACTCAGATGGCAGAAAATCATTCACTGCAACTTCCTTGTTCTCGTTTTTCTTGTAGTCTTCGAAGA","CCCFFDFFHHGDHGIIJJIGIIJIIIJJJGGGGIJ8GIJIJJJJJJJJGIJJJIIJJBDHEIIIHGHHHHFFFFFFEEDDBBDDDCDDDDDEEEED?ABD@",31602,{1,"83",1,0,40,0,40,null},31514,"+",3803}, {147,"1",31536,40,"67M90N34M","1",31520,-123,"GCTTCAACCGCAGACTCAGATGGCAGAAAATCATTCACTGCAACTTCCTTGTTCTCGTTTTTCTTGTAGTCTTCGAAGAAACGACGGATTTCAGAGAGACG","ACCCDB9DDDDDDDDDCCCCCDDCDCAEDEE@;FFDBEDHCCECDAIHGIHJIIIIJJIIGCJIIIJJHBGEJJIJIIIGGIHAIJJJGGGHHFFFFF@CC",31726,{1,"101",1,0,40,0,40,"45"},31726,"-",3702}] |
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:50174:+__1:50297:-" |
| "1:50547:+__1:50969:-" |
| "1:65849:+__1:66322:-" |
| "1:48341:+__1:48495:-" |
| "1:31514:+__1:31726:-" |
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[15]] | 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,{1,"101",1,0,40,0,40,null},10761,"+",3934}, {83,"1",10837,40,"101M","1",10761,-177,"ACATGTGTAAACTGTGTATATATAGGGTAGGTATATGTGTATATATATAGTAATTGACAAATGATTTAGGTTCTAACATATATTCTAAAAGTACTCATGAG","ADFFFFFFHHEHHJIIJJJIJJJJJJIIJIJGJJJJJJIJJJIJJJIJJIJJJIJIJJJJJJJJJJJJJJJJIJJJJJJJJJJJJIJJHHHHFFFFFFCC@",10937,{1,"101",1,0,40,0,40,null},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,{1,"2G55",1,1,40,0,40,"45"},1107,"+",3909}, {147,"1",8048,40,"101M","1",1150,-170,"ACTCAATATCCCCAAACTTGGAAATTAGTTTGTTGCTTGAGGTCTAAGATACTTCTATATATGGAAAAAGATTTTCAAAGCCAGATATTTCCACAAGTTTG","ADDEFDDDDDDFFFFEHEHHHJIJJJJJJJJIJJJJJJJJJIJJJJJJIHHJJJIJJIJJJJJJJJJJJJJJJJJJJJJJJJIJJJJJHHHHFFFFFFCCC",8148,{1,"101",1,0,40,0,40,null},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,{1,"100",1,0,40,0,40,null},11368,"+",3522}, {83,"1",11567,40,"90M1D11M","1",11369,-301,"AGTTATAGCGAATATTATGGATATAATTAGCTAACATTTGAGGCATGTGAACCTGTTATTTTATGTAAATTATATATATAGTTTATATACAAATTGAAAAG","DDDFFFFEEHHGIIJIHIGEGIJJJIIIJHGIHIEGHEGIGGEJJIHEIJIHEGGGIHIIHIJIIHIJIJJGGIJJHIJGJGIJIIHGAHHHFFFFFF@@C",11668,{1,"90^A11",1,1,40,0,40,null},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,{1,"101",1,0,40,0,40,null},11483,"+",3953}, {147,"1",11508,40,"101M","1",11483,-126,"CACCTAGTGTCAAGTATTAAGATTAACATATAACTTATGAAGATGATGAGTAACCAACGAGTTATAGCGAATATTATGGATATAATTAGCTAACATTTGAG","AADEEEEFFFFFFHHHHHHHJJJJIHJJJJJIGJIIIIJJJJJJIJJJIJIHFJIJIJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJHHHHHFFFFFCCC",11608,{1,"101",1,0,40,0,40,null},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,{1,"18^C82",1,1,40,0,40,null},11582,"+",4000}, {83,"1",11616,40,"101M","1",11583,-135,"AACCTGTTATTTTATGTAAATTATATATATAGTTTATATACAAAATTGAAAAGATGCGAGTTTCAACATGGTGACAAAAGCCTAATGATGATGAACATCAA","DCCEEEFFFFFFHHGHHHHHJJIJJJIJJJIJJJJJJJJIJJJIJJJIJJIJJJJJJJJJJIIIJJJJJJJJJHJIJJJJJJJJJJJJHHHGHFFFFFCCC",11716,{1,"101",1,0,40,0,40,null},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, 15)
| flag | rname | pos | mapq | cigar | rnext | pnext | tlen | seq | qual | end | tags | 5p_start | strand | total_quality |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| u16 | cat | i32 | u8 | str | cat | i32 | i32 | str | str | i32 | struct[8] | i64 | str | i64 |
| 163 | "1" | 10761 | 40 | "101M" | "1" | 10837 | 177 | "GAAAATTATGATCCGTAGAGACAGCATTTA… | "CCCFFFFFHHGHHJJGIJIJJIHHGGIJJJ… | 10861 | {1,"101",1,0,40,0,40,null} | 10761 | "+" | 3934 |
| 83 | "1" | 10837 | 40 | "101M" | "1" | 10761 | -177 | "ACATGTGTAAACTGTGTATATATAGGGTAG… | "ADFFFFFFHHEHHJIIJJJIJJJJJJIIJI… | 10937 | {1,"101",1,0,40,0,40,null} | 10937 | "-" | 4016 |
| 99 | "1" | 1150 | 40 | "43S13M6872N45M" | "1" | 8048 | 170 | "CAGACAGGAACTAGCAATGCTTGAAATCAA… | "CCCFFFFFHHHHHJIJJJJIJJJJJJJJJJ… | 8079 | {1,"2G55",1,1,40,0,40,"45"} | 1107 | "+" | 3909 |
| 147 | "1" | 8048 | 40 | "101M" | "1" | 1150 | -170 | "ACTCAATATCCCCAAACTTGGAAATTAGTT… | "ADDEFDDDDDDFFFFEHEHHHJIJJJJJJJ… | 8148 | {1,"101",1,0,40,0,40,null} | 8148 | "-" | 3977 |
| 163 | "1" | 11369 | 40 | "1S100M" | "1" | 11567 | 301 | "TGGCACAGAGAGTACAATTCATGAAATTTA… | "?B@ADDFFFDDF+A:EGGHEDF@HEHEHIG… | 11468 | {1,"100",1,0,40,0,40,null} | 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', 99981, '+')]
WARNING: read is missing pair: [('1', 99987, '+')]