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 :)
•
written
by
Matteo Schiavinato
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.
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"?
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.
I made a test myself this morning, so for future readers here's some confirmations:
Just tested HISAT2 reports CIGAR in reference coordinate as well
Log in to answer this question.
