mamba install -c bioconda methphaser
, conda install -c bioconda methphaser
methphasing
into ~/.bashrc
for meth_phaser_parallel
to call it. exact command: export "PATH=/path/to/methphaser/:$PATH"
methphaser.yaml
!!! NOTE: MethPhaser won't work on sex chromosomes because of random X inactivation! Please only consider autosomes. !!!
!!! We recommend using MethPhaser for at least 30X coverage ONT reads since most of the SNV phasing methods (HapCUT2, Whatshap) only work reasonably above that coverage !!!
Test data at https://zenodo.org/record/8360845
Step 1: Run MethPhaser to get block relationship. MethPhaser default ignores the largest unphased region (-ml -1
), if the input data does not contain any large poorly mapped region like centromere, please use -ml -2
to not ignore any gap.
./meth_phaser_parallel -b test_data/HLA.R10.haplotagged.bam -r test_data/GCA_000001405.15_GRCh38_no_alt_analysis_set.chr6.fna -g test_data/LSK.filtered.gtf -vc test_data/HLA.R10.phased.vcf.gz -o test_data/work
Step 2: Run post processing script to get modified reads and vcf file
./meth_phaser_post_processing -ib test_data/HLA.R10.haplotagged.bam -if test_data/work/ -ov test_data/output.vcf -ob test_data/output -vc test_data/HLA.R10.phased.vcf.gz
usage: meth_phaser_parallel [-h] -b -r -g -vc [-vt] [-t] [-ml] [-c] [-a] [-o] [-k]
methphaser: phase reads based on methlytion informaiton
optional arguments:
-h, --help show this help message and exit
-vt , --vcf_truth Truth vcf file for benchmarking
-t , --threads threads
-ml , --max_len maximum homozygous region length for phasing, default: -1 (ignore the largest homozygous
region, centromere), input -2 for not skipping anything
-c , --cut_off the minimum percentage of vote to determine a read's haplotype
-a , --assignment_min
minimum assigned read number for ranksum test
-o , --output_dir output_directory
-k , --k_iterations use at most k iterations, use -1 for unlimited iterations.
Required arguments:
-b , --bam_file input methylation annotated bam file
-r , --reference reference genome
-g , --gtf gtf file from whatshap visualization
-vc , --vcf_called called vcf file from HapCUT2
Note: The BAM file should only contain primary alignment (no secondary alignment and supplementary alignment), otherwise it will trigger pysam MM tag reading error. Suggested command:
samtools view -bF 2304 -o output.bam input.bam
Recomanded phasing flow: