My earlier post detailed several scripts to reprocess vendor supplied BAMs using a methodology based on Broad Institute’s Best Practices. FTDNA Big Y customers attempting to follow this process could have several issues present themselves. Each has a solution, but most result in a small percentage of reads being eliminated. Ideally, FTDNA would release your full set of reads instead of filtering to prevent this.

Here are a list of common issues and the work-around needed to get them through the GATK workflow. This post will be a living document that is updated as new issues are brought to light.

Recipe 1: Fix Broken Read Names

If you see:

INFO 2018-06-15 05:06:19 RevertSam Discarded 10069697 out of 10069697 (100.000%) reads in order to sanitize output.

Then do this:

samtools bam2fq 0629-c0fc0db1-8e0e-4548-aa9c-2098bcfc4e3a.bam | gzip > interleaved_reads.fq.gz

bbmap/repair.sh in=interleaved_reads.fq.gz out1=1.fq out2=2.fq outs=singletons.fq repair

What it does:

The root issue here is that FTDNA has discarded part of the read pair name in their processing. We need to start over from the fundamentals by reverting their BAM to an interleaved FASTQ file. Once we have the raw read data we can use Bbmap’s repair module to reapply the missing /1 and /2 portion of the names. The first two output files are suitable for running into fastq_to_sam.sh. The singletons.fq file contains the unpaired reads caused by FTDNA filtering your BAM to only the chrY regions. BWA MEM will not be able to align them as accurately without their mate, and clean BAMs should not contain them under best practices.

Recipe 2: Fix Broken CIGAR Encoding

If you see:

htsjdk.samtools.SAMFormatException: SAM validation error: ERROR: Read name SN7001368:263:HALJ0ADXX:1:2214:15631:92476, Read CIGAR M operator maps off end of reference

OR

htsjdk.samtools.SAMFormatException: SAM validation error: WARNING: Read name SN7001368:184:H9CKDADXX:1:2203:6512:94819, No M or N operator between pair of I operators in CIGAR

Then do this:

samtools bam2fq 3198-2f06198e-5c31-4ac1-aacf-ce7cce5cf646.bam | gzip > interleaved_reads.fq.gz

bbmap/repair.sh in=interleaved_reads.fq.gz out1=1.fq out2=2.fq outs=singletons.fq repair

What it does:

The original generation of Big Y BAMs, while containing all of your sequencing data, had an issue with CIGAR encoding in aengine. The reads are being encoded as falling off the reference or incorrectly tagging substitution and/or indel events. The solution is basically the same as before. Revert to the raw FASTQ files. The repair is likely not necessary but the overhead over simply splitting the interleaved file is minimal. Unfortuantely, GATK doesn’t seem to provide a command line option to simply ignore the error as the mapping information is being discarded by RevertSam anyway.

Broad Institute released a new major version of their Swiss Army knife of bio-informatics software known as GATK tools on January 9, 2018.  This means it is time to overhaul the Haplogroup-R.org variant discovery workflow to use the latest suite.  The key change as far as workflow is the introduction of GenomicsDB, which vastly improves the performance of joint calling.

Software Tools Needed

The requirements for the workflow include a POSIX operating system, a Java 8 run-time environment, and access to typical development tools in some cases.  The software packages included in the workflow’s processing include:

GATK 4.0.4.0 (Download here.)

BWA 0.7.17 (Available here.)

BBmap (Download here.)

The software also requires a bundle of reference files to align raw sequencing data.  Which reference you choose GRCh37 or GRCh38 is a matter of preference still, but most of the direct to consumer providers have migrated to GRCh38.  Previous versions of this workflow have used the 1000 Genome’s bundle, since it provided the ability to use their aligned samples in the analysis set without fuss.  However, those starting today may wish to use the bundle provided by the Broad Institute.  The following command will clone their hg38 reference to current directory:

wget -m ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38

Hardware Needed

The scripts assume 16GB of RAM is available on the hardware.  BWA MEM only requires 5.5 GB for a single thread, so one could run those processes on lesser hardware and trade a much longer processing time. The Java memory options come from GATK’s sample scripts.  It is possible less than 16GB is needed for a Big Y BAM, but no effort was spent on tuning these parameters.  WGS samples typically need 50 to 120 GB of disc space per intermediate step in the scripts.  Having 1TB of free space for a temporary directory is recommended.

What is new?

