Bam file fails to validate due to non-adjacent paired-end mates.

[WaltL]

Hi,

I have been struggling with this for over a week now, and have just about given up.  First, let me say that I have used a small test dataset (10,000 paired-end read subset) in the following process, and everything has worked as far as generating the isoform.result and gene.result files.


1) I assembled 130xe6 paired-end reads using Trinity to generate a de novo transcriptome comprised of 42,000 contigs. 

2)  Next I used the Trinity downstream tools to calculate expression using bowtie as follows:
/usr/local/trinity/latest/util/alignReads.pl --left input_001.fastq --right input_R2_001.fastq --seqType fq --target Trinity.sirex.fasta --aligner bowtie

I get the following error msg upon completion:

"The two mates of paired-end read HWI-ST1082:82:D0MARACXX:8:1101:11359:185247 are not adjacent!
The input file is not valid!"

Looking in early part of the stderr file, I see:

# reads processed: 130702809
# reads with at least one reported alignment: 104850716 (80.22%)
# reads that failed to align: 25852093 (19.78%)

3)  I then installed rsem and ran rsem-prepare-reference file using the Trinity contig set- this worked fine.

4)  I then ran rsem-sam-validator directly on the bowtie_out.nameSorted.PropMapPairsForRSEM.bam that was generated via Trinity alignment (before I knew that rsem-sam-validator was already being called by the downstream Trinity alignment utility), and I got the same "input file is invalid" msg.

5) I next checked the bowtie_out.nameSorted.PropMapPairsForRSEM.bam file with samtools to identify the problem mate pairs, and found that there were 4 pairs with this identical name grouped together as per the output below:

HWI-ST1082:82:D0MARACXX:8:1101:11359:185247     131     comp4093_c0_seq1        73      255     101M    =
       455     483     CTAGATATATAAAAATTTGAATAAAATTTCATTTCTATAAATAAATTATAAATAAATATACAACATTAACTAACATTTAATTTATAGCTCATTTTTTTTGT   111==,2=A?D42CFA?4CAC<@H9E<<CH9<EDE4+9CF499?DD9C@<CE*0*???<DBBFCH94?B8BF:8<7CD:@>474=D###############   XA:i:0  MD:Z:99T0C0     NM:i:2
HWI-ST1082:82:D0MARACXX:8:1101:11359:185247     147     comp4093_c0_seq1        455     255     101M    =
       73      483     ACAAAAAAAATGAGCTATAAATTAAATGTTAGTTAATGTTGTATATTTATTTATAATTTATTTATAGAAATGAAATTTTATTCAAATTTTTATATATCTAG   ###############D=474>@:DC7<8:FB8B?49HCFBBD<???*0*EC<@C9DD?994FC9+4EDE<9HC<<E9H@<CAC4?AFC24D?A=2,==111   XA:i:0  MD:Z:0G0A99     NM:i:2
HWI-ST1082:82:D0MARACXX:8:1101:11359:185247     67      comp4093_c0_seq1        73      255     101M    =
       455     483     CTAGATATATAAAAATTTGAATAAAATTTCATTTCTATAAATAAATTATAAATAAATATACAACATTAACTAACATTTAATTTATAGCTCATTTTTTTTTC   =??DABDD?HHA>?3BF+AA?BA9CHEF>9FGGGGI+*?CAEGIFGEIIEAF><B<FFG>DD<@FCA>F==@BC3@FFH4@@DED>@A@;ACAEE?>>B@?   XA:i:0  MD:Z:101        NM:i:0
HWI-ST1082:82:D0MARACXX:8:1101:11359:185247     83      comp4093_c0_seq1        455     255     101M    =
       73      483     GAAAAAAAAATGAGCTATAAATTAAATGTTAGTTAATGTTGTATATTTATTTATAATTTATTTATAGAAATGAAATTTTATTCAAATTTTTATATATCTAG   ?@B>>?EEACA;@A@>DED@@4HFF@3CB@==F>ACF@<DD>GFF<B<>FAEIIEGFIGEAC?*+IGGGGF9>FEHC9AB?AA+FB3?>AHH?DDBAD??=   XA:i:0  MD:Z:101        NM:i:0
