-
phased multiple alleles calling for Next Generation sequencing generated long sequencing reads (designed for), e.g. PacBio HiFi
-
call variants of STR/SSR, SNP, InDel at the same time
-
run faster with parallel computing
-
also works for short reads
-
Programmed and developed by Ph.D Xuewen Wang for those love genomics and science
Macrohaplotype caller for STR, SNP, and InDel from NGS long-read sequencing data, especially for PacBio HiFi reads The MacroHapCaller calls targeted STRs/SSRs, SNPs, and InDels sites simultaneously in a row from each NGS read and clusters the variants into phased haplotype strings. MacroHapCaller is the best tool to analyze all haplotypic variants from genetically inherited DNA. It suits diploid, polyploid, and complex DNA mixtures from many individuals, e.g. DNA forensics. MacroHapCaller is programmed in Java with parallel computing enabled so it can run on any computing platform.
The MacroHapCaller is also integrated into a pipeline with a shell script or Python3. The pipeline takes the BAM or fastq reads as input, align reads to the genome, sort, and index alignment, and then user can use the output to run MacroHapCaller. Additional tools are also available as util tools. For util usage : java -jar appName.jar
where appName should be replaced with a detailed tool name.
MacroHapCaller runs fast. It takes ~ 2 mins for merged data from four full runs of PacBio SMRTcell HiFi reads on a Linux machine with 12 threads, and a maximum allowed 120G RAM.
V1.0
V0.3, V0.4
the newest config for the panels: config_v0.4
updates: optimized the CSF1PO locus, removed variant sites (2 in the snp panel) in the noisy regions at ch5:150,083,350-150,083,403; chr5:150,084,343-150,084,480; corrected chr5 150076781 to CSF1PO
The new config is further optimized to have better robust MH alleles based on our real data testing with more samples.
The Java standard runtime environment (SE) is needed, which is installed in most computer. You may just need to update it to the latest version. The Latest LTS Java or Java SE 17 or higher is recommended.
Step 1. Go to the Orcale website, download and unzip the latest SE file from https://download.oracle.com/java/21/latest/jdk-21_linux-x64_bin.tar.gz ;
This binary version does not need installation, just unzip it.
Step 2. Add the path to your unzipped Java bin directory before running the command. e.g., the Java is unzipped into c/java21, then the path will be c/java21/jdk-21/bin, the running command will be :
path=c/java21/jdk-21/bin
$path/java -jar MacroHapCaller##.jar
Or you can install the Java JDK 21 in your computer, then you can run without $path/ before java , just like instructed below.
Or add the path into to environment.
The MacroHapCaller is programmed in Java and compiled. So just download the APP and run it.
git clone https://github.com/XuewenWangUGA/MacroHapCaller
cd MacroHapCaller
Or just download all folders from Github and uncomppress the zip file. then go to the unzipped folder.
cd MacroHapCaller
fastq-filter is the only needed dependency if want to filter the reads. Details on https://github.com/LUMC/fastq-filter. To install
pip install fastq-filter
Open a terminal window. Type the following command for help:
java -jar MacroHapCaller##.jar
where ## is the version number, e.g. : 0.4.
to change the version number as needed.
usage: java -jar -Xmx100G MacroHapCaller1.0.jar [options]
-a,--strAnchorFile configure file of anchors for STRs in tabular plain text
-c,--MinReadCount integer, minimum read count to report a macrohaplotype, default [2]
-d,--inDelBedFile configure file of InDels in BED format in tabular plain text
-i,--input input BAM file with a path
-l,--homopolymerLen integer, homopolymer length threshold to skip call variant,default [4]
-m,--MaxMismatches integer, maximum distance of edit distance in anchor match, default [1]
-n,--snpPanelPosFile configure file of SNPs in tabular plain text
-o,--output output file name for macrohaplotype result
-p,--MinReadProportion double, minimum read proportion of all reads for each locus to report a macrohaplotype, default [0.001]
-q,--LowQualThreshold integer, phred-scale quality threshold for variant bases, default [15]
-r,--refGenomeFastaFile Genome reference sequence file in .fasta format
-s,--span T|F, T to report the read spaning information, F not to report, default T. For v1.0 or above only
-t,--ThreadNumber integer, the number of computing threads, default [12]
Download the subfolders and put them inside the folder of "MacroHapCaller". Then users can use parameters and configure files as used in the demo listed below.
java -jar -Xmx120G MacroHapCaller0.4.jar -i demo_data/hg002.8kampl.Q30.4kreads.fastq.gz_GRCh38.bam -o hg002.8kampl.Q30.fastq.gz_GRCh38.HapVar.tsv -a config_v0.3/CODISSTR_anchor.XW.config_v0.3.txt -d config_v0.3/MHindels_v0.3.bed -n config_v0.3/MHsnps.pos_v0.3.txt -r GCA_000001405.15_GRCh38_no_alt_analysis_set.fasta -l 2 -m 1 -q 15 -c 100 -p 0.01 -t 12
the option -Xmx120G is to set allowed maximum memory when running this app. it can be ignored or changed to a small memory , e.g., -Xmx4G
E.g., human genome hg38:
Download the genome sequence of human from the 1000 Genome Project to the folder "VarSeqStitcher" ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa
wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa
rename the genome sequence file as a short name:
mv GRCh38_full_analysis_set_plus_decoy_hla.fa genome.fasta
Then index the genome sequence with samtools (tool link: https://www.htslib.org/) samtools faidx genome.fasta
This pipeline will preprocess the original PacBio HiFi reads in BAM to fastq, read quality control, umi-analysis (optional), statistical summary, read-reference alignment, sort, and index. The major results are a read-reference alignment in the BAM format and BAM index, as well as statistical report files. Edit the path and files name to your specified names.
MH_umi_dedup_map_pipe.sh
Assume the alignment output will be named as "alignment.bam"
java -jar -Xmx100G MacroHapCaller0.4.jar -i alignment.bam -o out.tsv -a config_v0.3/CODISSTR_anchor.XW.config_v0.3.txt -d config_v0.3/MHindels_v0.3.bed -n config_v0.3/MHsnps.pos_v0.3.txt -r genome.fasta -l 2 -m 1 -q 15 -c 100 -p 0.01 -t 12
Replace the real file name with your file name, and add file location path before the file if necessary. The above command uses the configure files coming with the software for human CODIS STR based HG38.
The output will be in file: out.tsv, which can be viewed in spreadsheet, Mircosoft Excel or any other text file editor.
-
a read-reference alignment file in BAM format. This alignment file could be generated from
MH_umi_dedup_map_pipe.sh
or user-generated alignment file in .Bam format -
Configure files of targeted variant positional information for each of STRs, SNPs, and InDels
The configure file listed the STR loci information in a tabular separated text file. An example file for forensic CODIS STR of human genome is given in "CODISSTR_anchor.XW.config_v0.2.txt" Each line for one STR locus. There is one fixed headline before starting the detailed STR locus. Every line has 11 columns. The integer value for the last column is the estimated maximum length of TR, where 500 is enough for most cases. The coordinate starts from 1.
Chrom ChromStartPos_Str ChromEndPos_Str Name Repeat_unit_length Motif_must_present Inner_offset Anchor_left Anchor_right LeftAnchorStartPos RightAnchorEndPos MaxRefLength
chr1 230765214 230765259 D1S1656@STR_118267 2 CACACA 0 CAGAAAATGAGAACACATG GGTTATGCCAAAAGGGC 230765190 230765281 500
chr1 230769616 230769683 D1S1656 4 TCTATCTATCTA 0 TCAGAGAAATAGAATCACTA TGAGCAACACAGGCTTGAAC 230769556 230769721 500
chr2 1489653 1489684 TPOX 4 AATGAATG 0 CAGAACAGGCACTTAGGGAA AACGCTGACAAGGACAGAAG 1489624 1489716 500
The config file for InDel is in a BED format tabular separated plain text file, with coordinates started from 0. The REF and ALT are the bases for InDel in the reference genome and variant. An example file is given for 8-kb fragments containing human CODIS STR. For Details, please read our paper. There is one fixed headline and 5 columns in each line. Then each subsequent line is for one InDel site. If multiple InDel for the same coordinate, list one of them is enough.
Chrom ChromStart ChromEnd REF ALT
chr1 230764747 230764748 GA G
chr1 230764749 230764750 T TC
chr2 1485704 1485705 C CT
chr2 1486795 1486803 TAGTGGGG T
chr2 1487522 1487526 CAGG C
chr2 1489924 1489975 GCACACAGGAGGAGTCACGACAGAGCAGTGTAAGAGCCGCCACGTGGGTCC G
The config file for SNP is in a tabular separated plain text file, with coordinates started from 1. There is one fixed headline and 5 columns in each line. Then each subsequent line is for one SNP site. An example file is given for 8-kb fragments containing 20 human CODIS STRs. the "ID" column is the name of a DNA fragment, which uses the same name as one of the STR in this 8-kb fragment as STR configure. All SNPs with the same "ID" will be put together with STR alleles in the same fragment once the ID of STR is shared between SNP and STR configure file.
#CHROM POS ID REF ALT
chr1 230765010 D1S1656 A G
chr1 230765280 D1S1656 G A
chr1 230765336 D1S1656 C T
chr1 230765352 D1S1656 G C
Macrohaplotypes and supported read count in tabular text file .tsv. E.g., the two macrohaplotypes around FBI's CODIS loci D2S441 in 8 kb PacBio HiFi reads of benchmark reference hg002. Locus details are at FBI
The first line shows the total reads for this locus D3S1358. The position lines present the coordinate of each targeted variant site on chromosomes in the human genome Hg38. Then the macrohaplotypes, including three parts of SNPs, InDels and STRs separated by ";". If there is no specififed variant site for each type, the report will use "null" to indicate this type variant is not specified.
#Total hapVar: D3S1358 47621
#Markername Counts HapVarLen hapVar(s)
# Position:45534683,45534822,45535067,45535094,45535770,45535776,45535998,45536129,45536613,45536683,45536800,45536841,45537042,45537158,45537183,45537457,45537777,45537899,45538041,45538207,45538590,45538748,45539142,45539154,45539315,45539413,45539452,45539580,45539887,45539923,45540056,45540467,45540750,45540754,45540758,45541155,45541805,45541806,45537166,45537617,45539736,45540709
D3S1358 13843 108 T,A,T,C,A,C,T,A,G,A,T,A,A,G,G,C,C,C,T,T,C,G,T,G,A,C,G,G,T,A,G,C,A,A,A,T,C,G;TG,CA,TA;TCTATCTGTCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTA
D3S1358 13600 101 G,G,T,T,A,A,T,G,C,G,C,A,A,T,G,T,G,T,T,T,C,A,T,G,A,T,A,A,T,A,G,C,G,G,A,T,C,G;T,C,T;TCTATCTGTCTGTCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTATCTA
The macrohaplotypes can be visualized or edited in a text processor, Microsoft excels, and spreadsheet.
The graphically view and compare the macrohaplotypes, please use our graphic tool: USAT
A colorful displaying the macrohaplotypes in a web browser is also developed. For this utility, please email me for a free copy.
These tools are designed to process bam in any computing platform, bam processing for windows. It can be a Windows version of some functions of samtools
. Please make sure you have the latest version of JAVA standard environment installed (Java 17 or higher version). This tool has been tested on Windows, Linux, and MacOS.
Java environment download link: https://www.oracle.com/java/technologies/downloads/.
for help: java -jar BamSubset.jar
Usage:
java -jar BamSubset.jar Bed_file Input.bam reduced_input.bam
for help:
java -jar BamAlignSortIndex.jar
Usage:
java -jar BamAlignSortIndex.jar input.bam
This work was sponsored in part by award 15PNIJ-21-GG-04159-RESS, awarded by the National Institute of Justice, Office of Justice Programs, U.S. Department of Justice.
Xuewen Wang, Melissa Muenzler, Jonathan King,Muyi Liu, Hongmin Li, Bruce Budowle, Jianye Ge. 2024. A complete pipeline enables haplotyping and phasing macrohaplotype in long sequencing reads for polyploidy samples and a multi-source DNA mixture. https://doi.org/10.1002/elps.202300143
or https://analyticalsciencejournals.onlinelibrary.wiley.com/doi/10.1002/elps.202300143