One of the major changes in GATK’s Best Practice documentation is in the pre-processing of FASTQ.  Their recommendation in tutorial #6483 is to produce a clean BAM & BAI to use as input into a variant discovery chain.  The process involves generating an unaligned BAM file and marking the sequencing adapters.  The adapters are trimmed and the remaining reads are streamed into BWA MEM to locate on the genome.  Finally, the alignment meta data is merged back with the original unaligned BAM.  The resulting file will have all the original read data and properly adjusted alignment meta data.  Read the linked tutorial for much more detail.

In older versions of this work flow a large number of gVCFs were batched together for genotyping.  This process could take over 24 hours with the full 1,000 gathered samples.  While most project administrators will not be dealing with batches that large, it’s reassuring to see the headroom gained with this optimization.

The Process 

 Script 1a:  Convert FASTQ to unaligned BAM

# USAGE: sh fastq_to_sam.sh <fastq1> <fastq2> <sample_name> <read_group> <platform_unit>
gatk=~/Genomics/gatk-4.0.4.0/gatk
$gatk --java-options "-Xmx8G" FastqToSam \
    -FASTQ=$1 \
    -FASTQ2=$2 \
    -OUTPUT=$3.unaligned.bam \
    -READ_GROUP_NAME=$4 \
    -SAMPLE_NAME=$3 \
    -LIBRARY_NAME=$3 \
    -PLATFORM_UNIT=$5 \
    -PLATFORM=illumina

Script 1b:  Revert BAM to unaligned BAM

# USAGE:  sh revert_bam.sh <sample name>
# Assumes GATK is on the path.  Based on https://gatkforums.broadinstitute.org/gatk/discussion/6484#latest%23top
gatk RevertSam \
    -I=$1.bam \
    -O=$1.unaligned.bam \
    -SANITIZE=true \
    -MAX_DISCARD_FRACTION=0.005 \
    -ATTRIBUTE_TO_CLEAR=XT \
    -ATTRIBUTE_TO_CLEAR=XN \
    -ATTRIBUTE_TO_CLEAR=AS \
    -ATTRIBUTE_TO_CLEAR=OC \
    -ATTRIBUTE_TO_CLEAR=OP \
    -SORT_ORDER=queryname \
    -RESTORE_ORIGINAL_QUALITIES=true \
    -REMOVE_DUPLICATE_INFORMATION=true \
    -REMOVE_ALIGNMENT_INFORMATION=true

Script 2: Prepare the Clean BAM for an Illumina Sample

# USAGE:  sh create_clean_bam.sh <sample name>
# Based on https://software.broadinstitute.org/gatk/documentation/article.php?id=6483
# CONFIG VARIABLES:  Update to match environment
gatk=~/Genomics/gatk-4.0.4.0/gatk
reference=~/Genomics/Reference/GRCh38_full_analysis_set_plus_decoy_hla.fa
tmp_dir=/Volumes/External/tmp

# Mark the Illumina adapters (if present.  The sequencing lab should have removed them
# prior to delivering the results.)
$gatk --java-options "-Xmx8G" MarkIlluminaAdapters \
-I=$1.unmapped.bam \
-O=$1.markilluminaadapters.bam \
-M=$1.markilluminaadapters.metrics.txt \
-TMP_DIR=$tmp_dir

# Perform a piped operation that trims the Illumina adapters from the reads, aligns them
# to the target reference, and creates a clean BAM for variant discovery.
$gatk --java-options "-Xmx8G" SamToFastq \
-I=$1.markilluminaadapters.bam \
-FASTQ=/dev/stdout \
-CLIPPING_ATTRIBUTE=XT -CLIPPING_ACTION=2 -INTERLEAVE=true -NON_PF=true \
-TMP_DIR=$tmp_dir | \
~/Genomics/bwa/bwa mem -M -t 7 -p $reference /dev/stdin | \
$gatk --java-options "-Xmx16G" MergeBamAlignment \
-ALIGNED_BAM=/dev/stdin \
-UNMAPPED_BAM=$1.unmapped.bam \
-OUTPUT=$1.bwa.clean.bam \
-R=$reference -CREATE_INDEX=true -ADD_MATE_CIGAR=true \
-CLIP_ADAPTERS=false -CLIP_OVERLAPPING_READS=true \
-INCLUDE_SECONDARY_ALIGNMENTS=true -MAX_INSERTIONS_OR_DELETIONS=-1 \
-PRIMARY_ALIGNMENT_STRATEGY=MostDistant -ATTRIBUTES_TO_RETAIN=XS \
-TMP_DIR=$tmp_dir

