diff --git a/docs/content/tools/alignmentSieve.rst b/docs/content/tools/alignmentSieve.rst
index f10b165fbd..510f8cfc92 100644
--- a/docs/content/tools/alignmentSieve.rst
+++ b/docs/content/tools/alignmentSieve.rst
@@ -84,3 +84,16 @@ As can be seen, such fragments are considered to be on the ``-`` strand, so nega
.. note::
If the ``--shift`` or ``--ATACshift`` options are used, then only properly-paired reads will be used.
+
+Shifted BAM output
+^^^^^^^^^^^^^^^^^^
+
+With ``--shift`` or ``--ATACshift`` and BAM output, each kept read is rewritten as a single match (``M``) covering the shifted read, positioned at the shifted start. SEQ and QUAL are empty and all auxiliary tags are dropped. Because no reference or read bases are needed, expanding shifts and shifts longer than the read do not fail: a read that would shrink below 1 bp is written as a 1 bp record, and is only dropped if it ends up outside the chromosome.
+
+If a downstream step needs a tag such as a single-cell barcode, copy it with ``--keepTags``:
+
+.. code:: bash
+
+ $ alignmentSieve -b possorted.bam --ATACshift --keepTags CB -o shifted.bam
+
+Only the listed tags are written, with their input values; reads that lack a listed tag are written without it. Each tag must be a valid SAM tag name (a letter followed by a letter or digit). Tags are copied verbatim and not checked for meaning: the CIGAR, positions, mate fields, SEQ and QUAL are rewritten, so tags that describe them (for example ``MD``, ``NM``, ``MC``, ``SA``, ``MM`` or ``ML``) will be stale on the output. ``--keepTags`` needs ``--shift`` or ``--ATACshift`` and BAM output. Without it, behaviour is unchanged.
diff --git a/galaxy/wrapper/alignmentSieve.xml b/galaxy/wrapper/alignmentSieve.xml
index ea8a697a7f..83d0e1e8db 100644
--- a/galaxy/wrapper/alignmentSieve.xml
+++ b/galaxy/wrapper/alignmentSieve.xml
@@ -69,6 +69,10 @@
#elif $ATACshift:
--ATACshift
#end if
+ #if ($shift or $ATACshift) and $keepTags and not $BED:
+ #set tags = " ".join(["'{}'".format(x) for x in $keepTags.split(" ")])
+ --keepTags $tags
+ #end if
#if $BED:
--BED
-o '$outFile'
@@ -98,6 +102,10 @@
take strand into consideration for strand-specific
protocols. Note that only properly paired reads are considered."/>
+
+ ^[A-Za-z][A-Za-z0-9]( [A-Za-z][A-Za-z0-9])*$
+
@@ -153,6 +161,21 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
None:
+ """Copy the fixture, adding a CB tag derived from the read name."""
+ with pysam.AlignmentFile(BAMFILE_IN, "rb") as fin, pysam.AlignmentFile(path, "wb", template=fin) as fout:
+ for read in fin:
+ read.set_tag("CB", f"CELL{hash(read.query_name) % 7}-1")
+ fout.write(read)
+
+
+def test_alsieve_shift_drops_tags_by_default():
+ """Shifting without --keepTags drops every tag, SEQ and QUAL (unchanged behaviour)."""
+ _, inbam = tempfile.mkstemp(suffix=".bam")
+ _, outbam = tempfile.mkstemp(suffix=".bam")
+ _bam_with_cell_barcodes(inbam)
+ aln_seive.main(f"--bam {inbam} -o {outbam} --ATACshift".split())
+ with pysam.AlignmentFile(outbam, "rb") as f:
+ reads = list(f)
+ assert reads
+ assert all(not r.get_tags() for r in reads)
+
+
+def test_alsieve_keeptags_cb_matches_default_alignments():
+ """--keepTags CB keeps only CB, with the input value, and does not change positions."""
+ _, inbam = tempfile.mkstemp(suffix=".bam")
+ _, plain = tempfile.mkstemp(suffix=".bam")
+ _, kept = tempfile.mkstemp(suffix=".bam")
+ _bam_with_cell_barcodes(inbam)
+ aln_seive.main(f"--bam {inbam} -o {plain} --ATACshift".split())
+ aln_seive.main(f"--bam {inbam} -o {kept} --ATACshift --keepTags CB".split())
+ with pysam.AlignmentFile(inbam, "rb") as fi:
+ expected_cb = {(r.query_name, r.flag): r.get_tag("CB") for r in fi}
+ with pysam.AlignmentFile(plain, "rb") as fp, pysam.AlignmentFile(kept, "rb") as fk:
+ plain_reads, kept_reads = list(fp), list(fk)
+ assert len(plain_reads) == len(kept_reads) > 0
+ for p, k in zip(plain_reads, kept_reads):
+ assert (p.query_name, p.flag, p.reference_start, p.cigarstring, p.template_length) == (
+ k.query_name, k.flag, k.reference_start, k.cigarstring, k.template_length)
+ assert [t for t, _ in k.get_tags()] == ["CB"]
+ assert k.get_tag("CB") == expected_cb[(k.query_name, k.flag)]
+
+
+def test_alsieve_keeptags_copies_alignment_tags_verbatim():
+ """Tags describing the original alignment are not refused; they are copied as-is (here an integer tag)."""
+ _, inbam = tempfile.mkstemp(suffix=".bam")
+ _, outbam = tempfile.mkstemp(suffix=".bam")
+ with pysam.AlignmentFile(BAMFILE_IN, "rb") as fin, pysam.AlignmentFile(inbam, "wb", template=fin) as fout:
+ for read in fin:
+ read.set_tag("NM", len(read.query_name) % 5)
+ fout.write(read)
+ aln_seive.main(f"--bam {inbam} -o {outbam} --ATACshift --keepTags NM".split())
+ with pysam.AlignmentFile(outbam, "rb") as f:
+ reads = list(f)
+ assert reads
+ for r in reads:
+ assert r.get_tags() == [("NM", len(r.query_name) % 5)]
+
+
+def test_alsieve_keeptags_rejects_malformed_tag_names():
+ """Non-SAM tag names fail as a CLI error, not as a Rust panic (e.g. multi-byte characters)."""
+ _, outbam = tempfile.mkstemp(suffix=".bam")
+ for tag in ("éA", "A-", "1A", "C", "CBX"):
+ try:
+ aln_seive.main(["--bam", BAMFILE_IN, "-o", outbam, "--ATACshift", "--keepTags", tag])
+ except SystemExit as e:
+ assert "SAM tags" in str(e.code)
+ else:
+ raise AssertionError(f"--keepTags {tag!r} was accepted")
+
+
+def test_alsieve_keeptags_deduplicates():
+ """A repeated tag is accepted and written once per record."""
+ _, inbam = tempfile.mkstemp(suffix=".bam")
+ _, outbam = tempfile.mkstemp(suffix=".bam")
+ _bam_with_cell_barcodes(inbam)
+ for tags in ("--keepTags CB CB", "--keepTags CB --keepTags CB"):
+ aln_seive.main(f"--bam {inbam} -o {outbam} --ATACshift {tags}".split())
+ with pysam.AlignmentFile(outbam, "rb") as f:
+ reads = list(f)
+ assert reads
+ assert all([t for t, _ in r.get_tags()] == ["CB"] for r in reads)
+
+def test_alsieve_keeptags_repeated_option_accumulates():
+ """--keepTags CB --keepTags UB keeps both tags, not just the last group."""
+ inbam = tempfile.mkstemp(suffix=".bam")[1]
+ outbam = tempfile.mkstemp(suffix=".bam")[1]
+ with pysam.AlignmentFile(BAMFILE_IN) as src, pysam.AlignmentFile(inbam, "wb", template=src) as dst:
+ for r in src:
+ r.set_tag("CB", "AAAC")
+ r.set_tag("UB", "TTTG")
+ dst.write(r)
+ aln_seive.main(f"--bam {inbam} -o {outbam} --ATACshift --keepTags CB --keepTags UB".split())
+ with pysam.AlignmentFile(outbam) as out:
+ for r in out:
+ assert sorted(t for t, _ in r.get_tags()) == ["CB", "UB"]
+
+def test_alsieve_keeptags_needs_shift_and_bam_output():
+ _, outbam = tempfile.mkstemp(suffix=".bam")
+ for extra in ("", "--ATACshift --BED"):
+ try:
+ aln_seive.main(f"--bam {BAMFILE_IN} -o {outbam} --keepTags CB {extra}".split())
+ except SystemExit:
+ continue
+ raise AssertionError(f"--keepTags CB {extra!r} was accepted")
diff --git a/src/alignmentsieve.rs b/src/alignmentsieve.rs
index 70e917077a..6028556e84 100644
--- a/src/alignmentsieve.rs
+++ b/src/alignmentsieve.rs
@@ -29,6 +29,7 @@ pub fn r_alignmentsieve(
maxfraglen: u32, // maximum fragment length.
label: &str, // user defined label
smartlabels: bool, // derive label from filename
+ keep_tags: Vec, // 2-character aux tags to copy onto shifted records.
) -> PyResult<()> {
// Open input BAM once and stream through in order
let mut bam = Reader::from_path(bamifile)
@@ -47,7 +48,14 @@ pub fn r_alignmentsieve(
let write_filters = filtered_out_readsfile != "None";
let readshift: Vec = shift.extract(py).expect("Failed to extract shift");
- // shift is of length 0, 2, or 4.
+ // shift is of length 0 or 4 (the CLI expands 2 values to 4).
+ let keep_tags: Vec<[u8; 2]> = keep_tags
+ .iter()
+ .map(|t| match t.as_bytes() {
+ [a, b] => [*a, *b],
+ _ => panic!("--keepTags expects 2-character tags, got '{}'", t),
+ })
+ .collect();
// Build chrom_names and chrom_sizes from BAM header — indexed by tid for O(1) lookup, zero allocations per record.
let h = bam.header();
@@ -194,7 +202,7 @@ pub fn r_alignmentsieve(
}
let tid = record.tid() as usize;
- if let Some(shifted) = apply_shift(&record, &readshift, chrom_sizes[tid]) {
+ if let Some(shifted) = apply_shift(&record, &readshift, chrom_sizes[tid], &keep_tags) {
if let Some(ref mut w) = obam {
w.write(&shifted).unwrap_or_else(|e| {
panic!(
@@ -489,7 +497,12 @@ fn query_alignment_end(record: &bam::Record) -> i64 {
read_pos
}
-fn apply_shift(record: &bam::Record, shift: &[i32], chrom_len: u64) -> Option {
+fn apply_shift(
+ record: &bam::Record,
+ shift: &[i32],
+ chrom_len: u64,
+ keep_tags: &[[u8; 2]],
+) -> Option {
if !record.is_proper_pair() {
return None;
}
@@ -554,7 +567,9 @@ fn apply_shift(record: &bam::Record, shift: &[i32], chrom_len: u64) -> Option> = newrec
.aux_iter()
.filter_map(|r| r.ok().map(|(t, _)| t.to_vec()))
@@ -562,6 +577,11 @@ fn apply_shift(record: &bam::Record, shift: &[i32], chrom_len: u64) -> Option