MergeBamAlignment Problem

Hi,

I'm trying to merge two bam files that I got going through the Drop-seq pipeline. I've check both inputs using ValidateSamFiles and have no errors from either. But when I run MergeBamAlignment, I keep getting the same error when it's running:

WARNING 2018-06-20 11:33:53     SamAlignmentMerger      Exception merging bam alignment - attempting to sort aligned reads and try again: Inappropriate call if not paired read

And this one after it finishes reading in all the records from the alignment SAM/BAM:

Exception in thread "main" java.lang.IllegalStateException: Inappropriate call if not paired read
        at htsjdk.samtools.SAMRecord.requireReadPaired(SAMRecord.java:866)
        at htsjdk.samtools.SAMRecord.getProperPairFlag(SAMRecord.java:874)
        at picard.sam.AbstractAlignmentMerger.setValuesFromAlignment(AbstractAlignmentMerger.java:655)
        at picard.sam.AbstractAlignmentMerger.transferAlignmentInfoToFragment(AbstractAlignmentMerger.java:548)
        at picard.sam.AbstractAlignmentMerger.transferAlignmentInfoToPairedRead(AbstractAlignmentMerger.java:578)
        at picard.sam.AbstractAlignmentMerger.mergeAlignment(AbstractAlignmentMerger.java:390)
        at picard.sam.SamAlignmentMerger.mergeAlignment(SamAlignmentMerger.java:157)
        at picard.sam.MergeBamAlignment.doWork(MergeBamAlignment.java:266)
        at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:209)
        at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:95)
        at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:105)

I'm not sure what the issue is. I've been successful using the same code on previous Drop-seq files but this newest batch is causing these issues. As far as I can tell, they are the same as older runs that worked. Any ideas would be greatly appreciated!

Thanks!
pschnepp

Answers

  • Here's my code

    java -jar $PICARD_JARS/picard.jar MergeBamAlignment UNMAPPED_BAM=unmapped_mc_tagged_polyA_filtered_sorted.bam ALIGNED_BAM=bwa_align_correct_sorted_hg19.bam OUTPUT=bwa_aligned_reads_hg19.merged.bam REFERENCE_SEQUENCE=../../BWAIndex/GSM1629193_hg19_ERCC.fasta INCLUDE_SECONDARY_ALIGNMENTS=false PAIRED_RUN=false
    
  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @pschnepp
    Hi pschnepp,

    Have a look at this thread which I hope will help.

    -Sheila

  • I am using picard to do my sorting. That doesn't seem to be the issue. Any other ideas?

    java -jar $PICARD_JARS/picard.jar SortSam INPUT=bwa_aligned_reads_hg19.bam OUTPUT=bwa_align_sorted_hg19.bam SORT_ORDER=queryname
    java -jar $PICARD_JARS/picard.jar SortSam INPUT=unmapped.bam OUTPUT=unmapped_mc_tagged_polyA_filtered_sorted.bam SORT_ORDER=queryname
    
  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭
    edited July 2018

    Hi @pschnepp,

    Exception in thread "main" java.lang.IllegalStateException: Inappropriate call if not paired read
    at htsjdk.samtools.SAMRecord.requireReadPaired(SAMRecord.java:866)

    indicates something is going on with the paired and single-ended reads in your file and the tool's expectation towards them.

    Can you run samtools flagstat or gatk FlagStat like the researcher in
    https://github.com/broadinstitute/pilon/issues/51?

    The error is from code lines here.

    Do your paired reads have the appropriate SAM flag values? You can check using https://broadinstitute.github.io/picard/explain-flags.html.

    Once you get back to us, I will ask a developer to help look into this issue.

    P.S. It may be helpful for us to have a snippet of your data. Please be prepared to provide one.

  • I ran the samtools flagstat. For the unaligned bam file:

    21063186 + 0 in total (QC-passed reads + QC-failed reads)
    0 + 0 secondary
    0 + 0 supplementary
    0 + 0 duplicates
    0 + 0 mapped (0.00% : N/A)
    0 + 0 paired in sequencing
    0 + 0 read1
    0 + 0 read2
    0 + 0 properly paired (N/A : N/A)
    0 + 0 with itself and mate mapped
    0 + 0 singletons (N/A : N/A)
    0 + 0 with mate mapped to a different chr
    0 + 0 with mate mapped to a different chr (mapQ>=5)
    
    

    And the aligned bam file:

    21231057 + 0 in total (QC-passed reads + QC-failed reads)
    0 + 0 secondary
    167871 + 0 supplementary
    0 + 0 duplicates
    7314991 + 0 mapped (34.45% : N/A)
    0 + 0 paired in sequencing
    0 + 0 read1
    0 + 0 read2
    0 + 0 properly paired (N/A : N/A)
    0 + 0 with itself and mate mapped
    0 + 0 singletons (N/A : N/A)
    0 + 0 with mate mapped to a different chr
    0 + 0 with mate mapped to a different chr (mapQ>=5)
    

    I can provide a portion of my sequencing file for your developers to test.

    Thanks!

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @pschnepp
    Hi pschnepp,

    Yes, I think you can go ahead and submit a bug report following instructions here.

    Thanks,
    Sheila

  • I've uploaded it to the ftp server under the zipped file 'pschnepp_mergebamfile_error'. I've included both a small portion of my bam file as well as reference files I use. I also have my code log and outputs when I run it in that zipped file as well. The .o# file is just a text file of the output log.

    Thanks! Hopefully one of the developers can figure this out.

    Issue · Github
    by Sheila

    Issue Number
    3138
    State
    open
    Last Updated
    Assignee
    Array
  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @pschnepp
    Hi pschnepp,

    Thanks. I will have a look and get back to you.

    -Sheila

  • I thought I would post this but I tried using a different aligner (STAR instead of bwa) and it seems to be working now. I don't know why it worked but I don't have the error message anymore using the same samples.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @pschnepp
    Hi pschnepp,

    Interesting. Are you using RNA-seq samples?

    Thanks,
    Sheila

  • yes. They are either mouse immune cells or cancer cells sequencing using the drop-seq method. so single cell RNA-seq

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @pschnepp
    Hi pschnepp,

    I see. Yes, then I guess that explains it. STAR is what we recommend for RNA data.

    -Sheila

Sign In or Register to comment.