[offspring]

Unfortunately I can't reproduce this with the datasets I have available. Since the datasets are probably too large to share, could you run some tests and report the results? The following would be interesting for a start:

Take the read recorded in the Picard exception, e.g. HISEQ01:514:H812AADXX:1:2105:9732:29000 in the example you posted earlier.

Look for the read in both accepted_hits.bam and unfixed unmapped.bam as well as fixed unmapped.bam:

samtools view accepted_hits.bam | grep HISEQ01:514:H812AADXX:1:2105:9732:29000

samtools view unmapped.bam | grep HISEQ01:514:H812AADXX:1:2105:9732:29000

samtools view unmapped_fixup.bam | grep HISEQ01:514:H812AADXX:1:2105:9732:29000

Hopefully we can pin down the problem.

[id0]

I downloaded the newest fix_tophat_unmapped_reads.py. I saw the changed order. I ran it on TopHat 2.0.9 data, but I am still getting the "Mapped mate should have mate reference name" errors.

Sorry about the continued complaints.

[offspring]

I am still getting Picard "Mapped mate should have mate reference name" errors with the new fixed unmapped.bam. Perhaps that particular step is not the culprit.

Yes, this error occurred for files generated with TopHat < 2.0.7. Basically I messed up the order in which the processing steps should be performed. The updated version should fix this.

[id0]

I am still getting Picard "Mapped mate should have mate reference name" errors with the new fixed unmapped.bam. Perhaps that particular step is not the culprit.

[id0]

Hmm, the mate unmapped problem looks like the culprit then indeed. You can check your unmapped.bam file with the following command:

samtools view -f 0x8 unmapped.bam

If the file is OK, there should be at least a couple of read pairs (e.g. two reads with the same name) in there. I'll re-check what TopHat does and follow up with a new bug report to the developers if needed.

I disabled this specific piece of fixup code for TopHat > 2.0.8. You can enable it manually by simply changing "False" to "True" on line 43 for now.

EDIT: I can confirm that the bug in TopHat is still there. I've updated the script to make the fixup code unconditional again.

Looks like you already checked yourself, but just to confirm, I am not getting any reads with that flag. The BAM file has over 100,000,000 reads, so that should definitely be plenty.

Thanks for the quick follow-up! I'll try the updated script and see how it goes.

Also, I tried running Picard on the faulty merged file without CREATE_INDEX argument. It still spits out "Mapped mate should have mate reference name" errors, but completes without terminating. Indexing with samtools works. The merged file still has the same number of reads before and after Picard processing, so the erroneous reads are not being thrown away. The old script might still be useable.

[offspring]

That's what it looked like. However, I tried running this earlier today on TopHat 2.0.9 data. The script ran successfully. When I ran Picard on the merged file, I got "Mapped mate should have mate reference name" errors.

Eventually, Picard terminated with this:

Exception in thread "main" net.sf.samtools.SAMException: Exception when processing alignment for BAM index HISEQ01:514:H812AADXX:1:2105:9732:29000 1/2 51b unmapped read.
	at net.sf.samtools.BAMFileWriter.writeAlignment(BAMFileWriter.java:120)
	at net.sf.samtools.SAMFileWriterImpl.addAlignment(SAMFileWriterImpl.java:168)
	at net.sf.picard.sam.AddOrReplaceReadGroups.doWork(AddOrReplaceReadGroups.java:100)
	at net.sf.picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:177)
	at net.sf.picard.cmdline.CommandLineProgram.instanceMainWithExit(CommandLineProgram.java:119)
	at net.sf.picard.sam.AddOrReplaceReadGroups.main(AddOrReplaceReadGroups.java:66)
Caused by: net.sf.samtools.SAMException: Exception creating BAM index for record HISEQ01:514:H812AADXX:1:2105:9732:29000 1/2 51b unmapped read.
	at net.sf.samtools.BAMIndexer.processAlignment(BAMIndexer.java:94)
	at net.sf.samtools.BAMFileWriter.writeAlignment(BAMFileWriter.java:117)
	... 5 more

I went back to this thread and the only problem I could think of was that last step. If that's not the issue, do you know what may be? This is my first attempt to merge and process TopHat BAM files.

Hmm, the mate unmapped problem looks like the culprit then indeed. You can check your unmapped.bam file with the following command:

samtools view -f 0x8 unmapped.bam

If the file is OK, there should be at least a couple of read pairs (e.g. two reads with the same name) in there. I'll re-check what TopHat does and follow up with a new bug report to the developers if needed.

I disabled this specific piece of fixup code for TopHat > 2.0.8. You can enable it manually by simply changing "False" to "True" on line 43 for now.

EDIT: I can confirm that the bug in TopHat is still there. I've updated the script to make the fixup code unconditional again.

[id0]

It's included, lines 85-89 in the current version. However, the bug has been fixed in TopHat 2.0.9.

That's what it looked like. However, I tried running this earlier today on TopHat 2.0.9 data. The script ran successfully. When I ran Picard on the merged file, I got "Mapped mate should have mate reference name" errors.

Eventually, Picard terminated with this:

Exception in thread "main" net.sf.samtools.SAMException: Exception when processing alignment for BAM index HISEQ01:514:H812AADXX:1:2105:9732:29000 1/2 51b unmapped read.
	at net.sf.samtools.BAMFileWriter.writeAlignment(BAMFileWriter.java:120)
	at net.sf.samtools.SAMFileWriterImpl.addAlignment(SAMFileWriterImpl.java:168)
	at net.sf.picard.sam.AddOrReplaceReadGroups.doWork(AddOrReplaceReadGroups.java:100)
	at net.sf.picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:177)
	at net.sf.picard.cmdline.CommandLineProgram.instanceMainWithExit(CommandLineProgram.java:119)
	at net.sf.picard.sam.AddOrReplaceReadGroups.main(AddOrReplaceReadGroups.java:66)
