GLIMPSE is a tool of software for imputation and phasing of low-coverage datasets in the form of genotype gikelihoods (GLs) at all variant positions. GLIMPSE requires HTSlib and BCFtools v1.7 (or later) as a requirement, since it makes use of indexed VCF/BCF files. Here we also used BCFtools to compute genotype likelihoods.

In this tutorial, we show how to run GLIMPSE from sequencing reads data (BAM/CRAM file) to obtain refined genotype and haplotype calls. A small dataset containing chromosome 22 downsampled 1x sequencing reads for one individual (NA12878) is provided with GLIMPSE, together with a reference genome FASTA file (for chromsome 22) needed for genotype likelihood computations.

A minimal set of data to run this pipeline can be found in the folder GLIMPSE/tutorial. All datasets used here are in GRCh38/hg38 genome assembly.

1.1. Setting the environment and binaries

All scripts have been written assuming  > cd GLIMPSE/tutorial/ as current directory. Therefore, we require to move to the tutorial directory:

 > cd GLIMPSE/tutorial/

From this folder we build all the GLIMPSE software needed (chunk, phase, ligate, sample and concordance). For this reason, we require that GLIMPSE can correctly compile on your machine and you might be required to edit the Makefiles manually. See installation instructions if there are any problems at this stage. We created a script to configure a directory containing symbolic links to the GLIMPSE binaries we will need later. To run the setup script, simply run:

 > ./step1_script_setup.sh

If the script runs correctly, the software compiles and the bin folder containing symbolic links to the binaries has been created, everything is correctly setup.

1.2. Low coverage reads

In this example we downsampled the publicly available 30x data available for NA12878 provided by the Genome In A Bottle consortium (GIAB). We downloaded the full dataset, kept only chromosome 22 reads and downsampled to 1x. The resulting BAM file can be found in:

GLIMPSE/tutorial/NA12878_1x_bam/NA12878.bam[.bai]

Details of how we downsampled the dataset are provided in appendix section A1.1.

In order to run accurate imputation we recommend to use a targeted reference with ancestry related to the target samples, if possible. For this tutorial we use the 1000 Genomes Project 30x NYGC reference panel. Since the reference panel contains data of our target sample and few relatives, we need to remove them from the reference panel. This steps downloads the 1000 genomes b38 data from the EBI ftp site.

2.1. Download of the reference panel files

The 1000 Genomes Project reference panel is publicly available at EBI 1000 genomes ftp site. We can download chromosome 22 data using:

 > wget -c http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20201028_3202_phased/CCDG_14151_B01_GRM_WGS_2020-08-05_chr22.filtered.shapeit2-duohmm-phased.vcf.gz{,.tbi}

The reference panel size is approximately 520 MB.
Update (20/05/2021): users reported that the download of the reference panel can be slow and it is often interrupted. Unfortunately, we do not have control on that. However, using wget with the -c option allows to restart the download from where it was interrupted. Additionally, we recommend to compare the checksum to what it was reported. Instructions for this can be found in the following script:

step2_script_reference_panel.sh

2.2. Remove NA12878 (and family) from the reference panel and perform basic QC

We used BCFtools to remove sample NA12878 and related samples from the reference panel and we export the dataset in BCF file format (for efficiency reasons). We also performs a basic QC step by keeping only SNPs and remove multiallelic records. The reason for this is multiple: GLIMPSE only handles bi-allelic SNPs and Bcftools does not perform a good job at calling indels, causing potential problems at neighboring SNPs, therefore we remove them completely from the analysis. Alternatively, it is possible to keep the indels in the reference panel and force missing-data imputation there.
The original reference panel files are then deleted from the main tutorial folder:

 > CHR=22
 > bcftools norm -m -any CCDG_14151_B01_GRM_WGS_2020-08-05_chr22.filtered.shapeit2-duohmm-phased.vcf.gz -Ou --threads 4 |
      bcftools view -m 2 -M 2 -v snps -s ^NA12878,NA12891,NA12892 --threads 4 -Ob -o reference_panel/1000GP.chr22.noNA12878.bcf
 > bcftools index -f reference_panel/1000GP.chr22.noNA12878.bcf --threads 4
 > rm CCDG_14151_B01_GRM_WGS_2020-08-05_chr22.filtered.shapeit2-duohmm-phased.vcf.gz*

The resulting reference panel files can be found in:

GLIMPSE/tutorial/reference_panel/1000GP.chr22.noNA12878.bcf[.csi]

The reference panel can now be used for genotype likelihoods calculations (only the sites, not the data) and imputation.

GLIMPSE requires input data to take the form of Genotype Likelihoods (GLs). GLs need to be computed at all target individuals and all variant sites present in the reference panel of haplotypes used for the imputation. In this section, we describe a simple procedure to compute GLs from sequencing data using BCFtools. However, other callers, such as GATK for instance, can be also considered as soon as the GLs are encoded in a VCF/BCF file using the FORMAT/PL field.

