mergeBamAlignment error: inappropriate call if not paired read
REQUIRED for all errors and issues:
a) GATK version used: 4.2.2.0
b) Exact command used:
/DATA/peeper_lab/gatk/gatk --java-options '-Xmx16G -XX:ParallelGCThreads=4' MergeBamAlignment R=/DATA/peeper_lab/reference_fasta_files/ensembl_release_106/Homo_sapiens.GRCh38.dna.primary_assembly.fa UNMAPPED_BAM=/DATA/a.vliet/BMS/Melanoma_patients/DNA/Set1/bams/reverted/PD41397a_1_reverted.bam ALIGNED_BAM=/DATA/a.vliet/BMS/Melanoma_patients/DNA/Set1/output/mapped/PD41397a_1_hg38_sorted.bam O=/DATA/a.vliet/BMS/Melanoma_patients/DNA/Set1/output/mapped/PD41397a_1_hg38_clean.bam CREATE_INDEX=true ADD_MATE_CIGAR=true CLIP_ADAPTERS=false CLIP_OVERLAPPING_READS=true
c) Entire program log:
Using GATK jar /DATA/peeper_lab/gatk/gatk-package-4.2.2.0-2-g0f89f46-SNAPSHOT-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -Xmx16G
-XX:ParallelGCThreads=4 -jar /DATA/peeper_lab/gatk/gatk-package-4.2.2.0-2-g0f89f46-SNAPSHOT-local.jar MergeBamAlignment R=/DATA/peeper_lab/reference_fasta_files/ensembl_r
elease_106/Homo_sapiens.GRCh38.dna.primary_assembly.fa UNMAPPED_BAM=/DATA/a.vliet/BMS/Melanoma_patients/DNA/Set1/bams/reverted/PD41397a_1_reverted.bam ALIGNED_BAM=/DATA/a.
vliet/BMS/Melanoma_patients/DNA/Set1/output/mapped/PD41397a_1_hg38_sorted.bam O=/DATA/a.vliet/BMS/Melanoma_patients/DNA/Set1/output/mapped/PD41397a_1_hg38_clean.bam CREATE
_INDEX=true ADD_MATE_CIGAR=true CLIP_ADAPTERS=false CLIP_OVERLAPPING_READS=true
INFO 2022-11-23 13:32:44 MergeBamAlignment
********** NOTE: Picard's command line syntax is changing.
**********
********** For more information, please see:
********** https://github.com/broadinstitute/picard/wiki/Command-Line-Syntax-Transition-For-Users-(Pre-Transition)
**********
********** The command line looks like this in the new syntax:
**********
********** MergeBamAlignment -R /DATA/peeper_lab/reference_fasta_files/ensembl_release_106/Homo_sapiens.GRCh38.dna.primary_assembly.fa -UNMAPPED_BAM /DATA/a.vliet/BMS/Melanoma_patients/DNA/Set1/bams/reverted/PD41397a_1_reverted.bam -ALIGNED_BAM /DATA/a.vliet/BMS/Melanoma_patients/DNA/Set1/output/mapped/PD41397a_1_hg38_sorted.bam -O /DATA/a.vliet/BMS/Melanoma_patients/DNA/Set1/output/mapped/PD41397a_1_hg38_clean.bam -CREATE_INDEX true -ADD_MATE_CIGAR true -CLIP_ADAPTERS false -CLIP_OVERLAPPING_READS true
**********
13:32:44.924 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/DATA/peeper_lab/gatk/gatk-package-4.2.2.0-2-g0f89f46-SNAPSHOT-local.jar!/com/intel/gkl/native/libgkl_compression.so
[Wed Nov 23 13:32:44 CET 2022] MergeBamAlignment UNMAPPED_BAM=/DATA/a.vliet/BMS/Melanoma_patients/DNA/Set1/bams/reverted/PD41397a_1_reverted.bam ALIGNED_BAM=[/DATA/a.vliet/BMS/Melanoma_patients/DNA/Set1/output/mapped/PD41397a_1_hg38_sorted.bam] OUTPUT=/DATA/a.vliet/BMS/Melanoma_patients/DNA/Set1/output/mapped/PD41397a_1_hg38_clean.bam CLIP_ADAPTERS=false CLIP_OVERLAPPING_READS=true ADD_MATE_CIGAR=true CREATE_INDEX=true REFERENCE_SEQUENCE=/DATA/peeper_lab/reference_fasta_files/ensembl_release_106/Homo_sapiens.GRCh38.dna.primary_assembly.fa ADD_PG_TAG_TO_READS=true PAIRED_RUN=true IS_BISULFITE_SEQUENCE=false ALIGNED_READS_ONLY=false MAX_INSERTIONS_OR_DELETIONS=1 ATTRIBUTES_TO_REVERSE=[OQ, U2] ATTRIBUTES_TO_REVERSE_COMPLEMENT=[E2, SQ] READ1_TRIM=0 READ2_TRIM=0 ALIGNER_PROPER_PAIR_FLAGS=false SORT_ORDER=coordinate PRIMARY_ALIGNMENT_STRATEGY=BestMapq HARD_CLIP_OVERLAPPING_READS=false INCLUDE_SECONDARY_ALIGNMENTS=true UNMAP_CONTAMINANT_READS=false MIN_UNCLIPPED_BASES=32 MATCHING_DICTIONARY_TAGS=[M5, LN] UNMAPPED_READ_STRATEGY=DO_NOT_CHANGE VERBOSITY=INFO QUIET=false VALIDATION_STRINGENCY=STRICT COMPRESSION_LEVEL=2 MAX_RECORDS_IN_RAM=500000 CREATE_MD5_FILE=false GA4GH_CLIENT_SECRETS=client_secrets.json USE_JDK_DEFLATER=false USE_JDK_INFLATER=false
Nov 23, 2022 1:32:45 PM shaded.cloud_nio.com.google.auth.oauth2.ComputeEngineCredentials runningOnComputeEngine
INFO: Failed to detect whether we are running on Google Compute Engine.
[Wed Nov 23 13:32:45 CET 2022] Executing as a.vliet@pasteur on Linux 5.4.0-131-generic amd64; OpenJDK 64-Bit Server VM 1.8.0_352-8u352-ga-1~20.04-b08; Deflater: Intel; Inflater: Intel; Provider GCS is available; Picard version: 4.2.2.0-2-g0f89f46-SNAPSHOT
INFO 2022-11-23 13:32:45 SamAlignmentMerger Processing SAM file(s): [/DATA/a.vliet/BMS/Melanoma_patients/DNA/Set1/output/mapped/PD41397a_1_hg38_sorted.bam]
WARNING 2022-11-23 13:32:45 SamAlignmentMerger Exception merging bam alignment - attempting to sort aligned reads and try again: Inappropriate call if not paired read
INFO 2022-11-23 13:32:49 SamAlignmentMerger Read 1000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:32:55 SamAlignmentMerger Read 2000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:33:00 SamAlignmentMerger Read 3000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:33:05 SamAlignmentMerger Read 4000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:33:10 SamAlignmentMerger Read 5000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:33:10 SamAlignmentMerger Read 5000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:33:15 SamAlignmentMerger Read 6000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:33:20 SamAlignmentMerger Read 7000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:33:25 SamAlignmentMerger Read 8000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:33:30 SamAlignmentMerger Read 9000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:33:34 SamAlignmentMerger Read 10000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:33:39 SamAlignmentMerger Read 11000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:33:44 SamAlignmentMerger Read 12000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:33:48 SamAlignmentMerger Read 13000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:33:52 SamAlignmentMerger Read 14000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:33:57 SamAlignmentMerger Read 15000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:34:02 SamAlignmentMerger Read 16000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:34:07 SamAlignmentMerger Read 17000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:34:12 SamAlignmentMerger Read 18000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:34:16 SamAlignmentMerger Read 19000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:34:21 SamAlignmentMerger Read 20000000 records from alignment SAM/BAM.
INFO 2022-11-23 13:34:22 SamAlignmentMerger Finished reading 20068062 total records from alignment SAM/BAM.
[Wed Nov 23 13:34:22 CET 2022] picard.sam.MergeBamAlignment done. Elapsed time: 1.63 minutes.
Runtime.totalMemory()=6422003712
To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp
java.lang.IllegalStateException: Inappropriate call if not paired read
at htsjdk.samtools.SAMRecord.requireReadPaired(SAMRecord.java:899)
at htsjdk.samtools.SAMRecord.getProperPairFlag(SAMRecord.java:907)
at picard.sam.AbstractAlignmentMerger.setValuesFromAlignment(AbstractAlignmentMerger.java:973)
at picard.sam.AbstractAlignmentMerger.transferAlignmentInfoToFragment(AbstractAlignmentMerger.java:680)
at picard.sam.AbstractAlignmentMerger.transferAlignmentInfoToPairedRead(AbstractAlignmentMerger.java:752)
at picard.sam.AbstractAlignmentMerger.mergeAlignment(AbstractAlignmentMerger.java:481)
at picard.sam.SamAlignmentMerger.mergeAlignment(SamAlignmentMerger.java:186)
at picard.sam.MergeBamAlignment.doWork(MergeBamAlignment.java:368)
at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:308)
at org.broadinstitute.hellbender.cmdline.PicardCommandLineProgramExecutor.instanceMain(PicardCommandLineProgramExecutor.java:37)
at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:160)
at org.broadinstitute.hellbender.Main.mainEntry(Main.java:203)
at org.broadinstitute.hellbender.Main.main(Main.java:289)
My issue:
Following best practices, I revert my demultiplexed uBAMs with RevertSam, run MarkIlluminaAdapters, do SamToFastq | bwa | queryname sort | MergeBamAlignment (| = piped into the next program).
This has not caused any issues until now-- I noticed a couple of samples hadn't been processed, probably because of snakemake issues. Everything before today produced nice merged bams with zero errors, using the same GATK version.
The portion of my snakemake pipeline that does what i decribed above is unchanged. I cannot figure out where this error is coming from. Both my reverted bam and my aligned bam are queryname sorted. Flagstat output for both:
samtools flagstat <ubam>
20286804 + 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)
20286804 + 0 paired in sequencing
10143402 + 0 read1
10143402 + 0 read2
0 + 0 properly paired (0.00% : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
samtools flagstat <aligned+sorted bam>
20068062 + 0 in total (QC-passed reads + QC-failed reads)
129462 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
20022206 + 0 mapped (99.77% : 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'm thinking maybe there is a problem with bwa-- I started getting an additional error:
java.lang.IllegalStateException: Aligned record iterator (HS31_27052:1:1101:10000:11450) is behind the unmapped reads (HS33_27105:2:1101:10000:33014)
although I can't figure out what would be going wrong for this sample and not the rest. I'm just going to run everything again from the top.
One read goes missing after alignment:
$ samtools view -h <original ubam> | grep "HS33_27105:2:1101:10000:33014"
HS33_27105:2:1101:10000:33014 163 2 182413796 60 75M = 182413847 126 ACTAAATTTTTCAGCTTATTGCTATTTCCTATAAAACATGGTCAAAGAAACGAAGTTAAAGTAGATAGGCAAGGT BBBBBFFFFFBFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF<BBFFFFFFFFFFFFFFFFFFFB AS:i:75 XS:i:21 MQ:i:60 ms:i:2745 mc:i:182413921 MC:Z:75M MD:Z:75 NM:i:0 RG:Z:27105_2#1
HS33_27105:2:1101:10000:33014 83 2 182413847 60 75M = 182413796 -126 GAAGTTAAAGTAGATAGGCAAGGTAAACACAGATTTAACAGTAGTTTAGTGTTACATATACCTGAGCTTATATGC FFFFF<FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFBBBBB AS:i:75 XS:i:0 BC:Z:TGTCTATC-TCTTTCCC QT:Z:BBBBBFFF /<BBBFFF MQ:i:60 ms:i:2729 mc:i:182413796 MC:Z:75M MD:Z:75 NM:i:0 RG:Z:27105_2#1
$ samtools view -h <ubam with adapters marked> | grep "HS33_27105:2:1101:10000:33014"
HS33_27105:2:1101:10000:33014 77 * 0 0 * * 0 0 GCATATAAGCTCAGGTATATGTAACACTAAACTACTGTTAAATCTGTGTTTACCTTGCCTATCTACTTTAACTTBBBBBFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF<FFFFF BC:Z:TGTCTATC-TCTTTCCC RG:Z:27105_2#1 XS:i:0 QT:Z:BBBBBFFF /<BBBFFF mc:i:182413796 ms:i:2729
HS33_27105:2:1101:10000:33014 141 * 0 0 * * 0 0 ACTAAATTTTTCAGCTTATTGCTATTTCCTATAAAACATGGTCAAAGAAACGAAGTTAAAGTAGATAGGCAAGGBBBBBFFFFFBFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF<BBFFFFFFFFFFFFFFFFFFFB RG:Z:27105_2#1 XS:i:21 mc:i:182413921 ms:i:2745
$ samtools view -h <aligned and sorted bam> | grep "HS33_27105:2:1101:10000:33014"
HS33_27105:2:1101:10000:33014 0 2 181549069 60 75M * 0 0 ACTAAATTTTTCAGCTTATTGCTATTTCCTATAAAACATGGTCAAAGAAACGAAGTTAAAGTAGATAGGCAAGGT 8988:979::8<;>=;;:;<>><;;<<>><;<<>>>==<?@<?=>>@=>>=9:<=<=<>=@<<@<<<@@>=>??8 NM:i:0 MD:Z:75 AS:i:75 XS:i:21 -
Hi Alex van Vliet,
Thank you for writing to the GATK forum! I hope that we can help you sort this out.
I reviewed your inquiry with our developers and received some feedback to share with you.
It appears that you are encountering an error with BWA rather than GATK.Could you please clarify which tool you are using to sort? There are differences between different sort tools that may be contributing to your troubles. You could use samtool's sort tool if you haven't already.
We suggest you try different tools for sorting or sort at a different place in your pipeline. You should also take a look at the fastq to see if both reads are there as well. If they both are, try running sort re-alignment separately and try to identify where the read gets dropped.
I hope this helps! Please let me know if this leads you to success. If you have further questions, please do not hesitate to reach out.
Best,
Anthony -
Hi Alex van Vliet,
It has been a while since we've heard from you, so we'll be closing out this ticket in our system. However, please note that if you still require assistance, you need only respond to this thread, and we can swiftly re-open your ticket and pick up where we left off!
Thank you for being a valued contributor to the GATK community!
Best,
Anthony
-
Hi Anthony,
I ended up solving this by figuring out my fastq was not interleaved. I have no idea how this happened for one file, I suspect it had something to do with an abrupt Snakemake stop and intermediate files not being removed.Anyone encountering this issue, I suggest you first run samtools flagstat on your aligned bam to check if the reads are paired or not. If not, there must be an issue with the bwa mem input files. If you are following GATK best practices, you know to use the -p option for interleaved fastq's in bwa mem. Generate your interleaved fastq with SamToFastq and make sure it's actually interleaved (I just used less -S to take a peek at the file and saw that every header line ended in /1, when it should be /1 and /2 alternating).
Please sign in to leave a comment.
4 comments