Tutorial for human WGS

About this page

Here, we will analyze two human WGS from the 1000GP as a test. We will use Singularity 3 to use MEGAnE in this tutorial, but MEGAnE also works on Docker.

Install Singularity 3

Please follow this official instruction to install Singularity 3.

Install MEGAnE

See the Installation page for more details.

sudo singularity build MEGAnE.sif docker://shoheikojima/megane:tagname
# or
singularity build --fakeroot MEGAnE.sif docker://shoheikojima/megane:tagname

Download test files

We are going to use two 30x WGS, NA12878 and NA18999, for test. The CRAM files that are available from the 1000GP are position-sorted, so these can be directly analyzed by MEGAnE. We will download the reference genome (GRCh38DH) for the CRAM file as well.

wget ftp://ftp.sra.ebi.ac.uk/vol1/run/ERR323/ERR3239334/NA12878.final.cram
wget ftp://ftp.sra.ebi.ac.uk/vol1/run/ERR323/ERR3239577/NA18999.final.cram
wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa

Prepare MEGAnE k-mer set

First, we will process the reference genome (i.e. GRCh38DH). This requires ~10 min and 50GB RAM usage.

• -fa: the path to the reference fasta file.
• -outdir: any name of a directory where the output files will be stored.

singularity exec MEGAnE.sif build_kmerset \
-fa ./GRCh38_full_analysis_set_plus_decoy_hla.fa \
-outdir megane_kmer_set

The command above should generate two files: ./megane_kmer_set/GRCh38_full_analysis_set_plus_decoy_hla.fa.mk and ./megane_kmer_set/GRCh38_full_analysis_set_plus_decoy_hla.fa.mi. Those two files are the MEGAnE k-mer set.

Analyze WGS

This is the main step of MEGAnE, i.e. MEV calling and genotyping. This requires 1-2 hours and up to 15GB RAM per WGS.

• Here, we use call_genotype_38 command, because the WGS is mapped on the GRCh38-related genome build (i.e. GRCh38DH).
• If your sample is mapped on a GRCh37-related genome build, please use call_genotype_37 command.
• -i: the path to the input CRAM file (can take BAM format as well).
• -fa: the path to the reference fasta file.
• -mk: the name of the .mk file that was generated above.
• -outdir: any name of a directory where the output files will be stored.
• -sample_name specifies any name. This name will be used in the output VCF file. Please avoid using wild cards (e.g. . / ' " | \ etc.).
• -p: Thread number. We recommend 2 to 4 threads.

singularity exec MEGAnE.sif call_genotype_38 \
-i ./NA12878.final.cram \
-fa ./GRCh38_full_analysis_set_plus_decoy_hla.fa \
-mk ./megane_kmer_set/GRCh38_full_analysis_set_plus_decoy_hla.fa.mk \
-outdir MEGAnE_result_NA12878 \
-sample_name NA12878 \
-p 2

singularity exec MEGAnE.sif call_genotype_38 \
-i ./NA18999.final.cram \
-fa ./GRCh38_full_analysis_set_plus_decoy_hla.fa \
-mk ./megane_kmer_set/GRCh38_full_analysis_set_plus_decoy_hla.fa.mk \
-outdir MEGAnE_result_NA18999 \
-sample_name NA18999 \
-p 2

After the run, you can check whether the analysis finished correctly by seeing the log file. If there is a message All analysis finished! Thank you for using MEGAnE! at the end of the log file, it means the analysis finished correctly.

cat ./MEGAnE_result_NA12878/for_debug.log
cat ./MEGAnE_result_NA18999/for_debug.log

Joint calling

Next, we will merge the VCF from NA12878 and NA18999. In this step, we need to run twice, one is for ME insertions, and one is for ME absences. Here, we only merges two samples, so it will finish in a few minutes, but would take longer than several hours and >10GB RAM when merging >1000 samples.

• -merge_mei: when this is specified, MEGAnE generates a joint callset of non-reference ME insertions.
• -merge_absent_me: when this is specified, MEGAnE generates a joint callset of reference ME variations.
• -f: the path to the file containing the paths to the individual VCF files (i.e. dirlist.txt).
• -fa: the path to the reference fasta file.
• -outdir: any name of a directory where the output files will be stored.
• -cohort_name: any name. This name will be used as ME IDs in the output VCF file. Please avoid using wild cards (e.g. . / ' " | \ etc.).

# first, list up the output directories from above in a file
ls -d $(pwd)/MEGAnE_result_* > dirlist.txt

# merge non-reference ME insertions
singularity exec MEGAnE.sif joint_calling_hs \
-merge_mei \
-f dirlist.txt \
-fa ./GRCh38_full_analysis_set_plus_decoy_hla.fa \
-outdir joint_calling_result \
-cohort_name test

# merge reference ME variations
singularity exec MEGAnE.sif joint_calling_hs \
-merge_absent_me \
-f dirlist.txt \
-fa ./GRCh38_full_analysis_set_plus_decoy_hla.fa \
-outdir joint_calling_result \
-cohort_name test

See the joint callset

In the directory joint_calling_result, there should be two files: test_MEI_jointcall.vcf.gz and test_MEA_jointcall.vcf.gz. The file with MEI contains the non-reference ME insertions, while that with MEA (which means ME absence) contains the reference ME variations.