Convert local alignments to spliced alignments in SAM file

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!

• I strongly encourage you to not do this with normal unix tools and to instead code something in python (or whatever language you prefer). – Devon Ryan Aug 30 '17 at 13:28
• I would use an aligner that outputs proper SAM, such as star, gmap, spaln or minimap2. – user172818 Aug 30 '17 at 16:00
• yes, I already run the same alignment with GMAP, STAR. we have suspect that LAST might work quite well, but as for now the output is not directly helpful. I will try minimap2 (nice to know you can use it for RNA also) – aechchiki Aug 30 '17 at 16:31
• @user172818 minimap2 works amazingly good... many thanks for redirecting me to it – aechchiki Aug 31 '17 at 10:03

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