Tophat2: prepare unmapped.bam file for input into a tophat run on alternative genome

[danielsbrewer]

I have some paired-end Illumina RNAseq data and have run tophat2 on it against the human genome. I would like now like to run tophat2 again to align the unmapped bams on some alternative genomes to check for contamination/infection. To do this I need to convert the unmapped.bam into fastq files.

To do this I do the following:
1) Remove any reads without a matching pair

samtools view -f1 -b unmapped.bam > unmapped_paired.bam

2) Sort the reads according to name

samtools sort -n unmapped_paired.bam unmapped_paired_sort.bam

3) Run tophat's bam2fastx to get fastq

bam2fastx -q -Q -A -P -o test unmapped_paired_sort.bam

Unfortunately this reports an error:

Error: couldn't retrieve both reads for pair HISEQ2500-01:110:H7AGVADXX:1:1101:1336:2967. Perhaps the input file is not sorted by name?

The problem is that the unmapped.bam file does not seem to have any information in the RNEXT column about the read name of the matched pair. Anyway three steps just to convert the data back to fastqs seems over the top.

Does anyone have any idea how to fix this problem, or provide a better way to do it?

Thanks

[danielsbrewer]

On further examination, it appears that the FLAGS in the unmapped.bam are inaccurate and even after filtering out the reads without the unpaired flag, there are still reads that are unpaired. I assume this is because the other read of the pair has been mapped.

[dpryan]

You might want to try the "--no-mixed" option for tophat2 next time.

[danielsbrewer]

Yes that would have done the trick. Still playing around with RNAseq data so I am definitely in the learning phase!

The script in the following looks like it will help:

http://seqanswers.com/forums/showthread.php?t=34520

Just giving it a go now.

[danielsbrewer]

Yes that would have done the trick. Still playing around with RNAseq data so I am definitely in the learning phase!

The script in the following looks like it will help:

http://seqanswers.com/forums/showthread.php?t=34520

Just giving it a go now.

[bpb9]

bam2fastx libz error

I too am trying to make a fast file out of the unmapped reads so that I can run top hat on an alternative genome. I get a different error:

samtools sort -n unmapped.bam unmapped_sort.bam
bam2fastx -q -Q -A -o outfile unmapped_sort.bam.bam

I get this error:
bam2fastx: /lib64/libz.so.1: no version information available (required by bam2fastx)

Anyone come across this error before?

[GenoMax]

One possibility is that you are running older versions of libz/libxml2. Are you able to get the bam2fastx to complete (that "error" is likely a warning) otherwise?

[bpb9]

Warning can be ignored

One possibility is that you are running older versions of libz/libxml2. Are you able to get the bam2fastx to complete (that "error" is likely a warning) otherwise?

Hm…sure enough, despite the warning, there is in fact a fastq file produced anyway.

But when I run the program from the cluster's login node (shame on me, I know) I don't get the error, and I still get the fast file. Could that be due to different versions of the program running on the login vs. compute nodes? Any idea?

[GenoMax]

But when I run the program from the cluster's login node (shame on me, I know) I don't get the error, and I still get the fast file. Could that be due to different versions of the program running on the login vs. compute nodes? Any idea?

That is certainly a possibility. On large clusters sometimes a few stray nodes don't get updated properly/fully. If you know which node gave you the error let the admins know. They should be able to manually update that node.

[offspring]

Just a note on this general topic, the script fix_tophat_unmapped_reads.py in https://github.com/cbrueffer/misc_bioinf/ fixes various issues in unmapped.bam files that prevent them from being used in downstream tools.

[fchatonnet]

It might be a very late answer, but apparently, tophat can even accept bam files as input. I tested it by error and it works perfectly, no differences with an alignment with a fastq file obtained after bam2fatsq transformation...
If anyone can confirm that I'm not doing anything wrong, it would be nice.