1
$\begingroup$

I was wondering how could I obtain the end position of an alignment using its start position and its CIGAR String.

If an alignment present some soft clipping, for example:


    20S55M

Then, would the start position indicated in the SAM file be the first soft clipped base, or the first correctly mapped base?

Thanks

$\endgroup$
2
$\begingroup$

The match is starting with the first mapped base ignoring any softclipped parts of the read.

As an example this are two reads from a bam file of mine:

start   cigar   sequence 
11710   17S83M  NNNNNNNNNNNNNNNNNGGCCTGGCCATGTGTATTTTTTTAAATTTCCACTGATGATTTTGCTGCATGGCCGGTGTTGAGAATGACTGCGCAAATTTGCC   
11713   100M    CTGGCCATGTGTATTTTTTTAAATTTCCACTGATGATTTTGCTGCATGGCCGGTGTTGAGAATGACTGCGCAAATTTGCCGGATTTCCTTTGCTGTTCCT
| improve this answer | |
$\endgroup$
2
$\begingroup$

The start position is the first correctly mapped base. One way to do it, is to use GenomicAlignments in R to read in the alignments.

For example the sam file looks like this:

cat test.sam 
@HD VN:1.3  SO:coordinate
@SQ SN:chrI LN:228861
D00689:415:CD4LUANXX:3:1316:7447:98083#CGGCTATGCAGGACGT 0   chrI    769 0   42M8S   *   0   0   CCATAGCCCGGTGGGGCCTGTGCTGCCACAGAGGGTGTGTCTAGATCGGA  B/B//FFFFFFBFFFF//7FF/<F/FFBFFFF//<F///<F/BFB//<F/  NM:i:1  MD:Z:21A20  AS:i:37 XS:i:37
D00689:415:CD4LUANXX:3:2308:9597:30445#CGGCTATGCAGGACGT 0   chrI    3965    0   3S47M   *   0   0   CCACACACCCACACACCACACCCACACCACACCCACACCACACCCACACC  BBBBBFFFFFFFFFFFFFFFFFFFFFFFFFFFFBFFFFFBF<FFFFF<FF  NM:i:0  MD:Z:47 AS:i:47 XS:i:47

We convert to bam:

samtools view -bS test.sam > test.bam

And read using R:

library(GenomicAlignments)
gr = readGAlignments("test.bam")
gr

GAlignments object with 2 alignments and 0 metadata columns:
      seqnames strand       cigar    qwidth     start       end     width
         <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
  [1]     chrI      +       42M8S        50       769       810        42
  [2]     chrI      +       3S47M        50      3965      4011        47
          njunc
      <integer>
  [1]         0
  [2]         0
  -------
  seqinfo: 1 sequence from an unspecified genome

The start and the end can be retrieved using:

start(gr)
[1]  769 3965
end(gr)
[1]  810 4011
| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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