HWI-ST1082:82:D0MARACXX:8:1101:11359:185247     131     comp4093_c0_seq2        73      255     101M    =
       455     483     CTAGATATATAAAAATTTGAATAAAATTTCATTTCTATAAATAAATTATAAATAAATATACAACATTAACTAACATTTAATTTATAGCTCATTTTTTTTGT   111==,2=A?D42CFA?4CAC<@H9E<<CH9<EDE4+9CF499?DD9C@<CE*0*???<DBBFCH94?B8BF:8<7CD:@>474=D###############   XA:i:0  MD:Z:99T0C0     NM:i:2
HWI-ST1082:82:D0MARACXX:8:1101:11359:185247     147     comp4093_c0_seq2        455     255     101M    =
       73      483     ACAAAAAAAATGAGCTATAAATTAAATGTTAGTTAATGTTGTATATTTATTTATAATTTATTTATAGAAATGAAATTTTATTCAAATTTTTATATATCTAG   ###############D=474>@:DC7<8:FB8B?49HCFBBD<???*0*EC<@C9DD?994FC9+4EDE<9HC<<E9H@<CAC4?AFC24D?A=2,==111   XA:i:0  MD:Z:0G0A99     NM:i:2
HWI-ST1082:82:D0MARACXX:8:1101:11359:185247     67      comp4093_c0_seq2        73      255     101M    =
       455     483     CTAGATATATAAAAATTTGAATAAAATTTCATTTCTATAAATAAATTATAAATAAATATACAACATTAACTAACATTTAATTTATAGCTCATTTTTTTTTC   =??DABDD?HHA>?3BF+AA?BA9CHEF>9FGGGGI+*?CAEGIFGEIIEAF><B<FFG>DD<@FCA>F==@BC3@FFH4@@DED>@A@;ACAEE?>>B@?   XA:i:0  MD:Z:101        NM:i:0
HWI-ST1082:82:D0MARACXX:8:1101:11359:185247     83      comp4093_c0_seq2        455     255     101M    =
       73      483     GAAAAAAAAATGAGCTATAAATTAAATGTTAGTTAATGTTGTATATTTATTTATAATTTATTTATAGAAATGAAATTTTATTCAAATTTTTATATATCTAG   ?@B>>?EEACA;@A@>DED@@4HFF@3CB@==F>ACF@<DD>GFF<B<>FAEIIEGFIGEAC?*+IGGGGF9>FEHC9AB?AA+FB3?>AHH?DDBAD??=   XA:i:0  MD:Z:101        NM:i:0

6)  Next, I tried running the samtools filter in hopes that maybe it would clean up the problem pairs (I also don't know if this is the only problem pair set, or just the first instance that rsem-validator "sees," since it seems to exit very quickly given the size of the input bam file (>15 GB).

The samtools command I used was:
samtools view -b -F 12 filename.bam > new_filename.filtered.bam

Unfortunately, the output bam file was identical in size to the input file, and it too gave the same error msg. when I ran rsem-sam-validator on it.

7)  So, I went ahead and tried to run rsem-calculate-expression as follows:
/usr/local/rsem/latest/rsem-calculate-expression --paired-end --bam bowtie_out.nameSorted.PropMapPairsForRSEM.filtered.bam rsem.reference/trinity.rsem_reference rsem_output

After ~48 hours, the program exited (but I can find no error files).  Only the rsem_output.transcript.bam file was created (way too small), and the .result files were not present.  Output directories are clearly incomplete, e.g.  there is an rsem_output.temp/ dir that was not present in the test set.

So, what now?  My questions are:

A)  Is there any way to clean up the bam file so that it will validate and can be used in the next step (I do not script, yet, and would need something "out of the box."

B)  Since the Trinity aligner utility is leveraging rsem and bowtie, does it make any sense to try and rerun the alignment with rsem directly? I have run the aligner via the Trinity downstream utility twice.  It takes ~ 3 days to run, and I get the same result- the bam output files are identical.

