CIGAR operations performed according to strand or to read?

Question

[Matteo Schiavinato]

Hi guys,

I know this question has been asked in similar ways other times, although not exactly the same way, hence I can't find an answer.

The ground for my question is this one: Sam And Reverse Complements

Now:

If I have this read:

>test
AATT

And I want to align it against this reference:

>ref
TTACGCGCGCGCG

The only way an aligned would have to perform this alignment is to reverse it and map it from position 2 on (that is: ATT, leaving the first A out).

The resulting CIGAR operations from this alignment would be:

1S3M

or:

3M1S

To better explain my question: would the CIGAR field show the operations performed on the reference, or on the read? I am finding trouble here, please help me out :)

Answer 1

[h.mon]

I don't know as well, but you teased my curiosity so I set to find out:

1) I grabbed the phiX174 sequence.

2) by hand, created 4 sequences with 70bp, modified the last 16bp do random ATCG

3) reverse-complemented the 4 sequences

4) mapped with BWA.

5) samtools view phiX.sam

READ1   0   NC_001422.1 3921    60  54M16S  *   0   0   TCAGGATTGACACCCTCCCAATTGTATGTTTTCATGCCTCCAAATCTTGGAGGCGCAGTCGTAGCTACGG  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA  NM:i:0  MD:Z:54 AS:i:54 XS:i:0
READ2   0   NC_001422.1 3991    60  54M16S  *   0   0   CTTATTACCCTTCTGAATGTCACGCTGATTATTTTGACTTTGAGCGTATCGAGGGCTACGTAGTCTAGCT  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA  NM:i:0  MD:Z:54 AS:i:54 XS:i:0
READ3   0   NC_001422.1 4061    60  54M16S  *   0   0   TGAGGCTTGTGGCATTTCTACTCTTTCTCAATCCCCAATGCTTGGCTTCCATAATCATCGACTACTATCG  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA  NM:i:0  MD:Z:54 AS:i:54 XS:i:0
READ4   0   NC_001422.1 4131    60  54M16S  *   0   0   ATCAAGCTCTTGGAAGAGATTCTGTCTTTTCGTATGCAGGGCGTTGAGTTCGATGCATGCTAGCTGCTAG  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA  NM:i:0  MD:Z:54 AS:i:54 XS:i:0
READ1rc 16  NC_001422.1 3921    60  54M16S  *   0   0   TCAGGATTGACACCCTCCCAATTGTATGTTTTCATGCCTCCAAATCTTGGAGGCGCAGTCGTAGCTACGG  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA  NM:i:0  MD:Z:54 AS:i:54 XS:i:0
READ2rc 16  NC_001422.1 3991    60  54M16S  *   0   0   CTTATTACCCTTCTGAATGTCACGCTGATTATTTTGACTTTGAGCGTATCGAGGGCTACGTAGTCTAGCT  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA  NM:i:0  MD:Z:54 AS:i:54 XS:i:0
READ3rc 16  NC_001422.1 4061    60  54M16S  *   0   0   TGAGGCTTGTGGCATTTCTACTCTTTCTCAATCCCCAATGCTTGGCTTCCATAATCATCGACTACTATCG  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA  NM:i:0  MD:Z:54 AS:i:54 XS:i:0
READ4rc 16  NC_001422.1 4131    60  54M16S  *   0   0   ATCAAGCTCTTGGAAGAGATTCTGTCTTTTCGTATGCAGGGCGTTGAGTTCGATGCATGCTAGCTGCTAG  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA  NM:i:0  MD:Z:54 AS:i:54 XS:i:0

So bwa mem reverse-complements the reads which map to the reverse strand, and reports the CIGAR of the reverse-complemented read. Is this true for all aligners? I don't know. but I just tested bbmap, and it does:

READ1   0   NC_001422.1 3921    17  54=4X1=6I2=1X2= *   0   0   TCAGGATTGACACCCTCCCAATTGTATGTTTTCATGCCTCCAAATCTTGGAGGCGCAGTCGTAGCTACGG  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA  NM:i:11 AM:i:17
READ2   0   NC_001422.1 3991    17  54=2I2X1=1X2=2X1=2X3=   *   0   0   CTTATTACCCTTCTGAATGTCACGCTGATTATTTTGACTTTGAGCGTATCGAGGGCTACGTAGTCTAGCT  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA  NM:i:9  AM:i:17
READ3   0   NC_001422.1 4061    15  54=1X2=13X  *   0   0   TGAGGCTTGTGGCATTTCTACTCTTTCTCAATCCCCAATGCTTGGCTTCCATAATCATCGACTACTATCG  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA  NM:i:14 AM:i:15
READ4   0   NC_001422.1 4131    16  54=4X1=2X1=2X2=3X1= *   0   0   ATCAAGCTCTTGGAAGAGATTCTGTCTTTTCGTATGCAGGGCGTTGAGTTCGATGCATGCTAGCTGCTAG  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA  NM:i:11 AM:i:16
READ1rc 16  NC_001422.1 3921    17  54=4X1=6I2=1X2= *   0   0   TCAGGATTGACACCCTCCCAATTGTATGTTTTCATGCCTCCAAATCTTGGAGGCGCAGTCGTAGCTACGG  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA  NM:i:11 AM:i:17
READ2rc 16  NC_001422.1 3991    17  54=2I2X1=1X2=2X1=2X3=   *   0   0   CTTATTACCCTTCTGAATGTCACGCTGATTATTTTGACTTTGAGCGTATCGAGGGCTACGTAGTCTAGCT  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA  NM:i:9  AM:i:17
READ3rc 16  NC_001422.1 4061    15  54=1X2=13X  *   0   0   TGAGGCTTGTGGCATTTCTACTCTTTCTCAATCCCCAATGCTTGGCTTCCATAATCATCGACTACTATCG  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA  NM:i:14 AM:i:15
READ4rc 16  NC_001422.1 4131    16  54=4X1=2X1=2X2=3X1= *   0   0   ATCAAGCTCTTGGAAGAGATTCTGTCTTTTCGTATGCAGGGCGTTGAGTTCGATGCATGCTAGCTGCTAG  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA  NM:i:11 AM:i:16

Other aligners may be tested similarly.

[Matteo Schiavinato]

Thanks for the test!

It answers partially to my question. What happens if you repeat the same identical experiment but targeting a portion of the phX174 sequence which is on the reverse strand? Is the CIGAR reported forward or it will be "reversed"?

[h.mon]

Looking (not very deeply) at the SAM specs, I think the CIGAR reports operations based on the reference, from left to right, starting from POS. The strand of the read is encoded on the FLAG, and has no influence on the CIGAR field.

POS: 1-based leftmost mapping POSition of the first CIGAR operation that “consumes” a reference base (see table below). The first base in a reference sequence has coordinate 1. POS is set as 0 for an unmapped read without coordinate. If POS is 0, no assumptions can be made about RNAME and CIGAR.

[Matteo Schiavinato]

I made a test myself this morning, so for future readers here's some confirmations:

• the CIGAR reports the operations performed on the reference FORWARD strand also if the read comes from the reverse strand
• the sequence reported in the SEQ field is the one on the FORWARD strand, even if the original read comes from the reverse strand (but the FLAG field has a 16 instead of 0)

Answer 2

[shouldsee.gem]

Just tested HISAT2 reports CIGAR in reference coordinate as well