6
$\begingroup$

I have a (small) BAM file with CIGAR and MD fields.

Question 1: What tools exists in Python and/or R to reconstruct the alignment between the reference and the read in a BAM? Given that this is a very standard analysis, I feel that there should be R packages or Python libraries which have functions which do this task given a BAM input...

My goal was to create a *tsv file with two columns, one with the reconstructed alignment and the other column with the corresponding sequence. This would be useful in my analysis. (I would either use R data.table, or pandas in Python.)

Here is read and reference:

CGGGCCGGTCCCCCCCCGCCGGGTCCGCCCCCGGC
|||||||||||||||| ||||||||||||||||| 
CGGGCCGGTCCCCCCC-GCCGGGTCCGCCCCCGGG

Here is the desired output dataframe (with only the one example row). We could use an R data.table/data.frame, or a pandas DataFrame. Column1 are the alignments (i.e. the reconstruction of the alignment between the reference and the read based on the CIGAR and MD strings), Column 2 are the corresponding sequences. :

ReadAligned    Reference   
CGGGCCGGTCCCCCCC-GCCGGGTCCGCCCCCGGG    CGGGCCGGTCCCCCCCCGCCGGGTCCGCCCCCGGC
...    ....

Question 2: Using the CIGAR/MD tags, one could then create columns with the number of insertions and the number of deletions. This problem feels much more straightforward, e.g. the python library cigar or the cigar parsing functions of pysam

$\endgroup$
5
  • $\begingroup$ Could you clarify your question: You want to store the alignments in a matrix? Or the CIGAR format in a matrix? $\endgroup$
    – llrs
    Commented Mar 1, 2018 at 8:56
  • $\begingroup$ @Llopis A tsv/csv file with two columns: the first column has the alignments, the second column has the sequence. Does this make sense? $\endgroup$ Commented Mar 1, 2018 at 14:28
  • 1
    $\begingroup$ Do you want sam2pairwise to generate the alignments for you, or do you want to parse the cigar strings yourself in R/python? Writing a cigar parser is tedious, but do-able. It might be easier if you just use sam2pairwise to output the alignments for each read, then you parse that into a table. $\endgroup$
    – conchoecia
    Commented Mar 1, 2018 at 14:51
  • 1
    $\begingroup$ Now your question is reduced to "how to count the number of insertions and deletions" in cigar, which is much easier. With perl, it is just an one-liner echo 78M2I10M5D30M|perl -ne 's/(\d+)([ID])/$h{$2}+=$1/eg; print "$h{I}\t$h{D}\n"'. Someone will probably post python answers. $\endgroup$
    – user172818
    Commented Mar 1, 2018 at 19:25
  • 1
    $\begingroup$ @user172818 I've tried to edit this to be more specific. Is this more organized? The second problem is easier; I'm primarily motivated by the first problem, which is related. $\endgroup$ Commented Mar 1, 2018 at 21:56

3 Answers 3

5
$\begingroup$

If all you want are the number of insertions/deletions, you can use this to extract that info:

from cigar import Cigar
c = c=Cigar('20M5I10M3I5M')
num_insertions = len(list(filter(lambda x: x[1] == 'I', c.items())))
insertion_length = sum([x[0] for x in list(filter(lambda x: x[1] == 'I', c.items()))])

Deletions would be the same thing with 'D' in place of 'I'

$\endgroup$
2
  • $\begingroup$ Is there a standard method to grab the string from the BAM within python? It looks like this is parsed automatically with pysam $\endgroup$ Commented Mar 9, 2018 at 17:58
  • $\begingroup$ I guess pysam would be the standard way $\endgroup$ Commented Mar 12, 2018 at 10:24
5
$\begingroup$

I am answering question 1. However, as I am unfamiliar with python, I am providing a javascript solution only to show the basic logic. To try the code at the end, you need to install node.js and run it with node this-script.js.

Given that this is a very standard analysis