Thanks for any suggestions and/or remedies for this problem.   

[Brian Haas]

Hi,

Are you using the latest version of Trinity's alignReads.pl script?
It should be running the included convert-sam-for-rsem (part of the
included RSEM distro), which should fix this type of problem.

This issue had appeared earlier and was somewhat rare - it appeared
when reads multiply mapped to different locations in the same
transcript, as opposed to just different transcripts.

best,

-brian

--
--
Brian J. Haas
Manager, Genome Annotation and Analysis, Research and Development
The Broad Institute
http://broad.mit.edu/~bhaas

[WaltL]

Brian,

Thanks for your reply.  I am actually using the previous version, r2012-05-18/ Trinity disto.  I will get our sys. admin to load the latest version on our cluster tomorrow.  Do you think this is where the problem lies?  One thing I noticed, both with my test and full data set, was that the /bowtie_out dir. contains a number of sam filenames with a .finished extension (listed below); however, all of them are empty.  Is this a normal output?  As I recall, there is a command to force sam files to be created.  Should I do this?

Do I need to re-run the alignment to generate sam file outputs, or can I  simply run convert-sam-for-rsem on my current, non-validating bam file?

Best,

Walt

-rw-r--r-- 1 wlorenz jddlab           0 Jun 26 20:41 bowtie_out.coordSorted.bam.finished
-rw-r--r-- 1 wlorenz jddlab           0 Jun 26 17:05 bowtie_out.coordSorted.sam.finished
-rw-r--r-- 1 wlorenz jddlab           0 Jun 26 16:30 bowtie_out.coordSorted.spliceAdjust.sam.finished
-rw-r--r-- 1 wlorenz jddlab 16749934134 Jun 26 23:57 bowtie_out.nameSorted.bam
-rw-r--r-- 1 wlorenz jddlab 15859830765 Jun 27 08:29 bowtie_out.nameSorted.PropMapPairsForRSEM.bam
-rw-r--r-- 1 wlorenz jddlab           0 Jun 27 02:35 bowtie_out.nameSorted.PropMapPairsForRSEM.sam.finished
-rw-r--r-- 1 wlorenz jddlab           0 Jun 26 22:16 bowtie_out.nameSorted.sam.finished
-rw-r--r-- 1 wlorenz jddlab           0 Jun 26 15:54 bowtie_out.pre.coordSorted.sam.finished
-rw-r--r-- 1 wlorenz jddlab           0 Jun 26 12:55 combined.nameSorted.sam.finished

 


On Sunday, July 8, 2012 10:39:28 AM UTC-4, Brian Haas wrote:

Hi,

Are you using the latest version of Trinity's alignReads.pl script?
It should be running the included convert-sam-for-rsem (part of the
included RSEM distro), which should fix this type of problem.

This issue had appeared earlier and was somewhat rare - it appeared
when reads multiply mapped to different locations in the same
transcript, as opposed to just different transcripts.

best,

-brian

[b...@cs.wisc.edu]

Hi Walt,

You can run convert-sam-for-rsem on your current BAM file. The purpose of
this script is to re-format SAM/BAM files so that RSEM can recognize them.

Best,
Bo

[WaltL]

Hi Bo,

I ran the convert-sam-for-rsem script yesterday, and I still have a file that will not validate, and I get the identical error msg. as before.  When I first ran it on our cluster, it failed due to insufficient disk write space (input sam file is 75 GB).  So I ran it on my Ubuntu box, and it completed after 4.5 hr; however, the output bam file (15.872 GB) is slightly larger than the input bam (15.859 GB), which seems odd to me, as I presumed that after "filtering" the file would be smaller.

Any suggestions as to what I should try next?

Best,

Walt

[Fanping Kong]