Please note that BCFtools does not call correctly genotype likelihoods for indels and this might affect the quality of the imputation even at non-indels. For this reason in this tutorial we only focus on bi-allelic SNPs.

For this tutorial we only follow step 3.1 and step 3.2, because step 3.3 is not necessary since we have only one target sample in our study panel. The code for sections 3.1 and 3.2 is included in the script step2_script_GL.sh.

3.1. Extracting variable positions in the reference panel

We a VCF/BCF file containing only sites to tell BCFtools at which positions to make a genotype call. Since BCFtools does not compute correctly genotype likelihood for indels, here we only focus on SNPs (however, GLIMPSE can impute any type of variants as soon it is bi-allelic and has GLs being properly defined). To perform the extraction from the chromosome 22 of a reference panel 1000GP.chr22.noNA12878.bcf, run first BCFtools as follows:

 > bcftools view -G -m 2 -M 2 -v snps reference_panel/1000GP.chr22.noNA12878.bcf -Oz -o reference_panel/1000GP.chr22.noNA12878.sites.vcf.gz
 > bcftools index -f reference_panel/1000GP.chr22.noNA12878.sites.vcf.gz

Then, convert the output file 1000GP.chr22.noNA12878.sites.vcf.gz into TSV format and index the file using tabix (requires htslib in PATH), using this command:

 > bcftools query -f'%CHROM\t%POS\t%REF,%ALT\n' reference_panel/1000GP.chr22.noNA12878.sites.vcf.gz | bgzip -c > reference_panel/1000GP.chr22.noNA12878.sites.tsv.gz
 > tabix -s1 -b2 -e2 reference_panel/1000GP.chr22.noNA12878.sites.tsv.gz

In a general pipeline, we should generate this pair of files (*.vcf.gz + *.tsv.gz) for all chromosomes separately.

3.2. Computing GLs for a single individual at specific positions

Once we have the VCF and TSV files ready, we can then start performing genotype calling for our low-coverage sample. We use the downsampled 1x BAM file generated in section 1.1 NA12878_1x_bam/NA12878.chr22.1x.bam.
BCFtools requires the reference genome in order to align the BAM file. We here use the reference genome version used by the 1000 Genomes Project (reference_genome/hs38DH.chr22.fa.gz[.fai,.gzi])
We can run the following command:

 > BAM=NA12878_1x_bam/NA12878.bam
 > VCF=reference_panel/1000GP.chr22.noNA12878.sites.vcf.gz
 > TSV=reference_panel/1000GP.chr22.noNA12878.sites.tsv.gz
 > REFGEN=reference-genome/hs38DH.chr22.fa.gz
 > OUT=NA12878_1x_vcf/NA12878.chr22.1x.vcf.gz
 > bcftools mpileup -f ${REFGEN} -I -E -a 'FORMAT/DP' -T ${VCF} -r chr22 ${BAM} -Ou | bcftools call -Aim -C alleles -T ${TSV} -Oz -o ${OUT}
 > bcftools index -f ${OUT}

Note here, that we use -T reference_panel/1000GP.chr22.noNA12878.sites.vcf.gz in the first part of the command line and
-T reference_panel/1000GP.chr22.noNA12878.sites.tsv.gz in the second part of the command line.

You may also tune the options of BCFtools to your specific needs, requirements and data (for example reducing the threshold for the minimum base/mapping quality to another positive number, in order to have less missing data).

The output of this step is a VCF file format containing genotype likelihoods at each site in the reference panel.

3.3. Merging genotype likelihoods of multiple individuals

In the case of multiple target individuals, an additional step, not required in this pipeline, is to generate a single file containing genotype likelihoods for all target individuals. This can be easily done using BCFtools, for example:

 > bcftools merge -m none -r chr22 -Oz -o merged.chr22.1x.vcf.gz -l list.txt

Where list.txt is a text file containing the full list of VCF/BCF files containing GLs of each target individual in the study, one individual file per line. The resulting file must be indexed and can be used in the subsequent steps of the pipeline.

One important step prior is to define the chunks where to run imputation and phasing. This step is not trivial because different long regions increase the running time, but small regions can drastically decrease accuracy. Also there are several variables to take into account: the amount of missingness in the defined region and the length in Mb. For these reasons we developed a tool in the GLIMPSE suite (GLIMPSE_chunk) that is able to quickly generate imputation chunks taking into account all this information.

4.1. Chunking a chromosome

We use GLIMPSE_chunk to generate imputation regions for the full chromosome 22, modifying few default parameters in this way:

 > bin/GLIMPSE_chunk --input reference_panel/1000GP.chr22.noNA12878.sites.vcf.gz --region chr22 --window-size 2000000 --buffer-size 200000 --output chunks.chr22.txt

This generates a file containing the imputation chunks and larger chunks including buffers that we will use to run GLIMPSE_phase.