How to Map reads to a reference with alternate contigs like GRCh38

IMPORTANT: This is the legacy GATK documentation. This information is only valid until Dec 31st 2019. For latest documentation and forum click here

created by shlee

on 2016-07-19

Document is in BETA. It may be incomplete and/or inaccurate. Post suggestions to the Comments section and be sure to read about updates also within the Comments section.

File header for -j alignment (alt-handling disabled) for example purposes. We do not perform this type of alignment in the tutorial. Notice the absence of any special tags in the header.

File header for alt-aware alignment after merging with MergeBamAlignment. We use this step in the next section. Again, notice the absence of any special tags in the header.

☞ What is the pa tag?

For BWA v0.7.15, but not v0.7.13, ALT loci alignment records that can align to both the primary assembly and alternate contig(s) will have a pa tag on the primary assembly alignment. For example, read chr19_KI270866v1_alt_4hetvars_26518_27047_0:0:0_0:0:0_931 of the ALTALT sample has five alignment records only three of which have the pa tag as shown below.

A brief description of each of the five alignments, in order:

1. First in pair, primary alignment on the primary assembly; AS=146, pa=0.967
2. First in pair, supplementary alignment on the alternate contig; AS=151
3. Second in pair, primary alignment on the primary assembly; AS=120; pa=0.795
4. Second in pair, supplementary alignment on the primary assembly; AS=54; pa=0.358
5. Second in pair, supplementary alignment on the alternate contig; AS=151

The pa tag measures how much better a read aligns to its best alternate contig alignment versus its primary assembly (pa) alignment. Specifically, it is the ratio of the primary assembly alignment score over the highest alternate contig alignment score. In our example we have primary assembly alignment scores of 146, 120 and 54 and alternate contig alignment scores of 151 and again 151. This gives us three different pa scores that tag the primary assembly alignments: 146/151=0.967, 120/151=0.795 and 54/151=0.358.

In our tutorial's workflow, MergeBamAlignment may either change an alignment's pa score or add a previously unassigned pa score to an alignment. The result of this is summarized as follows for the same alignments.

1. pa=0.967 --MergeBamAlignment--> same
2. none --MergeBamAlignment--> assigns pa=0.967
3. pa=0.795 --MergeBamAlignment--> same
4. pa=0.358 --MergeBamAlignment--> replaces with pa=0.795
5. none --MergeBamAlignment--> assigns pa=0.795

If you want to retain the BWA-assigned pa scores, then add the following options to the workflow commands in section 4.

• For RevertSam, add ATTRIBUTE_TO_CLEAR=pa.
• For MergeBamAlignment, add ATTRIBUTES_TO_RETAIN=pa.

In our sample set, after BWA-MEM alignment ALTALT has 1412 pa-tagged alignment records, PAALT has 805 pa-tagged alignment records and PAPA has zero pa-tagged records.

back to top

4. Add read group information, preprocess to make a clean BAM and call variants

