Best Of

Hi @jipvdinter,

We actually no longer use the shifted blacklisted sites. That's an older file that we just forgot to remove from the WDL. We now liftover the unfiltered VCF back to the regular reference before the filtering step which is why we no longer need the blacklisted sites in the shifted reference. I'll make a ticket to clean this up. But just so you know, the sites in the non-shifted blacklist are the correct ones. Thanks for catching this!

@jipvdinter You should be able to do that. The haplogroup is output from Haplochecker task (it's called major_haplogroup). You could use that haplogroup to filter your own list of variants. We use GATK's VariantFiltration tool with the mask argument, but I'm sure there are other tools out there that can filter sites based on an input VCF.

@deena_b

I believe Genotype Given Alleles mode was renamed at some point and the WDLs became out of date. That's why you couldn't find documentation on those arguments. I'll make a ticket to update the WDL, thanks for catching this!

The reason those arguments are inside of ~{"--argument " + value} is because they are optional inputs. This is a feature in WDL so that if you declare a value that whole section inside the curly braces will evaluate to the two strings concatenated together --argument declared_value. If no value is declared it won't put anything there for everything in the curly braces.

The bamout is an argument that you can set to be true or false. It will output a bam of the reads that Mutect2 has realigned around a variant. This can be expensive (because it generates a new BAM), but is very valuable for debugging or understanding what Mutect2 "saw" while making a call, which is why we don't have it on by default all the time.

@deena_b

5. I get that -L is used, for example, in CallShiftedMt to only output variants inside the control region, but what does this part of the code that defines the m2_extra_args variable do: select_first([m2_extra_args, ""])?
   Thank you

I noticed that as well. If you go to the preconfigured Terra pipeline you will see an input section where you can write additional custom GATK Mutect2 command line arguments. I assume this will apply those custom args to the CallMT/CallShiftedMT tasks of the AlignAndCall section where Mutect2 --mitochondria-mode is running.

Hi All,
I wrote the summary below to help myself understand the MitochondriaPipeline steps and thought it might be useful to others.

Note: the nomenclature in these scripts can get confusing, since there are so many steps. Here are three points that are important to grasp:

• To visualize anything in this script called "shifted", you can think of the 16569 nt linear MT DNA being cut near the middle and the latter half being concatonated in front of the first half. The cut was made after nt 8000, meaning that nt 8001 is the new beginning and what was nt #1 is now nt #8597 (16596 - 8000 + 1)
• You can think of the "control" region as the area around a bacteria's origin of replication (aka ori). The control region spans an area a little before and a little after the #1 position of the hg38 reference MT gene and that location is in the middle of the shifted reference. In this case a 'little' is ~500 nt on either side. The control region in this mito pipeline is 1,119 nt long. The "non-control" region is the rest of the MT gene
• Note that the control region is much smaller than the 'cut' that was used to make the MT_shifted reference

Pipeline

1. Use GATK's best practice to pre-process (including aligning) original fastqs with a WGS ref (hg38 is recomended for human)
2. Subset to only keep chrM mapped reads (use PrintReads)
3. Revert chrM mapped reads from aligned bam to unaligned bam (use RevertSam)
4. Convert unaligned bam to fastq (use SamToFastq)
5. Use bwa mem to separately map the fastq files to an MT_reference and an MT_shifted_reference fasta
6. Merge MT_reference aligned bam with unaligned bam from step 3, likewise, merge MT_shifted aligned bam with unaligned bam (use MergeBamAlignment)
7. MarkDuplicates in each aligned bam
8. Coordinate Sort each aligned bam (use SortSam)
9. CollectWgsMetrics for the MT aligned bam
10. Mark reads that are likely contaminants in the MT aligned bam (use haplocheck - note this is not a gatk tool)
    1. haplocheck leverages the mitochondrial phylogeny to detect contamination
11. Call variants in non-control region using the MT bam & in the control region using the MT_shifted bam (use Mutect2)
12. Use LiftoverVcf to return the MT_shifted vcf variant calls back to a standard numbering system (e.g. hg38 reference numbers)
13. Merge the variant calls from the non-control region with the variant calls of the control region (using MergeVcfs)
14. Generate merged statistics using MergeMutectStats
15. Use FilterMutectCalls with these filters
    1. --max-alt-allele-count
    2. --autosomal-coverage
    3. --min-allele-fractio
    4. --f-score-beta
    5. --contamination-estimate
    6. --mitochondria-mode
16. Filter out blacklisted_sites with VariantFiltration
17. Determine the read depth at every base (per_base_coverage) from the standard and shifted bam files using CollectHsMetrics and non_control or control interval_lists