We've moved!
This site is now read-only. You can find our new documentation site and support forum for posting questions here.
Be sure to read our welcome blog!
This site is now read-only. You can find our new documentation site and support forum for posting questions here.
Be sure to read our welcome blog!
MergeBamAlignment – Select primary alignment
micknudsen
DenmarkMember ✭✭
Hi,
In the current best practices workflow gatk4-data-processing, you recommend using uBAMs instead of FASTQ files. Great idea! However, when it comes to merging with the BWA alignment BAM, there is something that puzzles me.
Here is an example of a paired-end read mapped by BWA:
XXXXXXXX:412:YYYYYYYYY:1:11101:10001:10497 83 chr16 1229894 0 149M = 1229833 -210 GGGCCGCGTAGGCGCGGCTCGCCAGGACGGGCAGCGCCAGCAGCAGCAGATTCAGCATCTGGGGAGCAAGGAGGAGCATCGTGGGCCTGGCCGGGCCTCACAGGGCAGGGCTGGGGGCTACAGATTGTGGGGTGAAGAATGGAGCTGAG AAAAA/E<EEAA</A/<EA<<EEEEEEEE/EEEAAEEAEE/EAEAAEEEEEEEEEEEAEEAAEEAEAEAAEEEEEEEEEEEEAAEEEEAE6EAEEEEEEEE/EEEEEEE/EE/AEAAEEEEEEEEEAAEEEEEEEEEEEEEEEEAAAAA XA:Z:chr16,+1240848,149M,1;chr16,+1256211,149M,6; MC:Z:150M MD:Z:147G1 RG:Z:NS500158.1 NM:i:1 AS:i:147 XS:i:147 XXXXXXXX:412:YYYYYYYYY:1:11101:10001:10497 163 chr16 1229833 0 150M = 1229894 210 CCAGGCCCTGACCTGTGGAATGTGGTGAGGGGCAGGGTGGACCCCGGCTGGGACTCACCAGGGGCCGCGTAGGCGCGGCTCGCCAGGACGGGCAGCGCCAGCAGCAGCAGATTCAGCATCTGGGGAGCAAGGAGGAGCATCGTGGGCCTG AAAAAEEEEEEEEEEEEAE6EEEAEEEEEEEEEEEEEEAE/EEEEEEEEEEA/AEAEEEEEEEEEAEAE<EEE6A/EEAAAEEEA/EEAAEEAEEE/AAAAEEEEEEEAE/EEEEEEEEEEAEEEEEEAEEEAEE6EAEEAE<</AAA<6 XA:Z:chr16,-1240908,150M,0; MC:Z:149M MD:Z:150 RG:Z:NS500158.1 NM:i:0 AS:i:150 XS:i:150
Note that BWA has suggested an alternative alignment given in the XA tag. When using MergeBamAlignment as in the best practices pipeline, the alignment in XA is chosen. I have tried modifying the --PRIMARY_ALIGNMENT_STRATEGY parameter, but is doesn't change anything.
In the old days before uMAPs, you worked directly with FASTQ files and hence used the primary alignment selected by BWA. What is the motivation for changing that?
Tagged:
Answers
Hi @micknudsen, thanks for your question. Would you mind providing the output from MergeBamAligment for the problematic inputs, and the aligner's name and version (bwa mem vs bwa), with the full command line?
Hi @akovalsk. I lost track of where I got the example from, but I have found another run. When I run
bwa-memdirectly on the FASTQ files, I get the following:but when I use the
gatk4-data-processingversion withMergeBamAlignment, the following alignment is chosen:That is, one of the
XAreads is preferred. I now realize that in both this and the previous example, the read has0mapping quality, and I have not been able to find examples where this is not the case. Still, I find it a little strange why one would prefer one (ambiguous) mapping over the other.Thanks for sharing this info, @micknudsen
MergeBamAlignment uses BestMapqPrimaryAlignmentSelectionStrategy to choose which alignment to use. For paired-end reads, it chooses the alignment that gives the best sum of mapping qualities. In the event that there is a tie, one alignment is chosen randomly (with a fixed seed)
Hi @akovalsk
Looks like the best practices workflows still use MostDistant approach for selecting primary alignments.
Is this still true or has anything changed at best practices since recently?
Hi @SkyWarrior would you mind sharing the context where exactly that screenshot comes from?
Of course. Sorry my bad
https://github.com/gatk-workflows/broad-prod-wgs-germline-snps-indels?files=1
Hi @SkyWarrior that repo is indeed out of date, and users should be looking here instead: https://github.com/gatk-workflows/five-dollar-genome-analysis-pipeline/tree/master
The default alignment strategy for MergeBamAlignment has been BestMapq for some time.
Hi @akovalsk
In this git repo I see that
https://github.com/gatk-workflows/five-dollar-genome-analysis-pipeline/blob/master/tasks/Alignment.wdl
still uses MostDistant as the strategy?
Can you clarify that?
Also is this BestMapq for hg38 based alignments only or does that cover hg19 as well?
command <<< set -o pipefail set -e # set the bash variable needed for the command-line bash_ref_fasta=~{reference_fasta.ref_fasta} # if reference_fasta.ref_alt has data in it, if [ -s ~{reference_fasta.ref_alt} ]; then java -Xms1000m -Xmx1000m -jar /usr/gitc/picard.jar \ SamToFastq \ INPUT=~{input_bam} \ FASTQ=/dev/stdout \ INTERLEAVE=true \ NON_PF=true | \ /usr/gitc/~{bwa_commandline} /dev/stdin - 2> >(tee ~{output_bam_basename}.bwa.stderr.log >&2) | \ java -Dsamjdk.compression_level=~{compression_level} -Xms1000m -Xmx1000m -jar /usr/gitc/picard.jar \ MergeBamAlignment \ VALIDATION_STRINGENCY=SILENT \ EXPECTED_ORIENTATIONS=FR \ ATTRIBUTES_TO_RETAIN=X0 \ ATTRIBUTES_TO_REMOVE=NM \ ATTRIBUTES_TO_REMOVE=MD \ ALIGNED_BAM=/dev/stdin \ UNMAPPED_BAM=~{input_bam} \ OUTPUT=~{output_bam_basename}.bam \ REFERENCE_SEQUENCE=~{reference_fasta.ref_fasta} \ PAIRED_RUN=true \ SORT_ORDER="unsorted" \ IS_BISULFITE_SEQUENCE=false \ ALIGNED_READS_ONLY=false \ CLIP_ADAPTERS=false \ MAX_RECORDS_IN_RAM=2000000 \ ADD_MATE_CIGAR=true \ MAX_INSERTIONS_OR_DELETIONS=-1 \ PRIMARY_ALIGNMENT_STRATEGY=MostDistant \ PROGRAM_RECORD_ID="bwamem" \ PROGRAM_GROUP_VERSION="~{bwa_version}" \ PROGRAM_GROUP_COMMAND_LINE="~{bwa_commandline}" \ PROGRAM_GROUP_NAME="bwamem" \ UNMAPPED_READ_STRATEGY=COPY_TO_TAG \ ALIGNER_PROPER_PAIR_FLAGS=true \ UNMAP_CONTAMINANT_READS=true \ ADD_PG_TAG_TO_READS=false grep -m1 "read .* ALT contigs" ~{output_bam_basename}.bwa.stderr.log | \ grep -v "read 0 ALT contigs" # else reference_fasta.ref_alt is empty or could not be found else exit 1; fi >>> runtime { docker: "us.gcr.io/broad-gotc-prod/genomes-in-the-cloud:2.4.1-1540490856" preemptible: preemptible_tries memory: "14 GiB" cpu: "16" disks: "local-disk " + disk_size + " HDD" } output { File output_bam = "~{output_bam_basename}.bam" File bwa_stderr_log = "~{output_bam_basename}.bwa.stderr.log" } } task SamSplitter { input { File input_bam Int n_reads Int preemptible_tries Int compression_level } Float unmapped_bam_size = size(input_bam, "GiB") # Since the output bams are less compressed than the input bam we need a disk multiplier that's larger than 2. Float disk_multiplier = 2.5 Int disk_size = ceil(disk_multiplier * unmapped_bam_size + 20)