I got exactly the same error using the latest Trinity. I filtered reads in the error message and re-ran rsem-sam-validator. But still got similar error from some other reads. That makes me very frustrated. I would like to hear any suggestions!

At this moment, I would like to try run bowtie2 outside of trinity and select reads that properly aligned as input. Could I ask what's your opinion about this?

In trinity, at very beginning 40 maximum alignment can be reported during separated alignment. And the number reduces to 20 when merging both left and right reads. I cannot fully understand why we allow multiple alignment here. So it is difficult for me to choose the alignment number for bowtie2 to report. Could you please help me if you have any idea? Thanks in advance!

Best,
fanping

[b...@cs.wisc.edu]

Hi Fanping,

> I got exactly the same error using the latest Trinity. I filtered reads in
> the error message and re-ran rsem-sam-validator. But still got similar
> error from some other reads. That makes me very frustrated. I would like
> to
> hear any suggestions!

If you can give me a small SAM/BAM file which can trigger the error, I can
see what happens here.

>
> At this moment, I would like to try run bowtie2 outside of trinity and
> select reads that properly aligned as input. Could I ask what's your
> opinion about this?

Why do you want to run bowtie2? One main difference between bowtie 2 and
bowtie is the support of indel alignments. RSEM does not support indel
alignments currently. In addition, because the assembly is constructed
from the reads you are trying to align, I cannot see any reason for
finding indel alignments.

>
> In trinity, at very beginning 40 maximum alignment can be reported during
> separated alignment. And the number reduces to 20 when merging both left
> and right reads. I cannot fully understand why we allow multiple alignment
> here. So it is difficult for me to choose the alignment number for bowtie2
> to report. Could you please help me if you have any idea? Thanks in
> advance!

We allow multiple alignments for better expression estimation. RSEM has a
better model of how RNA-Seq reads are generated than the aligners. So it's
better to let RSEM decide how to deal with this multi-reads. For more
details, please refer to
http://bioinformatics.oxfordjournals.org/content/26/4/493.abstract

Best,
Bo

[Fanping Kong]

Thanks for the promptly reply!

fanping

On Wednesday, July 11, 2012 1:55:33 PM UTC-5, (unknown) wrote:

Hi Fanping,

> I got exactly the same error using the latest Trinity. I filtered reads in
> the error message and re-ran rsem-sam-validator. But still got similar
> error from some other reads. That makes me very frustrated. I would like
> to
> hear any suggestions!

If you can give me a small SAM/BAM file which can trigger the error, I can
see what happens here.

It is so nice that you offer the help. I will generate the SAM/BAM file that trigger the error and send it to you later.

>
> At this moment, I would like to try run bowtie2 outside of trinity and
> select reads that properly aligned as input. Could I ask what's your
> opinion about this?

Why do you want to run bowtie2? One main difference between bowtie 2 and
bowtie is the support of indel alignments. RSEM does not support indel
alignments currently. In addition, because the assembly is constructed
from the reads you are trying to align, I cannot see any reason for
finding indel alignments.

At least, in trinity, they align forward and reverse reads separately. And then merge all the alignment together to generate: 1) one file with all alignment read both proper paired and single for visualization; 2) another file with only properly paired-end aligned reads for further analysis for RSEM. I feel that the whole process took much longer time than running bowtie or bowtie2. So I am just wondering can we use bowtie2 to do the alignment first and then select out the properly aligned reads via the bitwise flag. I have do so but I think it will be much faster. Happy to hear your opinion.

I prefer bowtie2 because it more sensitive than bowtie rather than indel. Thanks.

>
> In trinity, at very beginning 40 maximum alignment can be reported during
> separated alignment. And the number reduces to 20 when merging both left
> and right reads. I cannot fully understand why we allow multiple alignment
> here. So it is difficult for me to choose the alignment number for bowtie2
> to report. Could you please help me if you have any idea? Thanks in
> advance!

We allow multiple alignments for better expression estimation. RSEM has a
better model of how RNA-Seq reads are generated than the aligners. So it's
better to let RSEM decide how to deal with this multi-reads. For more
details, please refer to
http://bioinformatics.oxfordjournals.org/content/26/4/493.abstract

