[offspring]

I've just release version 0.3 of tophat-recondition (https://github.com/cbrueffer/tophat-recondition), and it may fix the "Mapped mate should have mate reference name" problem id0 was seeing.

The details are in TopHat bug #16, https://github.com/infphilo/tophat/issues/16 and tophat-recondition bug #1 https://github.com/cbrueffer/tophat-...ition/issues/1.

[offspring]

Yes, this is seperate indeed. The point of tophat-recondition is to fix issues that prevent unmapped.bam from being used by downstream tools (i.e., to make it possible to merge accepted_hits.bam and unmapped.bam, and feed the merged file to a QC software).

[bpb9]

You are correct, there is indeed output for unmapped_fixup.bam but not unmapped.bam.

I still think it's really weird that there are no duplicates in the unmapped file, but that must be a separate issue.

[offspring]

Sounds like it might have an effect, I don't have experience with that option myself. To see the effect of the script itself, try running the following command on unmapped.bam and unmapped_fixup.bam:

samtools view -f 0x8 unmapped.bam

On the original file it should not print anything, but on the fixed file you should get output.

[bpb9]

So it got it to run if I stayed in the same directory as tophat-recondition.py and removed -hv command, and pointed it to the tophat output directories. I think the reason why is because on the cluster, HTseq is installed under one version of python, and pysam under another. If I use the version of python in my local directory, it runs. It produced an output file unmapped_fixup.bam. However when I run samtools flagstat on the unmapped.bam and the unmapped_fixup.bam, the flagstat summaries are the same…

samtools flagstat unmapped_fixup.bam
8903498 + 2048970 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
0 + 0 mapped (0.00%:0.00%)
8903498 + 2048970 paired in sequencing
4380074 + 1096160 read1
4523424 + 952810 read2
0 + 0 properly paired (0.00%:0.00%)
0 + 0 with itself and mate mapped
0 + 0 singletons (0.00%:0.00%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)



samtools flagstat unmapped.bam
8903498 + 2048970 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
0 + 0 mapped (0.00%:0.00%)
8903498 + 2048970 paired in sequencing
4380074 + 1096160 read1
4523424 + 952810 read2
0 + 0 properly paired (0.00%:0.00%)
0 + 0 with itself and mate mapped
0 + 0 singletons (0.00%:0.00%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)


somehow I just can't believe that none of the unmapped reads are duplicates…could this be a consequence of me running tophat with the --no-mix option turned on?

[offspring]

Does this program need to be installed in the same directory as the unmapped bams that need to be reconditioned? I can't seem to get the program to work. I have python and pysam installed, but the program doesn't run.

wget https://github.com/cbrueffer/tophat-...ive/master.zip
unzip master
set PATH to /home/user/tophat_recondition/tophat-recondition-master
set PATH to /cluster/software/python2.7.6/bin
#Try tophat-recondition
tophat-recondition.py -hv /home/scratch/user/Flowcell1/Raw_PE/sampleID/sampleBAM /home/scratch/user/Flowcell1/Raw_PE/sampleID/

Try without -hv, those are optional parameters that will terminate the program (i.e., -h print help and exit, -v print version and exit).
The script can be run from whereever you like, as long as python and pysam can be found.

[bpb9]

How to use?

Short update: the script posted earlier in this thread has seen some updates (Python 3 compatibility, better argument handling etc) and is now called tophat_recondition.

It can be found here: https://github.com/cbrueffer/tophat-recondition

Does this program need to be installed in the same directory as the unmapped bams that need to be reconditioned? I can't seem to get the program to work. I have python and pysam installed, but the program doesn't run.

wget https://github.com/cbrueffer/tophat-...ive/master.zip
unzip master
set PATH to /home/user/tophat_recondition/tophat-recondition-master
set PATH to /cluster/software/python2.7.6/bin
#Try tophat-recondition
tophat-recondition.py -hv /home/scratch/user/Flowcell1/Raw_PE/sampleID/sampleBAM /home/scratch/user/Flowcell1/Raw_PE/sampleID/

[offspring]

Short update: the script posted earlier in this thread has seen some updates (Python 3 compatibility, better argument handling etc) and is now called tophat_recondition.

It can be found here: https://github.com/cbrueffer/tophat-recondition

[offspring]

Yes, that's quite possible. We don't use any options like this in our TopHat commandline, as far as I see.

[id0]

The FASTQs were not filtered.

I am thinking this may be due to TopHat parameters like "--prefilter-multihits" to exclude multi-mapped reads. I believe those reads would not end up in either accepted_hits.bam or unmapped.bam.

[offspring]

The problematic reads are unmapped ones in unmapped.bam with the "mate is unmapped" flag not set. This means their mates would be expected to be mapped and in accepted_hits.bam; however they are not in there. This is a problem, as my script needs these mapped reads to update some fields in their unmapped mates. Because these fields cannot be updated, Picard barfs later on. The samtools output you posted earlier indicates this is also the problem in the original files.

Do you filter your fastq files in any way before handing them to TopHat?

[offspring]

You were right, with the mm10 genome I can reproduce the problem. I'll see what I can find.

[id0]

I think we are just using different reference sequences. I was using the iGenomes version of mm10.

I trimmed the original FASTQs to generate those BAMs.

Which Picard tool are you using? I just tried using SortSam, which does not require a reference sequence. That still gives me the same "Mapped mate should have mate reference name" errors for the merged BAM.

java -Xmx16G -jar ${PICARD_ROOT}/SortSam.jar \
VERBOSITY=WARNING QUIET=true VALIDATION_STRINGENCY=LENIENT \
SORT_ORDER=coordinate \
INPUT=all.bam \
OUTPUT=all.sortsam.bam

These are the specific errors (just the first few):

Ignoring SAM validation error: ERROR: Record 113550, Read name HISEQ01:267:H77T6ADXX:1:1101:15356:2352, Mapped mate should have mate reference name
Ignoring SAM validation error: ERROR: Record 114728, Read name HISEQ01:267:H77T6ADXX:1:1101:6588:2773, Mapped mate should have mate reference name
Ignoring SAM validation error: ERROR: Record 115145, Read name HISEQ01:267:H77T6ADXX:1:1101:15768:2826, Mapped mate should have mate reference name
Ignoring SAM validation error: ERROR: Record 115811, Read name HISEQ01:267:H77T6ADXX:1:1101:12628:3135, Mapped mate should have mate reference name
Ignoring SAM validation error: ERROR: Record 116103, Read name HISEQ01:267:H77T6ADXX:1:1101:19872:3119, Mapped mate should have mate reference name

[offspring]

Can you reproduce the problem you described with the downsampled files you provided?

I am getting the following Picard exception with them:

Exception in thread "main" net.sf.picard.PicardException: Discordant contig lengths: read chr1 LN=195471971, ref chr1 LN=249250621

To me this seems to indicate a problem with how those files were derived from the originals.

[id0]

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.

As expected, accepted_hits.bam does not match anything.

Column 5 in unmapped.bam is 255, but is changed to 0 in unmapped_fixup.bam.

I trimmed my starting FASTQs to just 100,000 reads for testing. The resulting BAMs are less than 10 MBs combined. I posted them here if you want to check for yourself:

https://dl.dropboxusercontent.com/u/5959741/TopHat-100K/accepted_hits.bam

https://dl.dropboxusercontent.com/u/5959741/TopHat-100K/unmapped.bam

https://dl.dropboxusercontent.com/u/5959741/TopHat-100K/unmapped_fixup.bam

(Whoever is reading this thread a few weeks from now, these will probably be broken links. Sorry.)