Once the second script has been applied, the unmapped and adapter marked BAMs are no longer necessary and can be removed.  You can always get the unmapped file back by applying the revert_bam.sh script.

Script 3: Recalibrate the read bases and generate a gVCF for the cleaned BAM

# USAGE:  sh prepare_gvcf.sh <sample name>
# CONFIG VARIABLES:  Update to match environmentgatk=~/Genomics/gatk-4.0.4.0/gatk
reference=~/Genomics/Reference/GRCh38/GRCh38_full_analysis_set_plus_decoy_hla.fa
known=~/Genomics/Reference/GRCh38/Mills_and_1000G_gold_standard.indels.b38.primary_assembly.vcf.gz

$gatk --java-options "-Xmx4G" \
	MarkDuplicates -I=$1.bwa.clean.bam -O=$1.dedup.bam -METRICS_FILE=metrics.txt
     
$gatk --java-options "-Xmx8G" \
	BaseRecalibrator -R $reference -O recal_data.table -I $1.dedup.bam \
	--known-sites $known

$gatk --java-options "-Xmx8G" ApplyBQSR \
  -R $reference \
  -I $1.dedup.bam \
  --bqsr-recal-file recal_data.table \
  -O $1.recal.bam
  
$gatk --java-options "-Xmx8G" HaplotypeCaller -R $reference \
-L chrY \
-L chrY_KI270740v1_random \
-O $1.b38.g.vcf.gz \
-I $1.recal.bam \
-ERC GVCF

All scripts are available at GitHubGist for convenience and collaboration.

Conclusion

This post outlines the base of haplogroup-r.org‘s new variant discovery workflow using the GATK toolkit.  We use BBmap to ensure empty or misnamed reads from the vendor do not adversely impact the alignment or variant calling processes, since GATK doesn’t have equivalent tools well documented at this current time.

Joint genotyping of the gVCFs is not discussed at this time as the GenomicsDB alterations are still in testing.  In the mean time you can reference the previous version of this workflow to use batches of gVCFs.

Prior versions of this workflow were able to create a callable loci report.  GATK4 lacks the walker tool that generated the BED file.  If this is a desired artifact in your analysis, it is recommended to run the GATK3 tool on the recalibrated BAM generated in Script 3.  The GATK team has stated a new version of CallableLoci is in the works, so the Gist for Script 3 will be updated at that time.

Future articles will look at evaluating alternative pipelines for alignment and variant calling.  Aligners such as Arioc and BarraCUDA offer optimizations to utilize GPU acceleration in the Smith-Waterman and Burrows-Wheeler algorithms.  Freebayes offers a Bayesian variant detection implementation written in C++.  The model eliminates GATK’s dependence on BaseRecalibrator preprocessing.

Updating a Y-DNA Variant Discovery Workflow

Last fall I shared the Y-DNA Variant Discovery Workflow – Pt. 1, which automated most of the process to apply for a consistent variant detection process in NGS data.  There have since been some new versions of the underlying software and changes in Broad Institute’s best practices.  This post will share the changes and discuss some of the why the change was necessary.

The scripts will take you into GRCh38 this time around, if you follow the steps exactly.  Those who wish to remain on GRCh37 for more convenient comparison with existing data sets should use the references in the older version of this post.

Getting the Tools

GATK 3.6 – Download page

Picard 2.6.0 – Download jar

Samtools – Download page

BWA – Github project

BBMap – Project page

Reference files available from 1000 Genomes.  The files are located in technical/reference/GRCh38_reference_genome in their FTP repository.  The known dbSNP and Mills_and_1000G indels files are in the other_mapping_resources subdirectory.

Setting Up Your Environment

The scripts are designed to work in a BASH shell on a *nix-like environment.  A Java 8 Run-time environment is required for many of the tools.  The setup is designed to have a Genomics directory in the user’s home directory.  This contains the unzipped version of the tools and a Reference directory containing the 1000 Genomes reference file.

Combining the Steps

gatk=~/Genomics/GenomeAnalysisTK-3.6/GenomeAnalysisTK.jar
picard=~/Genomics/picard-tools-2.6.0/picard.jar
reference=~/Genomics/Reference/GRCh38/GRCh38_full_analysis_set_plus_decoy_hla.fa
known=~/Genomics/Reference/GRCh38/Mills_and_1000G_gold_standard.indels.b38.primary_assembly.vcf
snpdb=~/Genomics/Reference/GRCh38/ALL_20141222.dbSNP142_human_GRCh38.snps.vcf
prefix_hla_hit=~/Genomics/Reference/GRCh38/resource-GRCh38/hs38DH-extra.fa
reference_fasta_file=~/Genomics/Reference/GRCh38/resource-GRCh38/hs38DH.fa.alt

