Segemehl -Z option doesnt work with split read mapping

Question

[mihir.metkar]

Hi,

My RNA-seq data has transcript fusions and novel splice-junctions. So, I am trying to use segemehl in local alignment (-S) mode to map it. I am using the -Z 20 (--minfraglen <n> min length of a spliced fragment). However, I get fragments of length <20 nt in my output. The command for mapping looks like this-

~/bin/segemehl/segemehl.x --silent -d reference.fa -q reads.fastq -i reference.idx -o file.sam -S -A 95 -Z 20

Example output,

M02096:114:000000000-AEH6V:1:1106:27332:13940 1:N:0:1   16      XIST    78      255     **9M**      *       0     0 **TTTCTCTTG**       IIIFIIIII       NM:i:0  MD:Z:9  NH:i:1  XI:i:0  XL:i:2  XA:Z:Q  XX:i:223        XY:i:231        XQ:i:1  XP:Z:XIST  XU:i:990        XS:i:0

Can someone please explain if I am doing something wrong? Why do I see alignments of <20 nt even after specifying -Z 20?

Answer 1

[David Langenberger]

Dear mihir.metkar,

I talked to one of the developer and he explained it like this:

Even though a min fragment length was specified in some cases the algorithm is shifting smaller parts of a read towards already known junctions! It is a feature!

[mihir.metkar]

Hi David,

Thank you for the answer. Is there a way to turn the feature off?

[Brian Bushnell]

You can postfilter the sam file with Reformat from the BBMap package to get rid of the short fragments:

reformat.sh in=mapped.sam out=filtered.sam minlen=20

[David Langenberger]

It can't be turned off. You can directly contact the developers to change it, or follow Brians idea.

[mihir.metkar]

Ok, thank you all for your answers.