Convert local alignments to spliced alignments in SAM file

Question

I mapped RNA reads to reference genome, using LAST in split mode, and converted the MAF alignment to SAM with maf-convert.

My problem is that the transcripts are not reported in a spliced manner, meaning that a transcript_ID is reported several times in the alignment SAM file with an identical bitwise flag in $2. From what I understand, this is due to the fact that only the exons are mapped (one exon per row) and reported as local alignments, as the software cannot deal with splicing patterns combination exon-intron (yet), which behaviour is clear from the CIGAR string.

For a more concrete and visual example let's consider the mapping of transcript FBtr0344900 to the reference genome as given by LAST:

$ cat last_aln.sam | grep FBtr0344900
FBtr0344900 0   4   42774   100 384=144H    *   0   0   TGCGACATTGTTCTACGATGACTACAAAAAATGACCAATAACTTCTATAAACCAATACGATATGTCAGGAGTTTCGGTCCCATACGAAGTCGCCGACTTAAGTATTTTATttttattttgatATGTGTTTGCTATTTTACCTTGTCGAATGCTTCCACACGCTATGAGAATACCATCGTGAGCGTAGCTTACTACTAGAATTTTGTTGAAGTTATTGACAAGCGATGTCTCAATATCTTCCGGACAGCCTCCAGCGTGACATTGCGGGGAATCATGTAACGGCCCAGTAACAGCCTCGGCCAGCACTCGAAGGTTTTCGTTAAGTTTAAGTATTTTATTTGTAGCACCCGCAAACAAAACATTGTGCATAAAGTCGAAGCTCAT    *   NM:i:0  AS:i:2304
FBtr0344900 0   4   43231   100 384H144=    *   0   0   CTGGAAGCTGTTGATTGAACTGGTATTGATGGCAAGTTAAACTGGGCGACTATGTCATTTAAGGGAGATAACGCCTGAGCCGGCAGTTCTTCAATGCAGTTAACGCAATAATGCTGAGAACCGAGTATGATAATAATACACAGT    *   NM:i:0  AS:i:864

And here is the mapping of the same transcript FBtr0344900 as given by STAR - software which reports the alignment right the way I need it:

cat star_aln.sam | grep FBtr0344900
FBtr0344900    0    4    42774    255    384M73N144M    *    0    0TGCGACATTGTTCTACGATGACTACAAAAAATGACCAATAACTTCTATAAACCAATACGATATGTCAGGAGTTTCGGTCCCATACGAAGTCGCCGACTTAAGTATTTTATTTTTATTTTGATATGTGTTTGCTATTTTACCTTGTCGAATGCTTCCACACGCTATGAGAATACCATCGTGAGCGTAGCTTACTACTAGAATTTTGTTGAAGTTATTGACAAGCGATGTCTCAATATCTTCCGGACAGCCTCCAGCGTGACATTGCGGGGAATCATGTAACGGCCCAGTAACAGCCTCGGCCAGCACTCGAAGGTTTTCGTTAAGTTTAAGTATTTTATTTGTAGCACCCGCAAACAAAACATTGTGCATAAAGTCGAAGCTCATCTGGAAGCTGTTGATTGAACTGGTATTGATGGCAAGTTAAACTGGGCGACTATGTCATTTAAGGGAGATAACGCCTGAGCCGGCAGTTCTTCAATGCAGTTAACGCAATAATGCTGAGAACCGAGTATGATAATAATACACAGT    *    NH:i:1    HI:i:1    NM:i:0    MD:Z:528

From discussion with the author, it seems that I cannot get what I need straight from the current LAST version, and is not a technical problem. So, I have to modify the output myself. The aim would be to get at least a CIGAR line representing a complete transcript.

My question is, do you know any software to do this? What I would need is a file containing one line per uniquely mapped transcripts, featuring three fields: transcript_ID, start position and CIGAR string.

From my side, I proceeded like this :

1) extract the fields that I need for the interested in from the SAM file:

$ cut -f1,4,6
FBtr0344900  42774  384=144H
FBtr0344900  43231  384H144=

2) split the CIGAR line to remove items I am not interested in - I simplify the command here, just assuming I only have perfect matches (which I am interested in) and hard-clipping (which I am not interested in):

$ cut -f3 | sed 's/H/_ /g'| sed 's/=/= /g' | sed 's/\w*_\s*//' | sed 's/=//g'
384
144

3) paste the modified CIGAR in the original file, with paste, resulting in:

$ paste (1) (2) | cut -f1,2,4
FBtr0344900  42774  384
FBtr0344900  43231  144

4) merge rows which start ba the same transcript_id:

$ awk -F'\t' -v OFS='\t' '{ x=$1 ; $1="" ; a[x]=a[x]$0 } END { for(x in a) print x,a[x] }' | 
FBtr0344900    42774    384    43231   144 

5) compute the new cigar, calculating the intron length as the arithmetic formula intron_length= (next_exon_start_coordinate - exon_length - previous_exon_start_coordinate), in this simple case above: intron_length= 43231-384-42774

$ awk '{ printf("%s", $1) }; { for (i=4; i<=NF; i+=2) { printf("%s%d%s%d", OFS, $(i-1), OFS, $i - $(i-1) - $(i-2)) } }; { printf("%s%d%s%s", OFS, $NF, OFS, RS) }'
FBtr0344900 384 73 144

6) ideally, with some method that I will figure out, I will modify the CIGAR string (adding alternate M,N at the end of each field except first one), this is how the final file should look like:

FBtr0344900 42774  384M73N144M

The issue of my basic approach are that:

1. I am not sure how to account for the 1-based SAM : should I add +1 at every exon_start_coordinate? Does not look like, as the STAR output has the exact same cigar string that I calculated on the STAR output.
2. BIG issue: This will only work for the reads mapped in forward strand: how to make it feasible with reads mapped in the reverse strand? If I keep my current approach, I will have negative intron sizes...

Any suggestion is welcome!

Answer 1

Just use minimap2 in split alignment mode to realign the reads.

If that is not an option, then you could try using pysam to modify the CIGAR strings. I do not recommend this, as there are many opportunities for subtle bugs because the SAM spec is complex. You would need to:

1. Sort the BAM on read ID so that you can efficiently retrieve reads that you want to merge
2. Iterate through the BAM with pysam while loading all reads of the same read ID
3. For each set of reads with the same read ID:
   4. Order the reads of the same read ID by the number of hard-clipped bases at the start of the CIGAR string.
   5. Build a new CIGAR string from the ordered CIGAR strings
   6. Merge the meta tags from the ordered reads
   7. Write out a single, merged read