SeqTrimMap mapping issue

[rndouglas]

I have Illumina smallRNA-RNAseq reads that have already had their adapters removed and have been filtered for tRNA/rRNA.

I'm trying to use SeqTrimMap to align these reads to my reference genome with the following command:

$SeqTrimMap -b -S -i -v 1 -o ~/seqtrimmap_B_output_496_1 -s -u 18 -t 28 ~/496_1_filtered_18_28nt.fa  ~/zma_b.fa

The program runs well and gives me output in SAM format as requested, but the program will only map reads at lengths 20-24 nt. Does anyone have a suggestion that might allow me to map reads at additional lengths?

[mbayer]

Hi rndouglas,

it sounds like SeqTrimMap is suppressing anything that isn't small RNA.

Here is what I would do:

1. Convert the SAM file that SeqTrimMap produced to BAM format.

2. Extract the unmapped reads from this. You can do that with samtools like so:

samtools view -f4 sample.bam > sample.unmapped.sam

3. Extract the unmapped reads from the sam file you just created and convert these into FASTQ format. I don't know if there is a tool for this so you may have to script it yourself.

4. Map those reads with another mapping tool that will map regardless of length, e.g. Bowtie, BWA etc. Convert your output into BAM if it isn't that already.

5. Merge the two BAM files with samtools merge (see http://samtools.sourceforge.net/samtools.shtml) using the -r flag. This will produce a single BAM file which contains both the shorter and the longer reads, but they will be distinguishable because samtools merge adds a read group tag to each read. This will give you a trail of where each read came from.

6. You can visualise this quite nicely in Tablet (http://bioinf.scri.ac.uk/tablet/). It includes a "colour by read group" colouring scheme which colours each read according to which read group it belongs to. If you select a contig and then click on the 4th tab in the left hand panel you will see a list of the read groups. The check boxes allow you to switch colour for each read group on and off individually so you can grey out one lot and make the other one stand out.

Hope this helps.

Micha

[rndouglas]

Thanks for the suggestion!

I also contacted the program author directly and was told that this sort of thing has popped up with Illumina data before.

With Illumina input the program uses the length of the first line in the input FASTA as the maximum length. In my case, the first line always happened to be a 24nt sRNA.

When I moved a 28nt sRNA to the top of the list the program ran as desired.

[toni_marco]

I'm glad you got around it. Let me know if any other issue appears.
Toni