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:
- I am not sure how to account for the 1-based SAM : should I add
+1
at everyexon_start_coordinate
? Does not look like, as the STAR output has the exact same cigar string that I calculated on the STAR output. - 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!