The initial alignment file is missing read group information. One way to add that information, which we use in production, is to use MergeBamAlignment. MergeBamAlignment adds back read group information contained in an unaligned BAM and adjusts meta information to produce a clean BAM ready for pre-processing (see Tutorial#6483 for details on our use of MergeBamAlignment). Given the focus here is to showcase BWA-MEM's alt-handling, we refrain from going into the details of all this additional processing. They follow, with some variation, the PairedEndSingleSampleWf pipeline detailed here.

Remember these are simulated reads with simulated base qualities. We simulated the reads in a manner that only introduces the planned mismatches, without any errors. Coverage is good at roughly 35x. All of the base qualities for all of the reads are at I, which is, according to this page and this site, an excellent base quality score equivalent to a Sanger Phred+33 score of 40. We can therefore skip base quality score recalibration (BQSR) since the reads are simulated and the dataset is not large enough for recalibration anyway.

Here are the commands to obtain a final multisample variant callset. The commands are given for one of the samples. Process each of the three samples independently in the same manner [4.1–4.6] until the last GenotypeGVCFs command [4.7].

[4.1] Create unmapped uBAM

java -jar picard.jar RevertSam \ I=altalt_bwamem.sam O=altalt_u.bam \ ATTRIBUTE_TO_CLEAR=XS ATTRIBUTE_TO_CLEAR=XA

[4.2] Add read group information to uBAM

java -jar picard.jar AddOrReplaceReadGroups \ I=altalt_u.bam O=altalt_rg.bam \ RGID=altalt RGSM=altalt RGLB=wgsim RGPU=shlee RGPL=illumina

[4.3] Merge uBAM with aligned BAM

java -jar picard.jar MergeBamAlignment \ ALIGNED=altalt_bwamem.sam UNMAPPED=altalt_rg.bam O=altalt_m.bam \ R=chr19_chr19_KI270866v1_alt.fasta \ SORT_ORDER=unsorted CLIP_ADAPTERS=false \ ADD_MATE_CIGAR=true MAX_INSERTIONS_OR_DELETIONS=-1 \ PRIMARY_ALIGNMENT_STRATEGY=MostDistant \ UNMAP_CONTAMINANT_READS=false \ ATTRIBUTES_TO_RETAIN=XS ATTRIBUTES_TO_RETAIN=XA

[4.4] Flag duplicate reads

java -jar picard.jar MarkDuplicates \ INPUT=altalt_m.bam OUTPUT=altalt_md.bam METRICS_FILE=altalt_md.bam.txt \ OPTICAL_DUPLICATE_PIXEL_DISTANCE=2500 ASSUME_SORT_ORDER=queryname

[4.5] Coordinate sort, fix NM and UQ tags and index for clean BAM As of Picard v2.7.0, released October 17, 2016, SetNmAndUqTags is no longer available. Use SetNmMdAndUqTags instead.

set -o pipefail java -jar picard.jar SortSam \ INPUT=altalt_md.bam OUTPUT=/dev/stdout SORT_ORDER=coordinate | \ java -jar $PICARD SetNmAndUqTags \ INPUT=/dev/stdin OUTPUT=altalt_snaut.bam \ CREATE_INDEX=true R=chr19_chr19_KI270866v1_alt.fasta

[4.6] Call SNP and indel variants in emit reference confidence (ERC) mode per sample using HaplotypeCaller

java -jar GenomeAnalysisTK.jar -T HaplotypeCaller \ -R chr19_chr19_KI270866v1_alt.fasta \ -o altalt.g.vcf -I altalt_snaut.bam \ -ERC GVCF --max_alternate_alleles 3 --read_filter OverclippedRead \ --emitDroppedReads -bamout altalt_hc.bam

[4.7] Call genotypes on three samples

java -jar GenomeAnalysisTK.jar -T GenotypeGVCFs \ -R chr19_chr19_KI270866v1_alt.fasta -o multisample.vcf \ --variant altalt.g.vcf --variant altpa.g.vcf --variant papa.g.vcf

The altalt_snaut.bam, HaplotypeCaller's altalt_hc.bam and the multisample multisample.vcf are ready for viewing on IGV.

Before getting into the results in the next section, we have minor comments on two filtering options.

In our tutorial workflows, we turn off MergeBamAlignment's UNMAP_CONTAMINANT_READS option. If set to true, 68 reads become unmapped for PAPA and 40 reads become unmapped for PAALT. These unmapped reads are those reads caught by the UNMAP_CONTAMINANT_READS filter and their mates. MergeBamAlignment defines contaminant reads as those alignments that are overclipped, i.e. that are softclipped on both ends, and that align with less than 32 bases. Changing the MIN_UNCLIPPED_BASES option from the default of 32 to 22 and 23 restores all of these reads for PAPA and PAALT, respectively. Contaminants are obviously absent for these simulated reads. And so we set UNMAP_CONTAMINANT_READS to false to disable this filtering.

HaplotypeCaller's --read_filter OverclippedRead option similarly looks for both-end-softclipped alignments, then filters reads aligning with less than 30 bases. The difference is that HaplotypeCaller only excludes the overclipped alignments from its calling and does not remove mapping information nor does it act on the mate of the filtered alignment. Thus, we keep this read filter for the first workflow. However, for the second and third workflows in section 6, tutorial8017toSE and tutorial8017postalt, we omit the --read_filter Overclipped option from the HaplotypeCaller command. We also omit the --max_alternate_alleles 3 option for simplicity.

back to top

5. How can I tell whether I should consider an alternate haplotype?

Here is a command to select out the intentional variant sites that uses SelectVariants:

java -jar GenomeAnalysisTK.jar -T SelectVariants \ -R chr19_chr19_KI270866v1_alt.fasta \ -V multisample.vcf -o multisample_selectvariants.vcf \ -L chr19:34,383,500 -L chr19:34,389,485 -L chr19:34,391,800 -L chr19:34,392,600 \ -L chr19_KI270866v1_alt:32,700 -L chr19_KI270866v1_alt:38,700 \ -L chr19_KI270866v1_alt:41,700 -L chr19_KI270866v1_alt:42,700 \ -L chr19:34,383,486 -L chr19_KI270866v1_alt:32,714

back to top

6. My locus includes an alternate haplotype. How can I call variants on alt contigs?

If you want to call variants on alternate contigs, consider additional data processing that overcome the following problems.

• Loss of alignments from filtering of overclipped reads.
• HaplotypeCaller's filtering of alignments whose mates map to another contig. Alt-handling produces many of these types of reads on the alternate contigs.
• Zero MAPQ scores for alignments that map to two or more alternate contigs. HaplotypeCaller excludes these types of reads from contributing to evidence for variation.

Let us talk about these in more detail.

Ideally, if we are interested in alternate haplotypes, then we would have ensured we were using the most up-to-date analysis reference genome sequence with the latest patch fixes. Also, whatever approach we take to align and preprocess alignments, if we filter any reads as putative contaminants, e.g. with MergeBamAlignment's option to unmap cross-species contamination, then at this point we would want to fish back into the unmapped reads pool and pull out those reads. Specifically, these would have an SA tag indicating mapping to the alternate contig of interest and an FT tag indicating the reason for unmapping was because MergeBamAlignment's UNMAP_CONTAMINANT_READS option identified them as cross-species contamination. Similarly, we want to make sure not to include HaplotypeCaller's --read_filter OverclippedRead option that we use in the first workflow.

Including the mate-mapped-to-other-contig alignments bolsters the variant call qualities for the two SNPs HaplotypeCaller calls on the alternate locus. We see the AD allele depths much improved for ALTALT and PAALT. Corresponding to the increase in reads, the GQ genotype quality and the QUAL score (highlighted in red) indicate higher qualities. For example, the QUAL scores increase from 332 and 289 to 2166 and 1764, respectively. We also see that one of the genotype calls changes. For sample ALTALT, we see a previous no call is now a homozygous reference call (highlighted in blue). This hom-ref call is further from the truth than not having a call as the ALTALT sample should not have coverage for this region in the primary assembly.

For our example data, tutorial8017's callset subset for the primary assembly and tutorial8017_toSE's callset subset for the alternate contigs together appear to make for a better callset.

What solution did we apply? As the workflow's name toSE implies, this approach converts paired reads to single end reads. Specifically, this approach takes the processed and coordinate-sorted BAM from the first workflow and removes the 0x1 paired flag from the alignments. Removing the 0x1 flag from the reads allows HaplotypeCaller to consider alignments whose mates map to a different contig. We accomplish this using a modified script of that presented in Biostars post https://www.biostars.org/p/106668/, indexing with Samtools and then calling with HaplotypeCaller as follows. Note this workaround creates an invalid BAM according to ValidateSamFile. Also, another caveat is that because HaplotypeCaller uses softclipped sequences, any overlapping regions of read pairs will count twice towards variation instead of once. Thus, this step may lead to overconfident calls in such regions.

Remove the 0x1 bitwise flag from alignments

samtools view -h altalt_snaut.bam | gawk '{printf "%s\t", $1; if(and($2,0x1)) {t=$2-0x1}else{t=$2}; printf "%s\t" , t; for (i=3; i<NF; i++){printf "%s\t", $i} ; printf "%s\n",$NF}'| samtools view -Sb - > altalt_se.bam

Index the resulting BAM

samtools index altalt_se.bam

Call variants in -ERC GVCF mode with HaplotypeCaller for each sample

java -jar GenomeAnalysisTK.jar -T HaplotypeCaller \ -R chr19_chr19_KI270866v1_alt.fasta \ -I altalt_se.bam -o altalt_hc.g.vcf \ -ERC GVCF --emitDroppedReads -bamout altalt_hc.bam

Finally, use GenotypeGVCFs as shown in section 4's command [4.7] for a multisample variant callset. Tutorial8017toSE calls 68 variant sites--66 on the primary assembly and two on the alternate contig.

6.2 Variant calls for tutorial8017postalt

BWA's postalt-processing requires the query-grouped output of BWA-MEM. Piping an alignment step with postalt-processing is possible. However, to be able to compare variant calls from an identical alignment, we present the postalt-processing as an add-on workflow that takes the alignment from the first workflow.

The command uses the bwa-postalt.js script, which we run through k8, a Javascript execution shell. It then lists the ALT index, the aligned SAM altalt.sam and names the resulting file > altalt_postalt.sam.

k8 bwa-postalt.js \ chr19_chr19_KI270866v1_alt.fasta.alt \ altalt.sam > altalt_postalt.sam

Postalt-processing lowers the MAPQ of primary locus alignments in the highly divergent region that map better to the alt locus. You can see this as a subtle change in the IGV screenshot. After postalt-processing we see an increase in white zero MAPQ reads in the highly divergent region of the primary locus for ALTALT and PAALT. For ALTALT, this effectively cleans up the variant calls in this region at chr19:34,391,800 and chr19:34,392,600. Previously for ALTALT, these calls contained some reads: 4 and 25 for the first workflow and 0 and 28 for the second workflow. After postalt-processing, no reads are considered in this region giving us ./.:0,0:0:.:0,0,0 calls for both sites.

What we omit from examination are the effects of postalt-processing on decoy contig alignments. Namely, if an alignment on the primary assembly aligns better on a decoy contig, then postalt-processing discounts the alignment on the primary assembly by assigning it a zero MAPQ score.

To wrap up, here are the number of variant sites called for the three workflows. As you can see, this last workflow calls the most variants at 95 variant sites, with 62 on the primary assembly and 33 on the alternate contig.

Workflow total on primary assembly on alternate contig tutorial_8017 56 54 2 tutorial_8017_toSE 68 66 2 tutorial_8017_postalt 95 62 33

back to top

7. Related resources

• For WDL scripts of the workflows represented in this tutorial, see the GATK WDL scripts repository.
• To revert an aligned BAM to unaligned BAM, see Section B of Tutorial#6484.
• To simulate reads from a reference contig, see Tutorial#7859.
• Dictionary entry Reference Genome Components reviews terminology that describe reference genome components.
• The GATK resource bundle provides an analysis set GRCh38 reference FASTA as well as several other related resource files.
• As of this writing (August 8, 2016), the SAM format ALT index file for GRCh38 is available only in the x86_64-linux bwakit download as stated in this bwakit README. The hs38DH.fa.alt file is in the resource-GRCh38 folder. Rename this file's basename to match that of the corresponding reference FASTA.
• For more details on MergeBamAlignment features, see Section 3C of Tutorial#6483.
• For details on the PairedEndSingleSampleWorkflow that uses GRCh38, see here.
• See here for VCF specifications.

back to top

Updated on 2017-01-18

From RMulet on 2016-10-19

First of all, thanks for this wonderful, comprehensive tutorial. I want to bring to your attention that the post-processing step in 6.2 does not work in the form presented here: the script expects the .bam alignment following the ALT index, not a .fasta file containing the reference genome. Thus, it should be as follows: bwa-postalt.js chr19_chr19_KI270866v1_alt.fasta.alt altalt.sam > altalt_postalt.sam

As stupid as it seems, it has taken me a while to figure out why it did not work…

From Sheila on 2016-10-20

@RMulet

Hi,

Thanks for the suggestion and praise. I will let my colleague @shlee know :smile:

-Sheila

From shlee on 2016-10-21

Hi @RMulet,

Thanks for the compliment and thanks for bringing this up. I just want to confirm, are you also using BWA v0.7.15?

From shlee on 2016-10-27

Thanks @RMulet for pointing this out. I corrected this in the body of the tutorial on this page. I should mention here that the [matching WDL script](https://github.com/broadinstitute/wdl/blob/develop/scripts/tutorials/gatk/tutorial_8017_postalt.wdl) remains correct and does not include the extraneous fasta file in the command.

From RMulet on 2016-10-28

You are welcome! Sorry for not replying to the previous comment before, but I did not receive any notification. I guess it’s no longer necessary, as it has been modified now, but I am indeed using BWA v0.7.15.

From alexey_larionov on 2016-11-07

Thank you for the tutorial.

Could you please explain in more details what exactly is done by SetNmAndUqTags ?

This tool is not described in the current Picard page at http://broadinstitute.github.io/picard/index.html

Alexey

From shlee on 2016-11-09

Hi @alexey_larionov,

SetNmAndUqTags is a new tool. I can refer you to two sources of information.

—-

(1) Although the website documentation does not yet reflect the tool specs, you can call them up from the commandline with the `-h` option.

```

$ java -jar picard.jar SetNmAndUqTags -h

```

This gives:

```

USAGE: SetNmAndUqTags [options]

Fixes the UQ and NM tags in a SAM file. This tool takes in a SAM or BAM file (sorted by coordinate) and calculates the NM and UQ tags by comparing with the reference.

This may be needed when MergeBamAlignment was run with SORT_ORDER different from ‘coordinate’ and thus could not fix

these tags then.

Usage example:

java -jar picard.jar SetNmAndUqTags \ I=sorted.bam \ O=fixed.bam \

Version: 2.5.0(2c370988aefe41f579920c8a6a678a201c5261c1_1466708365)

```

—-

(2) Also, I’ve given an overview of what the tool does in the context of a WDL pipeline. The documentation section pertaining to this tool is at [http://broadinstitute.github.io/picard/command-line-overview.html#SetNmAndUqTags](http://broadinstitute.github.io/picard/command-line-overview.html#SetNmAndUqTags). In this pipeline, because of the way it is setup, after the MergeBamAlignment step,

> Note the NM tag values will be incorrect at this point and the UQ tag is absent. Update and addition of these are dependent on coordinate sort order. Specifically, the script uses a separate SortAndFixTags task to fix NM tags and add UQ tags.

MergeBamAlignment adds the UQ tag to records given a coordinate sorted BAM. The tag gives the

> Phred likelihood of the segment, conditional on the mapping being correct.

If MergeBamAlignment is given a query-group sorted BAM, then the tool omits addition of the UQ tag. Also, the NM tag values will be incorrect. In this situation, SetNmAndUqTags corrects the NM tag values and adds the UQ tags.

I hope this clarifies your question.

From EADG on 2016-11-30

HI, @shlee

The Documentation about that tool is still missing at the picard docs :(.

From shlee on 2016-11-30

Hi @EADG, we have on our to do list to update the Picard docs. I hope you have enough information from the thread above to continue with your work.

From shlee on 2016-11-30

Actually @EADG, as of Picard v2.7.0, SetNmAndUqTags is now SetNm**Md**AndUqTags.

```

java -jar $PICARD SetNmMdAndUqTags -h

```

Gives:

```

USAGE: SetNmMdAndUqTags [options]

Documentation: http://broadinstitute.github.io/picard/command-line-overview.html#SetNmMdAndUqTags

Fixes the NM, MD, and UQ tags in a SAM file. This tool takes in a SAM or BAM file (sorted by coordinate) and calculates the NM, MD, and UQ tags by comparing with the reference.

This may be needed when MergeBamAlignment was run with SORT_ORDER different from ‘coordinate’ and thus could not fix

these tags then.

Usage example:

java -jar picard.jar SetNmMDAndUqTags \ I=sorted.bam \ O=fixed.bam \

Version: 2.7.1-SNAPSHOT

Options:

—help

-h Displays options specific to this tool.

—stdhelp

-H Displays options specific to this tool AND options common to all Picard command line tools.

—version Displays program version.

INPUT=File

I=File The BAM or SAM file to fix. Required.

OUTPUT=File

O=File The fixed BAM or SAM output file. Required.

IS_BISULFITE_SEQUENCE=Boolean Whether the file contains bisulfite sequence (used when calculating the NM tag). Default value: false. This option can be set to ‘null’ to clear the default value. Possible values: {true, false}

```

From EADG on 2016-12-01

@shlee Thank you!

From shlee on 2017-01-18

Updates as of January 18, 2017

• Newer versions of Picard going forward will now copy over AH tags as well as allow you to specify which @SQ lines header--from the uBAM or aligned BAM--to use in the merged output.
• Newer versions of GATK will allow disabling of the BadMateFilter in HaplotypeCaller and thereby enable calling variants in GRCh38 alternate contigs after BWA-MEM alt-aware alignment without having to strip the 0x1 SAM flag off of reads.

From Wang_GC on 2017-07-14

Hi,

Can I replace the “Mark Duplicates” part of best practices with 4.1-4.5 and follow the rest of best practices after mapping reads to reference?

From shlee on 2017-07-14

Hi @Wang_GC,

The best practice germline workflow for GRCh38 is outlined in .

My workflow above is missing BQSR (because the reads are simulated to have no errors) and you will need to include this step in your pipeline. This tutorial is meant to illustrate the effects of alt-aware alignment to a reference with alternate haplotypes. I’m sure there are inefficiencies that could be done better for real analyses. For example, you should already have in hand a uBAM that you can add read groups to and you would not need to revert your alignment. I’m told that the decision to include MergeBamAlignment in your workflow depends on your downstream analyses.

From Wang_GC on 2017-07-28

@shlee

Thanks for your reply! I will rerun my scripts with the recommended pipeline.

From ac67479 on 2017-10-31

Dear Mrs Shlee,

I was wondering whether there is any pipeline proposed for RNA-seq differential expression analysis which takes into account these alt loci on GRCh38 ? We are using hisat2 on a regular basis in the lab, so we would like to find a pipeline to use the new reference and take into account these alt loci, preferentially with hisat2.

Regards,

Alexandre Coudray

From shlee on 2017-10-31

Hi @ac67479,

The GATK focuses on variant discovery and not differential expression. We are currently developing, in GATK4, functionality towards differential copy number variant detection. However, these are not applicable to detecting differential RNA expression.

I’ve heard that HISAT2 is an aligner that can incorporate given variants into its graph-based assembly. This certainly sounds promising towards taking into account ALT loci in GRCh38. I am not familiar with the details of using HISAT2 as our workflows recommend STAR for RNA-Seq alignment and BWA for DNA-Seq alignment.

If you are interested in calling variants from HISAT2 alignments with GATK, then you should first test to see if the output format is compatible with the GATK workflow for germline variant discovery. I am actually interested in knowing what the result of this is.

Keep in mind there are limitations to calling variants from RNA-Seq data compared to DNA-Seq data.

P.S. You can just call me @shlee. If you insist on adding a title to my moniker, I’d prefer `Dr`. Thanks.

From shlee on 2017-10-31

Actually @ac67479, I wonder if you’d be willing to share your HISAT2 alignment data (assuming non-human) and I could check it out on my spare time. We have an FTP upload, with instructions at , if you would be willing.

From KyuriChoi on 2018-02-07

@shlee Why new bundle is composed with ALT contig fasta file, not using GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz? using “no alt” fasta is more easy way, isn’t it?

From SkyWarrior on 2018-02-07

@KyuriChoi

I guess you totally missed the point.

From Sheila on 2018-02-13

@KyuriChoi

Hi,

The reasons for using alt contig version are listed in this tutorial, but if you feel strongly about using the non-alt version, you can certainly do so.

-Sheila

From wilma on 2018-03-09

Thanks for this very useful and comprehensive article!

Assume for a second I needed a unified approach for processing samples to BAM, which leaves the decision to call variants on the ALT contigs or not for later. In that case 0×1 removal can’t be part of that pipeline, because it changes the BAM files irreversibly. My question is: is it always safe to include `k8 bwa-postalt.js`, i.e. can the output be used for “ALT calling” (after 0×1 removal) and “non-ALT calling” (no further processing needed)?

Thanks

From Sheila on 2018-03-14

@wilma

Hi,

Someone will get back to you soon.

-Sheila

From shlee on 2018-03-14

Hi @wilma,

In GATK4, you can disable a filter and call on the alt reads just as if you had removed the pair flag. This filter is disabled with the `—disable-read-filter MateOnSameContigOrNoMappedMateReadFilter` option. See my March 5 reply to KlausNZ at https://gatkforums.broadinstitute.org/gatk/discussion/comment/46649#Comment_46649.

The decision to include `k8 bwa-postalt.js` processing is left to researchers. In the tutorial, I tried to illustrate the effects of including it or not, to illustrate two modes in which our production pipelines were run at some point. As I note in the tutorial, our production currently generates data WITHOUT the `k8 bwa-postalt.js` processing and only makes calls on the primary assembly. Again, it is left to researchers to decide what type of calling is best for their research aims.

From gilh on 2018-03-28

Hi,

I see that this section is still in beta for over a year, and the hg38 bundle in the google drive is labelled as “v0”. I am just wondering if this protocol is still the correct one, or whether there are updates.

Specifically, I would like to use the mapping with ALT reads to tumor exome data in order to detect somatic variants, and I want to know whether this protocol is fitting.

Thank you,

Gil

From Sheila on 2018-04-03

@gilh

Hi Gil,

This document is the one to use. I will check and see why it is in beta, but there is no updated/new document that I am aware of.

-Sheila

From sergey_ko13 on 2018-05-24

dear Shlee,

is BWA postalt processing obligatory when alt-aware alignment is done?

(we are not searching for variants in alternative contigs)

if we skip postalt processing – does mutect’s disable of MateOnSameContigOrNoMappedMateReadFilter still make sense?

Thank You!

From shlee on 2018-05-24

Hi @sergey_ko13,

As I say above in other replies, whether or not you perform post-alt processing is up to you and will depend on your research aims. Perhaps you will find skimming [this post](https://software.broadinstitute.org/gatk/blog?id=8180) helpful.

From Phuong_Nguyen_1 on 2018-05-30

> “As of this writing (August 8, 2016), the SAM format ALT index file for GRCh38 is available only in the x86_64-linux bwakit download as stated in this bwakit README. The hs38DH.fa.alt file is in the resource-GRCh38 folder“

>

Hi,

I just want to know is there any ALT index available for current genome version (GRCh38p12) ?

Is there anyway to create that ALT index file on our own? If yes, could you please share the detail steps.

Thank you very much!

Phuong.

From shlee on 2018-05-30

Hi @Phuong_Nguyen_1. I checked with Heng Li and he says:

> If you don’t run the javascript, you can put all the ALT contig names into a file, e.g. `prefix.alt`, one name per line. If you want to run the javascript, you have to convert NCBI’s alignment in the gff3 format to SAM, which is complicated.

The javascript refers to the `bwa-postalt.js script` from the above tutorial that gets run through k8. I hope this is helpful.

From Phuong_Nguyen_1 on 2018-06-05

Thank you very much, I guess that I should stick with .alt index from hs38DH., I hope to see some scripts or tutorial which show us how to make .alt index for any version of GRCh38.

Additionally, I would like to ask another point: Are those 3 mapping processes applicable for Mutect2 somatic variant calling, is there any difference? In particular, did I need to perform “remove 0×1 flag” step before running Mutect2?

Thank you very much in advance.

From shlee on 2018-06-07

Hi @Phuong_Nguyen_1,

I hope to see such scripts too. One of the features of the .alt file that might be of interest to you is that it is in SAM format and therefore parseable with samtools.

Whichever mapping process you decide on, just note that you should try to match the analysis pipeline for your panel of normals samples, if you are using a PoN.

As for your question

> did I need to perform “remove 0×1 flag” step before running Mutect2?

For GATK4 this is no longer necessary as you can now disable the read filter that I was trying to hack around with the “remove 0×1 flag” step. In GATK3 you could not disable this read filter. In GATK4 you need only specify `—disable-read-filter MateOnSameContigOrNoMappedMateReadFilter` in your command.

From sergey_ko13 on 2018-06-18

Hi @shlee ,

I am a bit lost with altalt.sam appers first time in this tutorial only when we talk about postalt processing. Is that just same file as altalt_bwamem.sam? Do I run bwa-postalt.js on alighned sam right before MergeBamAlignment?

Thank You!

Sergey

From shlee on 2018-06-18

Hi @sergey_ko13,

The two files are different. The `altalt.sam` file is the input to k8 bwa-postalt processing. The output is `altalt_postalt.sam`:

```

k8 bwa-postalt.js \ chr19_chr19_KI270866v1_alt.fasta.alt \ altalt.sam > altalt_postalt.sam

```

In the tutorial, yes, postalt processing is run immediately after BWA MEM alignment and before MergeBamAlignment. I’d like to mention that since writing this tutorial I’ve learned that the 1000 Genomes Project , when postalt processing GRCh38 alignments, does NOT perform MergeBamAlignment.

From sergey_ko13 on 2018-06-19

Thank You @shlee . does MergeBamAlignment conflict with postalt alignments or 1000 genomes project just does not comment such choice?

From shlee on 2018-06-19

@sergey_ko13, the project makes no comment and I do not know of anyone who has investigated. I suspect there are some implications for pair assignments. If you find out exactly what is the difference, it would be great for you to share.

From shlee on 2018-09-26

@Adam_U0, I moved your question to https://gatkforums.broadinstitute.org/gatk/discussion/13055/question-on-article-8017-mapping-reads-to-grch38.

From blueskypy on 2019-01-07

hi, @shlee

Thanks for the article! May I ask three questions?

1. if the [Broad public pipeline](https://github.com/gatk-workflows/broad-prod-wgs-germline-snps-indels “Broad public pipeline”) does not do ALT calling, why does it specially says “Reference genome must be Hg38 with ALT contigs”?

2. do you know if the current variant annotation tools would be able to annotate variants in those ALT contigs?

3. Are most of the ALT contigs in MHC regions? If so, I suppose variant calling on ALT contigs would be beneficial for study of immunity-related diseases. Do you know any actual studies where variant calling on ALT contigs is used?

Happy New Year, @shlee !

From shlee on 2019-01-14

Happy New Year to you too @blueskypy!

1. Last I checked, which was a while ago, the Broad GRCh38 pipeline performs alt-aware alignment but does not perform post-alt processing. Alt-aware alignment still requires a reference with ALT contigs.

2. As far as I’m aware, GATK4 tools, including Funcotator, should be able to annotate variants in ALT contigs given annotation files such as GTFs that contain information for features in the ALTs. If you find this is not the case, please let us know.

3. Please see [Blog#8180](https://gatkforums.broadinstitute.org/gatk/discussion/8180/9-takeaways-to-help-you-get-started-with-grch38) for a breakdown of what the ALT contigs represent. The table in section 1 and section 5 should answer your question. Let me recapitulate the table I had made for you here:


I believe there are HLA-specific aligners and callers out there to aid in the study of such loci. Heng Li comments on HLA typing in the BWA repo at https://github.com/lh3/bwa/blob/master/README-alt.md#hla-typing. Also, check out his thoughts in the Preliminary Evaluation section and the Problems and Future Development section.

If you know of any studies where variant calling on ALT contigs is used, please do share.

From blueskypy on 2019-01-18

Thanks so much, @shlee !

From blueskypy on 2019-02-22

hi, @shlee

I’m working on WES. If the target regions of the capturing kit do not include ALT contigs, and I apply `-L target_regions` in BQSR, I suppose the bam files after BQSR won’t contain off-target mapping. It won’t make sense to apply this alt-mapping and calling in this case, and in general to targeted seq where targets do not include ALT contigs. Is my understanding correct?

From blueskypy on 2019-02-28

hi, @shlee, could you please take a look at my question? I really appreciate!

From shlee on 2019-02-28

Hi @blueskypy,

Take a look at Geraldine’s post [here](https://software.broadinstitute.org/gatk/blog?id=23384). As you can see, the GATK Comms Team has reorganized and also expanded (doubled in size in fact!). So, I just wanted to point out that I’m not a part of the Frontline Support Team, who started fielding forum questions since last summer. Sheila was the singular point person previously and so I helped her as I could as the only other Comms member on the GATK side of things. Now I am actually a part of the User Education Team, who develop and write content. I am still active on the forum, but these days typically at the request of the Frontline Support Team for help with some question or another. So I hope you understand I’m not on the forum as much. This is due to the [time-pie constant](https://software.broadinstitute.org/gatk/blog?id=9063). Even if you ping my handle personally, it is the Frontline Support Team who will typically answer.

So, the new world order is you have to wait until someone from the Frontline Support tickets your question for response. One tip is that you’ll have a better chance of someone answering your question quickly if you post to a new thread, instead of to an existing article.

You are asking two questions. For targeted exomes:

> I apply -L target_regions in BQSR, I suppose the bam files after BQSR won’t contain off-target mapping.

and

> It won’t make sense to apply this alt-mapping and calling in this case, and in general to targeted seq where targets do not include ALT contigs. Is my understanding correct?

For your first question, please apply base recalibration to the entirety of your data, including your unmapped reads. I think you will find the gatk-workflows repo’s [gatk4-exome-analysis-pipeline](https://github.com/gatk-workflows/gatk4-exome-analysis-pipeline) helpful here. If you dissect the pipeline, you see BQSR is performed on the entirety of data, not just on targets.

Can you rephrase your second question? It’s not clear to me what you are asking.

From blueskypy on 2019-02-28

Thanks so much @shlee! I’ll address my future questions to the Frontline team.

Regarding my first question, [this post](https://gatkforums.broadinstitute.org/gatk/discussion/11062/when-should-i-restrict-my-analysis-to-specific-intervals “this post”) says the following:

>In addition, there are some processing steps, such as BQSR, that should be restricted to the capture targets in order to eliminate off-target sequencing data, which is uninformative and is a source of noise.

Based on your comment, seems the above statement is wrong. I was wondering that, if it were right, after BQSR, the bam file does not contain reads mapping to ALT contigs anymore if the target region does not include ALT contigs; then there is no need to apply this ALT-mapping and calling any more either.

From blueskypy on 2019-02-28

btw, @shlee, does [gatk4-exome-analysis-pipeline](https://github.com/gatk-workflows/gatk4-exome-analysis-pipeline “gatk4-exome-analysis-pipeline”) perform ALT-mapping and calling?

From shlee on 2019-02-28

@blueskypy,

It makes sense to me as the linked article states to generate the BQSR report from just the padded targets. However, it’s my understanding you should apply the base recalibration itself (the second step) to the entirety of the data. I only quickly glanced at the gatk-workflow’s exome pipeline and it appeared to me that both steps used all data, presumably because (unchecked) the workflow limits downstream variant calling to the targets. I only took a quick glance so I could be mistaken. I hope @bshifaw can clarify. These pipelined workflows are meant to represent GATK recommendations.

BQSR keeps all reads as its function is to adjust base quality scores. To reiterate, you should have all the reads you started with.

Again, your second question is unclear to me. I think the above point about keeping all reads should help you.

From blueskypy on 2019-02-28

hi, @shlee, thanks for the quick response again! I really aprpeciate! I’ve posted a [new question](https://gatkforums.broadinstitute.org/gatk/discussion/23568/how-to-do-bqsr-on-wes-or-other-targeted-sequencing/p1?new=1 “new question”) for the Frontline team.

Regarding my 2nd question, here is what I mean:

If the mapping or variant calling is limited to targets only, and the `target.bed` has only the primary chromosomes, i.e. does not include ALT-contigs, there is no need to apply the approach (ALT-mapping and calling) in this article. is it correct?

From shlee on 2019-02-28

@blueskypy,

It’s my understanding exome target kit targets are empirically designed. If a kit’s targets have been calibrated to work well across populations so that there is some acceptable level of capture across the board no matter the sample, say 95%, and this sort of cutoff suits your research project and your project’s samples, then you STILL need to perform alt-aware alignment if you are in GRCh38. Please see the figure I made for section 4 of [Blog#8180](https://gatkforums.broadinstitute.org/gatk/discussion/8180/9-takeaways-to-help-you-get-started-with-grch38). The white reads are secondary alignments and GATK ignore these. If your analysis is limited to the primary assembly, you should perform your own due-diligence for your project aims for the value-add of the additional post-alt processing as there are differences in the data (that [Tutorial#8017](https://gatkforums.broadinstitute.org/gatk/discussion/8017/how-to-map-reads-to-a-reference-with-alternate-contigs-like-grch38) outlines) even on the primary assembly, e.g. the handling of reads that map to decoys. I touch upon these briefly in section 6.2.

The alternate contigs represent alternate haplotypes of corresponding regions of the primary assembly. It is not unreasonable to assume that exome target kits accommodate some extent of the alternate haplotypes, which can be common in the population, and in this case you might find it preferable to analyze the corresponding target region on the alternate contig instead of the primary assembly. You could do a quick back-of-the-envelope calculation on how much of your target kit’s target baits has corresponding alternate regions by performing a liftover of your bait intervals. You can do this at the [UCSC LiftOver site](https://genome.ucsc.edu/cgi-bin/hgLiftOver), the [NCBI remap site](https://www.ncbi.nlm.nih.gov/genome/tools/remap) and using [Picard LiftOverIntervalList](https://software.broadinstitute.org/gatk/documentation/tooldocs/current/picard_util_LiftOverIntervalList.php). Last I conversed with Heng Li on this topic, the recommendation is to perform all and compare as there are differences between the methods.

P.S. That your target intervals are limited to the primary assembly is an artificial construct of past methods and resources. To take full advantage of all the improvements GRCh38 brings, you should request from your target kit manufacturer corresponding intervals for the alternate haplotypes.

From blueskypy on 2019-02-28

hi, @shlee, thanks again for the reply! I’m interested in both primary chromosomes and ALT-contigs, the problem is the `target.bed` from those commercial exome capturing kits, e.g. Agilent SureSelect V7, specifically designed for hg38, only contains chr1-23, X, Y, M, no other contigs.

Therefore, to not use ALT-mapping and calling, I can use a version of hg38 with no ALT contigs based on [Heng Li’s blog](http://lh3.github.io/2017/11/13/which-human-reference-genome-to-use “Heng Li’s blog”).

But to use ALT-mapping and calling for WES, I wonder if the method 2 in [my post](https://gatkforums.broadinstitute.org/gatk/discussion/23568/how-to-do-bqsr-on-wes-or-other-targeted-sequencing/p1?new=1 “my post”) is the option?