Run for human WGS

Usage for human WGS (GRCh37 and GRCh38-related genomes)

Step 0. Prepare MEGAnE k-mer file

• Before analyzing your BAM/CRAM files, you need to make MEGAnE k-mer files from your human reference genome (e.g. GRCh38DH, hs37d5, etc).
• Usually, this takes ~10 min and requires ~50GB RAM.

sif=/path/to/MEGAnE.sif

singularity exec ${sif} build_kmerset \
-fa /path/to/reference_human_genome.fa \
-prefix reference_human_genome \
-outdir megane_kmer_set

Step 1. Call and genotype polymorphic MEs

• MEGAnE can analyze both BAM and CRAM format.
• MEGAnE supports BAM/CRAM mapping to both GRCh37- and GRCh38-related genomes.
• One 30x human WGS takes ~1 hour using 4 threads.
• MEGAnE currently supports BAM/CRAM files generated by BWA MEM, BWA MEM2, and DRAGEN. In the case of DRAGEN, please also specify the -skip_unmapped option.
• You may want to use further options (e.g. -low_dep option and -skip_unmapped option). Please also read "Things needed to be checked before analysis" page and "Options for MEGAnE step 2" page for details.

In the case of BAM/CRAM file mapping to GRCh37-related genome

• Please use the call_genotype_37 command.

sif=/path/to/MEGAnE.sif

singularity exec ${sif} call_genotype_37 \
-i /path/to/input.bam \
-fa /path/to/reference_human_genome.fa \
-mk /path/to/megane_kmer_set/reference_human_genome.mk \
-outdir MEGAnE_result_test \
-sample_name test_sample \
-p 4

In the case of BAM/CRAM file mapping to GRCh38-related genome

• Please use the call_genotype_38 command.

singularity exec ${sif} call_genotype_38 \
-i /path/to/input.cram \
-fa /path/to/reference_human_genome.fa \
-mk /path/to/megane_kmer_set/reference_human_genome.mk \
-outdir MEGAnE_result_test \
-sample_name test_sample \
-p 4

Step 2. Joint calling

• After the analysis of multiple BAM/CRAM files, you can make a joint call.
• This will take several hours when merging 1000s samples.
• The discovery of MEVs depends on sequencing platform (e.g. read depth, read length, mapping software, etc). We do NOT recommend making a joint callset that includes WGS datasets sequenced by different platforms, since some MEVs may be genotyped incorrectly across samples. If you choose to make such joint callset, we strongly recommend intense QC after joint calling.

sif=/path/to/MEGAnE.sif

# first, list up samples (output directories from Step 1) you are going to merge
ls -d /path/to/[all_output_directories] > dirlist.txt

# merge non-reference ME insertions
singularity exec ${sif} joint_calling_hs \
-merge_mei \
-f dirlist.txt \
-fa /path/to/reference_human_genome.fa \
-cohort_name test

# merge reference ME polymorphisms
singularity exec ${sif} joint_calling_hs \
-merge_absent_me \
-f dirlist.txt \
-fa /path/to/reference_human_genome.fa \
-cohort_name test

• MEGAnE supports joint calling from massive WGS (e.g. more than 10,000 WGS).
• In our computational environment, joint calling of 10,000 WGS is not a problem, but when merging more than that number, it may require massive memory space (> 50GB). In such case, we recommend -chr option. This option only merges a specified chromosome, so it does not require a large memory. User need to run for all chromosomes (i.e. chr1 to chr22, chrX, chrY). Please also see "Options for MEGAnE step 2" page as well.

# merge non-reference ME insertions on chr1
singularity exec ${sif} joint_calling_hs \
-merge_mei \
-f dirlist.txt \
-fa /path/to/reference_human_genome.fa \
-chr chr1 \
-cohort_name test

(Optional) Step 3. Make a joint call for haplotype phasing

• The step 2 above will generate two VCF files. This step merges the two files generated in the step 2.
• When merging the two VCF files, MEGAnE removes multi-allelic variants.
• This is particularly useful when you will do haplotype phasing using MEGAnE's results.

sif=/path/to/MEGAnE.sif

singularity exec ${sif} reshape_vcf \
-i /path/to/jointcall_out/[cohort_name]_MEI_jointcall.vcf \
-a /path/to/jointcall_out/[cohort_name]_MEA_jointcall.vcf \
-cohort_name test

(Optional) Example haplotype phasing of MEGAnE's result

• Here is one example of how to phase MEGAnE's result.
• This step requires external software.

threads=6
memory=30000

# first, generate a SNP VCF
snp_vcf=/path/to/SNP.vcf chr=1
exclude=/path/to/vcf_for_phasing/cohort_name_biallelic.bed.gz
chr=1

plink2 \
--threads ${threads} \
--memory ${memory} \
--vcf ${snp_vcf} \
--make-pgen \
--max-alleles 2 \
--mac 2 \
--hwe 1e-6 \
--chr ${chr} \
--indiv-sort natural \
--exclude bed0 ${exclude} \
--export vcf bgz \
--out SNP


# next, generate a ME VCF
me_vcf=/path/to/vcf_for_phasing/[cohort_name]_biallelic.vcf.gz

plink2 \
--threads ${threads} \
--memory ${memory} \
--vcf ${me_vcf} \
--make-pgen \
--max-alleles 2 \
--mac 2 \
--hwe 1e-6 \
--chr ${chr} \
--indiv-sort natural \
--vcf-half-call missing \
--export vcf bgz \
--out ME


# merge SNP and ME VCFs
cat SNP.vcf.gz > _SNP_ME.vcf.gz
zcat ME.vcf.gz | grep -v '#' | pigz -c >> _SNP_ME.vcf.gz
zcat _SNP_ME.vcf.gz | pigz -c > SNP_ME.vcf.gz
rm  _SNP_ME.vcf.gz

plink2 \
--threads ${threads} \
--memory ${memory} \
--vcf SNP_ME.vcf.gz \
--make-pgen \
--sort-vars \
--vcf-half-call missing \
--out SNP_ME

plink2 \
--threads ${threads} \
--memory ${memory} \
--pfile SNP_ME \
--export vcf bgz \
--out SNP_ME_sorted

tabix SNP_ME_sorted.vcf.gz


# haplotype-phasing
map=/path/to/genetic_map

shapeit4 \
--input SNP_ME_sorted.vcf.gz \
--map ${map} \
--region ${chr} \
--thread ${threads} \
--output SNP_ME.phased.vcf.gz