# Parse out the actual Kit# from the local naming convention [kit#_surname_origin_testtype.bam]
read KIT1 KIT2 KIT3 <<<$(IFS="_"; echo $1)

# Shuffle the source BAM so that the alignment process is not colored by the source
samtools bamshuf -Ou $2 - | samtools bam2fq - > interleaved_reads.fq

# Big Y BAMs released for most of 2016 have had the read-pairs named incorrectly.
# BBMap includes a tool which will correct the problems.  We have to specify an
# unpaired file as well, since the post filtering removes a small percentage of
# pairs that were aligned to a different chromosome.
~/Genomics/bbmap/repair.sh -Xmx8g in=interleaved_reads.fq out=fixed.fq outsingle=single.fq

# Align the FASTQ using b38.  The pipe to bwa-postalt.js is not needed if your
# reference does not contain alt-contigs.
~/Genomics/bwa/bwa mem -t 4 -M -p $reference fixed.fq | ~/Applications/k8-0.2.1/k8-darwin ~/Genomics/bwa/bwakit/bwa-postalt.js -p $prefix_hla_hit $reference_fasta_file > aligned_reads.sam
 
# Order the resulting SAM to a BAM.
java -jar $picard SortSam INPUT=aligned_reads.sam OUTPUT=sorted_reads.bam SORT_ORDER=coordinate
 
# Mark the duplicate reads
java -jar $picard MarkDuplicates INPUT=sorted_reads.bam OUTPUT=dedup_reads.bam METRICS_FILE=metrics.txt
 
# Create read groups using dummy values, we use the sample id# from the naming convention as this will become visible in the g.vcf
java -jar $picard AddOrReplaceReadGroups INPUT=dedup_reads.bam OUTPUT=$1.b38.bam RGLB=lib1 RGPL=illumina RGPU=unit1 RGSM=$KIT1
 
# Index the resulting b38 aligned BAM
samtools index $1.b38.bam

# Sequencers have systematic technical errors which impact the quality scores.  This
# process applies machine learning to empirically model the errors.
java -Xmx4g \
 -jar $gatk \
 -T BaseRecalibrator -nt 1 -l INFO -cov ReadGroupCovariate -cov QualityScoreCovariate \
 -cov CycleCovariate -cov ContextCovariate -R $reference -o recal_data.table -I $1.b38.bam \
 -knownSites $known

# Now that the error model has been learned adjust quality scores.  Most will be
# reduced but some will increase.  None of the scores will be changed so much that
# the calling step will fail.
java -Xmx4g \
 -jar $gatk -T PrintReads -l INFO -R $reference -o $1.recalibrated.b38.bam \
 -I $1.b38.bam -BQSR recal_data.table --disable_bam_indexing

samtools index $1.recalibrated.b38.bam

# Call the realigned b38 BAM to gVCF for stage 2.
# Note: Remove the -L if variants other than chrY are desired
java -Xmx8g -jar $gatk -R $reference -T HaplotypeCaller \
-L chrY \
 -L chrY_KI270740v1_random \
 --dbsnp $snpdb \
 -o $1.b38.g.vcf \
 -I $1.recalibrated.b38.bam \
 --emitRefConfidence GVCF \
 -stand_call_conf 30 \
 -stand_emit_conf 0 -mmq 0

# Compute some statistics about how many bases can be called
java -jar $gatk \
 -T CallableLoci \
 -L chrY \
 -L chrY_KI270740v1_random \
 -R $reference \
 -I $1.recalibrated.b38.bam \
 -summary table.txt \
 -o callable_status.bed
 
# Save some disc space by compressing the VCF and creating an index
bgzip $1.b38.g.vcf
tabix -p vcf $1.b38.g.vcf.gz
rm $1.b38.g.vcf.idx
 
# POST PROCESS - REMOVE THE TEMP FILES
rm *.fq
rm interleaved_reads.fq.gz
rm aligned_reads.sam
rm sorted_reads.bam
rm dedup_reads.bam
rm $1.b38.bam
rm $1.b38.bam.bai