Caused by: net.sf.samtools.SAMException: Exception creating BAM index for record HISEQ01:514:H812AADXX:1:2105:9732:29000 1/2 51b unmapped read.
	at net.sf.samtools.BAMIndexer.processAlignment(BAMIndexer.java:94)
	at net.sf.samtools.BAMFileWriter.writeAlignment(BAMFileWriter.java:117)
	... 5 more

I went back to this thread and the only problem I could think of was that last step. If that's not the issue, do you know what may be? This is my first attempt to merge and process TopHat BAM files.

[offspring]

What is a good way to fix this manually?

How come this step was not included in the linked fix_tophat_unmapped_reads.py script? Is there some caveat?

It's included, lines 85-89 in the current version. However, the bug has been fixed in TopHat 2.0.9.

[id0]

4. For paired-end reads where both ends are unmapped, Tophat does not set the "mate is unmapped" (0x8) SAM flag on either read. This has to be fixed manually.

What is a good way to fix this manually?

How come this step was not included in the linked fix_tophat_unmapped_reads.py script? Is there some caveat?

[offspring]

Hi Chris,
I tried the ReorderSam just now, but the mistake is the same as before.
Vivienne

I can't think of anything else that might cause this at the moment... I have used this workflow on more than 200 tophat runs without issues.

To be clear, I use the following:

1. fix up unmapped.bam
2. merge accepted_hits.bam and unmapped_fixup.bam (samtools merge)
3. reorder merged file (Picard ReorderSam)
4. add RG header to merged file (Picard AddOrReplaceReadGroups)
5. create SAM index for merged file (samtools index)
6. run RNA-SeQC

Steps 1-3 can be varied, depending on the original sorting of the files. I don't think the order you have done things in should make a difference.

Chris

[vivienne_lovely]

Hi Chris,
I tried the ReorderSam just now, but the mistake is the same as before.
Vivienne

Hi Vivienne,

try running ReorderSam on your merged file (same as with accepted_hits.bam before).

Chris

[offspring]

Hi Vivienne,

try running ReorderSam on your merged file (same as with accepted_hits.bam before).

Chris

[vivienne_lovely]

Good to hear! I use "samtools merge" afterwards, indeed. I've never used MergeBamAlignment, so I'm not sure what the error means. It would be interesting to find out though, I'll test it with one of our datafiles in the next couple of days.

If you figure out what the problem is in the meantime, please let me know.

Chris

Hi Chris,
I processed the tophat output with your instructions. I merged the accepted_hits.bam and unmapped_fixup.bam with "samtools merge", then I gave the merged bam file to "AddOrReplaceReadGroups.jar" to replace the header. But I received the error like follow:

INFO	2013-05-06 12:31:53	AddOrReplaceReadGroups	Processed    11,000,000 records.  Elapsed time: 00:03:33s.  Time for last 1,000,000:   19s.  Last read position: chr11:56,940,830
INFO	2013-05-06 12:32:13	AddOrReplaceReadGroups	Processed    12,000,000 records.  Elapsed time: 00:03:53s.  Time for last 1,000,000:   19s.  Last read position: chr11:64,878,216
[Mon May 06 12:32:14 CST 2013] net.sf.picard.sam.AddOrReplaceReadGroups done. Elapsed time: 3.93 minutes.
Runtime.totalMemory()=813891584
To get help, see http://picard.sourceforge.net/index.shtml#GettingHelp
Exception in thread "main" java.lang.NullPointerException
	at net.sf.samtools.SAMRecord.isValid(SAMRecord.java:1646)
	at net.sf.samtools.BAMFileReader$BAMFileIterator.advance(BAMFileReader.java:540)
	at net.sf.samtools.BAMFileReader$BAMFileIterator.next(BAMFileReader.java:522)
	at net.sf.samtools.BAMFileReader$BAMFileIterator.next(BAMFileReader.java:481)
	at net.sf.samtools.SAMFileReader$AssertableIterator.next(SAMFileReader.java:672)
	at net.sf.samtools.SAMFileReader$AssertableIterator.next(SAMFileReader.java:650)
	at net.sf.picard.sam.AddOrReplaceReadGroups.doWork(AddOrReplaceReadGroups.java:98)
	at net.sf.picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:177)
	at net.sf.picard.cmdline.CommandLineProgram.instanceMainWithExit(CommandLineProgram.java:119)
	at net.sf.picard.sam.AddOrReplaceReadGroups.main(AddOrReplaceReadGroups.java:66)

Is this pipeline right? what should I do with this error?
Thank you for your help.
Vivienne

[offspring]

Good to hear! I use "samtools merge" afterwards, indeed. I've never used MergeBamAlignment, so I'm not sure what the error means. It would be interesting to find out though, I'll test it with one of our datafiles in the next couple of days.

If you figure out what the problem is in the meantime, please let me know.

Chris

[zaki]

Thanks Chris!

Your script worked wonderfully!!

I am assuming you then proceed to merge the unmapped_fixed.bam & mapped.bam using samtools merge?

Wondering if you have encountered picard MergeBAMAliginment tool?

Thanks to your lovely script I managed to merge using samtools and proceed using picard tools. However I think MergeBAMAlignment returns an error if I attempt to merge unmapped.bam + accepted_hits.bam... and also returns the same error if using unmapped_fixed.bam + accepted_hits.bam

Exception in thread "main" net.sf.picard.PicardException: Program Record ID already in use in unmapped BAM file.

Would you have any idea why this might be?