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

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

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'

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!

  • 1
    $\begingroup$ I strongly encourage you to not do this with normal unix tools and to instead code something in python (or whatever language you prefer). $\endgroup$ – Devon Ryan Aug 30 '17 at 13:28
  • 2
    $\begingroup$ I would use an aligner that outputs proper SAM, such as star, gmap, spaln or minimap2. $\endgroup$ – user172818 Aug 30 '17 at 16:00
  • $\begingroup$ 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) $\endgroup$ – aechchiki Aug 30 '17 at 16:31
  • 1
    $\begingroup$ @user172818 minimap2 works amazingly good... many thanks for redirecting me to it $\endgroup$ – 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
  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:
    1. Order the reads of the same read ID by the number of hard-clipped bases at the start of the CIGAR string.
    2. Build a new CIGAR string from the ordered CIGAR strings
    3. Merge the meta tags from the ordered reads
    4. Write out a single, merged read
  • $\begingroup$ I did that at the end. At the time, I just wanted to reformat the non-spliced alignments from LAST (for example) into spliced alignments. Why? Benchmark the most tools possible. But it's ok, I gave up at the end and just used tools that already gave that spliced alignment directly. $\endgroup$ – aechchiki Nov 18 '18 at 10:28
  • 2
    $\begingroup$ That's great. However, the point of SE is to also help others that might have similar issues in the future by providing answers to questions so that users can upvote good answers. If you feel that my answer is a good one to the question you posed, then please consider up voting and accepting it. $\endgroup$ – winni2k Nov 18 '18 at 10:53
  • 1
    $\begingroup$ If you feel that this answer does not answer the question you really had, then please consider rewriting your question for clarity. $\endgroup$ – winni2k Nov 18 '18 at 10:55

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.