I will read the paper again. Thanks.

 

[b...@cs.wisc.edu]

Hi Fanping,

I see. As far as you can make sure that there are no indel alignments or
local alignments in your bowtie 2 SAM/BAM file, I think RSEM can deal
with it.

Best,
Bo

[Fanping Kong]

This information are very important for me. Thanks.

[b...@cs.wisc.edu]

Hi all,

I have looked into the small SAM file Fanping sent me and fixed the
problem. You can download the package under:

https://github.com/bli25wisc/RSEM/zipball/commit_convert-sam-for-rsem

If it is still not work, please feel free to email me.

Best,
Bo

[WaltL]

Bo,

This script did not work for my bam file.  I e-mailed you some more details about what I've tried so far, and also a test bam error file that will give the error.  Did you get it?

Best,
Walt

[Stuart G]

Bo,

after experiencing the same errors, the patch worked for me.

thanks for looking into it.

Regards,

Stuart

[senhao]

Hi Bo,

Thank you very much.

I encountered the same error msg using convert-sam-for-rsem in RSEM v1.1.21. But it works using the one you provided here (https://github.com/bli25wisc/RSEM/zipball/commit_convert-sam-for-rsem ). It took about 4.5h to complete this process for my 20000000 (100bp paired-end) reads.

Regards,
senhao

[JenW]

Hi Bo,

Does the latest version of RSEM (v1.2.0 Sept 11) include this patch? I've been having a similar problem when I try to use my bam alignment file generated by Trinity in RSEM using v.1.2.0 (error msg: The two mates of paired-end read D3NH4HQ1:134:C10YPACXX:5:1101:2024:31627 are not adjacent! The input file is not valid!)

The convert-sam-for-RSEM is unable to take my bam file as input.

Thanks

Jennifer

[Bo Li]

Hi Jennifer,

What do you mean by "The convert-sam-for-RSEM is unable to take my bam
file as input."?

Best,
Bo

> <http://broad.mit.edu/%7Ebhaas>
> >> >> >>> >>
> >> >> >>> >
> >> >> >>>
> >> >> >>>
> >> >> >>>
> >> >>
> >> >>
> >>
> >>
> >>
>

[Jennifer Wygoda]

Hi Bo,

I use the command:

convert-sam-for-rsem bowtie_out.nameSorted.PropMapPairsForRSEM.bam Converted_PropMapPairsForRSEM

and this is the message I get back:

/usr/local/bin/sam/samtools view -h -o Converted_PropMapPairsForRSEM.bam.temp/input.sam bowtie_out.nameSorted.PropMapPairsForRSEM.bam
"/usr/local/bin/sam/samtools view -h -o Converted_PropMapPairsForRSEM.bam.temp/input.sam bowtie_out.nameSorted.PropMapPairsForRSEM.bam" failed! Plase check if you provide correct parameters/options for the pipeline!


Thanks again,

Jennifer

[b...@cs.wisc.edu]

Hi Jennifer,

Can you check if samtools is successfully compiled? To check it, just try
sam/samtools under RSEM directory. If it is compiled correctly, you'll see
samtools prompt.

Best,
Bo

[b...@cs.wisc.edu]

Hi Jennifer,

Can you check if samtools is successfully compiled? To check it, just try
sam/samtools under RSEM directory. If it is compiled correctly, you'll see
samtools prompt.

Best,
Bo

[Jennifer Wygoda]

Hi Bo,

If I look in the RSEM directory, and move into the sam directory, I see the following:

-rwxr-xr-x+ 1 wadmin com.apple.local.ard_interact 400K Sep 12 12:40 samtools

Is this what you meant?

Thanks

Jennifer

[b...@cs.wisc.edu]

Hi Jennifer,

Yes. It seems that samtools is there. Can you send me a small portion of
your sam/bam file which can trigger the error?

Thanks,
Bo