I tend to believe this is a very rare analysis. I have only seen it done twice. Converting CIGAR+MD to padded alignment is complicated in coding, as you have to keep sequence, cigar and MD in sync all the time. The code below generates two strings in one go. It may be logically easier to generate the reference and the query strings separately.

In general, MD is ill defined. It is too difficult to work with. If possible, avoid this tag. You can do a lot of things with the combination of CIGAR and NM.

function md2pad(CIGAR, MD, SEQ) {
    var m, re_cigar = /(\d+)([MIDSHNX=])/g, re_MD = /(\d+)|(\^[A-Za-z]+)|([A-Za-z])/g;
    var cigar = [];
    while ((m = re_cigar.exec(CIGAR)) != null) // parse CIGAR into an array
        if (m[2] != 'H') // do nothing for hard clipping, as it has no effect
            cigar.push([parseInt(m[1]), m[2]]);
    var k = 0, sx = '', sy = ''; // k: k-th cigar operator; x: ref; y: query
    // cx/cy: start coordinate of cigar[k] on ref/query; mx/my: start of current MD op
    var cx = 0, cy = 0, mx = 0, my = 0;
    while ((m = re_MD.exec(MD)) != null) {
        if (m[2] != null) { // deletion from the reference
            var len = m[2].length - 1; // length of deletion
            sx += m[2].substr(1), sy += Array(len+1).join("-");
            mx += len, cx += len, ++k; // consume a deletion D from ref
        } else { // copy or mismatch
            var ml = m[1] != null? parseInt(m[1]) : 1; // length of this MD operation
            while (k < cigar.length && cigar[k][1] != 'D') {
                var cl = cigar[k][0], op = cigar[k][1];
                if (op == 'M' || op == 'X' || op == '=') {
                    if (my + ml < cy + cl) { // an MD ends in the middle of cigar M
                        if (ml > 0) {
                            if (m[3] != null) sx += m[3], sy += SEQ[my];
                            else sx += SEQ.substr(my, ml), sy += SEQ.substr(my, ml);
                        }
                        mx += ml, my += ml, ml = 0;
                        break;
                    } else {
                        var dl = cy + cl - my;
                        sx += SEQ.substr(my, dl), sy += SEQ.substr(my, dl);
                        cx += cl, cy += cl, ++k; // consume an M operator
                        mx += dl, my += dl, ml -= dl;
                    }
                } else if (op == 'I') {
                    sy += SEQ.substr(cy, cl), sx += Array(cl+1).join("-");
                    cy += cl, my += cl, ++k; // consume an insertion I to ref
                } else if (op == 'S') {
                    cy += cl, my += cl, ++k; // consume a soft-clipping S
                } else throw Error("inconsistent MD"); // not working with N
            }
            if (ml != 0) throw Error("inconsistent MD");
        }
    }
    if (cx != mx || cy != my) throw Error("inconsistent MD");
    return [sx, sy];
}

var cigar = '16025H3M2I55M1D108M306H';
var MD = '7T15C24T9^C0G2T6T16T3A1A7G21T17G6C1T6T5A4';
var seq = 'GATCACAGGCCTATCACCCTATTAATCACTCACGGGAGCTCTCCATGCATCTGGTATTTTTTCGGGGGGGGATGCACGCGATAGCATCGCGGGCCGCTGGAACCGGAGCACCCTATGTCGCAGGATCTGTCTTTGATTCCTACCTCATGCCATTATTAATCGCGCCTA';
var s = md2pad(cigar, MD, seq);
console.log(s[0]);
console.log(s[1]);
$\endgroup$
1
  • $\begingroup$ great reference code with the caveat that MD is indeed really annoying to work with. oftentimes easier to just recreate the mismatches from the sequence itself $\endgroup$
    – Colin D
    Commented Dec 20, 2023 at 5:36
0
$\begingroup$

This program in C++ is able to perform the task https://github.com/mlafave/sam2pairwise

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.