Batching the Files

Once your gVCFs have been processed it is advantageous to combine the files into batches of 100 to 200 files.  This is accomplished with the CombineGVCFs tool.

gatk=~/Genomics/GenomeAnalysisTK-3.6/GenomeAnalysisTK.jar
picard=~/Genomics/picard-tools-1.135/picard.jar
reference=~/Genomics/Reference/GRCh38/GRCh38_full_analysis_set_plus_decoy_hla.fa
known=~/Genomics/Reference/GRCh38/Mills_and_1000G_gold_standard.indels.b38.primary_assembly.vcf
snpdb=~/Genomics/Reference/GRCh38/ALL_20141222.dbSNP142_human_GRCh38.snps.vcf

target=Complete.b38.g.vcf

java -Xmx4g \
 -jar $gatk \
 -R $reference \
 -T CombineGVCFs \
 -o ~/Genomics/work/$target \
-V ./A.b38.g.vcf.gz \
-V ./Q.b38.g.vcf.gz \
-V ./R.b38.g.vcf.gz

bgzip ~/Genomics/work/$target
tabix -p vcf ~/Genomics/work/$target.gz
rm ~/Genomics/work/$target.idx

Discovering the Variants

Once the batches have been created, it’s time to genotype the set with GenotypeGVCFs.  This process allows every test to be sampled when a variant is discovered in anyone of the tests.  This is useful for generating a call matrix to feed into a parsimony tree algorithm later.

gatk=~/Genomics/GenomeAnalysisTK-3.6/GenomeAnalysisTK.jar
reference=~/Genomics/Reference/GRCh38/GRCh38_full_analysis_set_plus_decoy_hla.fa

java -Xmx4g \
-jar $gatk \
-T GenotypeGVCFs \
-R $reference \
-V $1 \
-o gatksnps.b38.vcf -nt 8

java -Xmx4g \
-jar $gatk \
-T VariantFiltration \
-R $reference \
--filterExpression "DP < 5 " --filterName "LowCoverage" \
--filterExpression "QD < 1.5 " --filterName "LowQD" \
--filterExpression "QUAL < 30.0 && QUAL < 50.0 " --filterName "LowQual" \
--filterExpression "QUAL <= 30.0" --filterName "VeryLowQual" \
--filterExpression "QUAL < 10.0" --filterName "REJECTED" \
--clusterWindowSize 10 \
-o gatksnps.b38.filtered.vcf \
-V gatksnps.b38.vcf

In Closing

The workflow above forms the central processing performed to produce the  CSVs on Haplogroup-R’s matrix.  While not fully tuned the results are comparable to similar comparisons available.  Processing time depends on the type of NGS test you are aligning and calling depends on the total number of samples.

Should you see glaring errors or minor omissions feel free to contact me.

A Simple Procedure for Combining FTDNA Variant Compare Files

One of the common themes among project administrators and analysts is a desire to easily compare multiple VCFs in a spreadsheet format.  Many have created custom scripts to perform this functionality.  The concept is simple read the tab-separated files into a collection and group the reads for each VCF by location.  Not everyone has the time to maintain or ability to create these utilities.

Combining the VCFs

The GATK suite of utilities mentioned several times before in this blog also comes with a tool for that.  CombineVariants allows a researcher to align the calls in multiple VCF files to begin creating a matrix.  The kits are arranged after the ninth column in the familiar fashion.  Reference reads yield a 0 genotype.  Alternate alleles yield a 0/1 or 1 genotype when actually covered.  Reads omitted from one of the samples simply use a ‘.’ or ‘./.’ symbol.

Unfortunately, FTDNA is rather incomplete in how the the VCF header is populated.  The REJECTED filter is quite common, but the criteria are not qualified.  Defining the meaning of the filter values is a common practice to allow analysts to make educated decisions around including these calls or not.  By default GATK considers this to be an error in FTDNA’s application of the file standard.  Fortunately, there is an –unsafe argument that can be used to ignore the error.

java -jar ~/Genomics/GenomeAnalysisTK-3.6/GenomeAnalysisTK.jar \
-T CombineVariants \
-R ~/Genomics/Reference/hg19/hg19.fasta \
-V Sample1_BigY_RawData_20160117/variants.vcf \
-V Sample2_BigY_RawData_20160712/variants.vcf \
-o output.vcf \
-genotypeMergeOptions UNIQUIFY \
--unsafe

Simplifying the Table

