Skip to content

Alignment and reference sequence from CIGAR and MD #110

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 17 commits into
base: master
Choose a base branch
from
Open

Alignment and reference sequence from CIGAR and MD #110

Changes from 1 commit
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
07bbab5
Use CIGAR field and MD aux field to reconstruct alignments from BAM r…
ingolia Sep 1, 2018
2704703
rustfmt
ingolia Sep 4, 2018
0e98f8c
Unmangle table
ingolia Sep 5, 2018
93c268a
Moved MD accessor convenience function onto record
ingolia Sep 23, 2018
8e35160
Fix up tests and fmt
ingolia Sep 23, 2018
5d19951
Simplified Cigar/MD reconstruction interface
ingolia Sep 24, 2018
5ece794
Reference sequence reconstruction method on Record
ingolia Sep 24, 2018
1a0ee99
Better MDString newtype using standard traits
ingolia Sep 28, 2018
32086f0
Use Iter over MatchDesc rather than reified vector
ingolia Sep 28, 2018
401c58c
Generic over MatchDesc iterator type
ingolia Sep 28, 2018
c419e6a
Moved to iterator access to Cigar information
ingolia Sep 28, 2018
27939c6
fmt
ingolia Sep 28, 2018
c4510e5
Moved all read sequence lookup functions onto CigarMDPos and removed …
ingolia Sep 29, 2018
a6d8fcf
fmt
ingolia Sep 29, 2018
1d14ba7
Collecting into an Alignment now goes through FromIterator
ingolia Sep 29, 2018
8a5e3ac
Alignment reconstruction on the Record
ingolia Sep 29, 2018
5072011
Merge branch 'master' into master
johanneskoester Nov 12, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Use Iter over MatchDesc rather than reified vector
  • Loading branch information
@ingolia
ingolia committed Sep 28, 2018
commit 32086f0a889ec76c76a28ee03bc64cc56f3a3cff
63 changes: 33 additions & 30 deletions src/bam/md_align.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
use std::fmt;
use std::fmt::Write;
use std::str;
use std::vec::IntoIter;

use bio_types::alignment::{Alignment, AlignmentMode, AlignmentOperation};

Expand Down Expand Up @@ -415,7 +416,8 @@ impl CigarMDPos {
/// the first sequenced base.
#[derive(Debug)]
pub struct CigarMDPosIter {
md_stack: Vec<MatchDesc>,
md_iter: IntoIter<MatchDesc>,
md_curr: Option<MatchDesc>,
cigar_stack: Vec<Cigar>,
ref_pos_curr: u32,
read_pos_curr: u32,
Expand All @@ -428,14 +430,15 @@ impl CigarMDPosIter {
///
/// * `record` is the BAM record whose alignment will be extracted
pub fn new_from_record(record: &bam::record::Record) -> Result<Self, MDAlignError> {
let mut md_stack = MDString::new_from_record(record)?.0;
md_stack.reverse();
let mut md_iter = MDString::new_from_record(record)?.into_iter();
let mut md_curr = md_iter.next();

let CigarString(mut cigar_stack) = (*record.cigar()).clone();
cigar_stack.reverse();

Ok(CigarMDPosIter {
md_stack: md_stack,
md_iter: md_iter,
md_curr: md_curr,
cigar_stack: cigar_stack,
ref_pos_curr: record.pos() as u32,
read_pos_curr: 0,
Expand All @@ -451,10 +454,10 @@ impl CigarMDPosIter {

let res = match self.cigar_stack.last_mut().unwrap() {
Cigar::Match(ref mut ciglen) => {
let pop_md;
let next_md;
let mmm = match self
.md_stack
.last_mut()
.md_curr
.as_mut()
.ok_or_else(|| MDAlignError::MDvsCIGAR)?
{
MatchDesc::Matches(ref mut mdlen) => {
Expand All @@ -465,7 +468,7 @@ impl CigarMDPosIter {
self.read_pos_curr += 1;
self.ref_pos_curr += 1;
*mdlen -= 1;
pop_md = *mdlen == 0;
next_md = *mdlen == 0;
*ciglen -= 1;
pop_cigar = *ciglen == 0;
pos
Expand All @@ -478,7 +481,7 @@ impl CigarMDPosIter {
};
self.read_pos_curr += 1;
self.ref_pos_curr += 1;
pop_md = true;
next_md = true;
*ciglen -= 1;
pop_cigar = *ciglen == 0;
pos
Expand All @@ -487,16 +490,16 @@ impl CigarMDPosIter {
return Err(MDAlignError::MDvsCIGAR);
}
};
if pop_md {
self.md_stack.pop();
if next_md {
self.md_curr = self.md_iter.next();
}
mmm
}
Cigar::Equal(ref mut ciglen) => {
let pop_md;
let next_md;
let m = match self
.md_stack
.last_mut()
.md_curr
.as_mut()
.ok_or_else(|| MDAlignError::MDvsCIGAR)?
{
MatchDesc::Matches(ref mut mdlen) => {
Expand All @@ -508,23 +511,23 @@ impl CigarMDPosIter {
self.ref_pos_curr += 1;

*mdlen -= 1;
pop_md = *mdlen == 0;
next_md = *mdlen == 0;
*ciglen -= 1;
pop_cigar = *ciglen == 0;
pos
}
_ => return Err(MDAlignError::MDvsCIGAR),
};
if pop_md {
self.md_stack.pop();
if next_md {
self.md_curr = self.md_iter.next();
}
m
}
Cigar::Diff(ref mut ciglen) => {
let pop_md;
let next_md;
let mm = match self
.md_stack
.last_mut()
.md_curr
.as_mut()
.ok_or_else(|| MDAlignError::MDvsCIGAR)?
{
MatchDesc::Mismatch(ref_nt) => {
Expand All @@ -535,15 +538,15 @@ impl CigarMDPosIter {
};
self.read_pos_curr += 1;
self.ref_pos_curr += 1;
pop_md = true;
next_md = true;
*ciglen -= 1;
pop_cigar = *ciglen == 0;
pos
}
_ => return Err(MDAlignError::MDvsCIGAR),
};
if pop_md {
self.md_stack.pop();
if next_md {
self.md_curr = self.md_iter.next();
}
mm
}
Expand All @@ -558,13 +561,13 @@ impl CigarMDPosIter {
pos
}
Cigar::Del(ref mut len) => {
let pop_md;
let next_md;

let del = match self.md_stack.last_mut() {
let del = match self.md_curr.as_mut() {
Some(MatchDesc::Deletion(ref mut ref_nts)) => {
if ref_nts.len() > 0 {
let ref_nt = ref_nts.remove(0);
pop_md = ref_nts.is_empty();
next_md = ref_nts.is_empty();
let pos = CigarMDPos::Delete {
ref_nt: ref_nt,
ref_pos: self.ref_pos_curr,
Expand All @@ -580,8 +583,8 @@ impl CigarMDPosIter {
};
*len -= 1;
pop_cigar = *len == 0;
if pop_md {
self.md_stack.pop();
if next_md {
self.md_curr = self.md_iter.next();
}
del
}
Expand Down Expand Up @@ -730,8 +733,8 @@ impl Iterator for CigarMDPosIter {
}

// Consume 0-length match
while self.md_stack.last() == Some(&MatchDesc::Matches(0)) {
self.md_stack.pop();
while self.md_curr == Some(MatchDesc::Matches(0)) {
self.md_curr = self.md_iter.next();
}

if self.cigar_stack.is_empty() {
Expand Down