Once all of the VCFs have been assembled in this manner it’s a simple matter of using VariantsToTable.  This manipulation tool allows you to specify which components of the VCF to write to a table.  Consult the documentation for all the options, but a typical run would look similar to the following.

java -jar ~/Genomics/GenomeAnalysisTK-3.6/GenomeAnalysisTK.jar \
-T VariantsToTable \
-R ~/Genomics/Reference/hg19/hg19.fasta \
-V output.vcf \
-F CHROM -F POS -F REF -F ALT \
-GF GT \
-o out.txt

The command produces a table with the chromosome, start location, reference allele, alternate alleles, and the allele present in each of the samples.

CHROM   POS     REF     ALT     aa827c54-4bfa-4a65-9d42-f74002917ac5.variant.GT cc85f1c9-ed2e-4ea1-8b35-698e9517d7a7.variant2.GT
chrY    2650265 A       .       A       A
chrY    2650371 G       .       G       G
chrY    2650701 G       .       G       G
chrY    2650853 G       .       G       G

The unfortunate effect is the samples are named with their test unique identifier and not their more natural kit #.  These could be overcome using Excel’s vlookup() or similar function from a real database.

Going Further

There are a good number of other tools included, which can improve your quality of life.  If you notice that INDELs are inconsistently formatted between tests, LeftAlignAndTrimVariants will provide the minimal representation.  Named variants can be annotated using VariantAnnotator.  This tool typically uses NCBI’s dbSNP as the known variant input.  That does not stop one from using an export of ybrowse.org‘s GFF export and creating their own reference VCF.  There are further tools for filtering and selecting that one can use as well.

Big Y BAMs have returned!

FTDNA temporarily suspended issuing BAMs for Big Y in early 2016 to rework their hosting solution.  The good news is that the feature was restored in early April 2016.  Customers who have been patiently waiting to gain access now have the ability to retrieve the data for third party analysis by companies like YFULL and FGC.  The workflow previously shared on this site can also be used to generate a VCF with higher sensitivity settings than those provided by FTDNA.

But something is missing!

As noted in the correction to the previous workflow post, BWA MEM benefits from using paired-end read markers.  Enabling this functionality results in better alignment scoring.  The new BAMs have the reads named in a way in which BWA is unable to match the pairs.  This results in realignment being done as if all reads are single-end sequences.  There have also been reports of this causing issues with Ugene, a popular tool for viewing alignments graphically.

BBtools to correct the paired-read names

BBmap contains a large number of utility scripts designed to normalize and fix data issues.  One such tool is repair.sh.  Designed to repair reads that are disordered or have missing mates, the Java-based shell application includes the functionality needed to rename the reads as expected by BWA MEM.  Those using a similar version of the variant discovery workflow can introduce a new step prior to re-alignment quite simply.
bbmap/repair.sh in=interleaved_reads.fq out=fixed.fq outsingle=single.fq

The fixed.fq file contains the Big Y paired-read names needed to allow aligner packages to use the paired-end modes.

With my recent decision to migrate my Y DNA analysis efforts to use GRCh38 a lurking problem was hidden under the surface.  Over the last few years Big Y has been an amazing success for the R haplogroup.  There are 1954 results collected using reference build 19 currently sitting on my hard disc.  With only 18% of the BAMs available to realign a migration plan is needed.

Fortunately University of California, Santa Cruz provides a set of mapping files that allows results to be converted to another reference.  The hg19 to hg38 chain is found in their repository.

Using a program like CrossMap and the UCSC chain files, you can fairly easily Convert Big Y VCF/BEDs to GRCh38.  CrossMap requires a copy of the target reference sequence for the VCF.  For compatibility with the 1000 Genomes dataset I am using theirs.

Once everything is installed converting the files is straightforward:

python #{crossMapPath} bed #{chainPath} regions.bed regions.hg38.bed

python #{crossMapPath} vcf #{chainPath} variants.vcf #{referencePath} variants.hg38.vcf

There is a gotcha in the resulting VCF’s.  Roughly one in ten files exhibits the reference allele listed in the alternate alleles in positions above hg38 chrY 56700000 like so:

chrY 56707268 . T T,C 74.7188 REJECTED . GT 1/2

This makes the context invalid in htsjdk, which is the library used in my database load process.  I believe the root cause is an incorrect ancestral allele in the hg19 VCFs.  My fix is to simply remove the reference allele from the alternate list and deduct one from the genotype values using a short Ruby script prior to import.