Strand Alignment in ShapeIt -- "Reference and Main panels are not well aligned"

Question

[whit18]

I followed the code on ShapeIt's website (http://www.shapeit.fr/pages/m03_phasing/imputation.html) to check strand alignment, and the results don't make sense, with the vast majority of the SNPs being listed as missing in the reference panel. The output is below in case it's helpful. I confirmed they're both hg37 build, and I'm at a loss as to what else to try. Any suggestions would be greatly appreciated!!

---

Segmented HAPlotype Estimation & Imputation Tool * Authors : Olivier Delaneau, Jared O'Connell, Jean-François Zagury, Jonathan Marchini * Contact : send an email to the OXSTATGEN mail list https://www.jiscmail.ac.uk/cgi-bin/webadmin?A0=OXSTATGEN * Webpage : https://mathgen.stats.ox.ac.uk/shapeit * Version : v2.r837 * Date : 08/01/2018 22:13:07 * LOGfile : [180108_chr10.alignments.log]

MODE -summarise : GENERATING SUMMARY STATISTICS OF THE INPUT DATA * Autosome (chr1 ... chr22) * Reference panel of haplotypes used

Parameters : * Seed : 1515467587 * Parallelisation: 1 threads * Ref allele is NOT aligned on the reference genome

Reading site list in [lifted010818.10.map] * 25649 sites included

Reading sample list and genotypes in [lifted010818.10.ped] where missing-code = [0] * 4028 samples included * 4028 unrelateds / 0 duos / 0 trios in 4028 different families

Reading sample list [/ysm-gpfs/datasets/genomes/1000Genomes/1000GP_Phase3/1000GP_Phase3/1000GP_Phase3.sample] * 5008 reference haplotypes included

Reading SNPs in [/ysm-gpfs/datasets/genomes/1000Genomes/1000GP_Phase3/1000GP_Phase3/1000GP_Phase3_chr10.legend.gz] * 224 reference panel sites included * 4013234 reference panel sites excluded

ERROR: Reference and Main panels are not well aligned: * #Missing sites in reference panel = 24873 * #Misaligned sites between panels = 552 * #Multiple alignments between panels = 0

[Runen]

Have you solved the problem? I also encountered the same problem. look forward to your reply

[Dustin]

I have been trying to solve this same issue for about 2 weeks now. I have done some troubleshooting and the positions in my sample file are in the haplotype reference panel. I can’t figure how my reference and main panels can be misaligned when they have 80% of the same sites between the two of them (according to my tests with my own scripts).

[Gennady Khvorykh]

I had a data set with unknown human genome assembly. I figure out it by intersecting the positions taken from dbSNP130 (hg18) and dbSNP150 (hg19). With incorrect assembly only 0.9% of positions were in common. According to SHAPEIT -check, you have 0.9% of SNPs in reference. Why not to double check the assemblies of data and reference panel?

Answer 1

[Kevin Blighe]

Edit: November 24, 2019

Zoom down to code for SHAPEIT2, here:

• C: Strand Alignment in ShapeIt -- "Reference and Main panels are not well aligned"

Another answer posted here:

• C: Phasing with SHAPEIT

The genetic maps can be download via the links, HERE.

-----------------------------------------------

Original answer:

You should first check alignment of the datasets:

Strand alignment

This is a crucial step of prephasing/imputation to make sure that the GWAS dataset is well aligned with the reference panel of haplotypes. You can check SNP alignment in two steps with Plink (step1 and step2) or with GTOOL.

[source: https://mathgen.stats.ox.ac.uk/genetics_software/shapeit/shapeit.html]

After you do that, you may have to remove variants from both datasets.

Then, follow the steps here: Phasing with a reference panel.

Basically, you have to run the program once with the -check flag so that it can check how aligned are the datasets. This initial check will produce a file with an 'exclude' extension, which comprises problematic variants. This file should then be passed to the --exclude-snp parameter when you run the actual pre-phasing step.

Kevin

[Dustin]

Hey Kevin,

Thank you for your response. I've tried it with my data prior to alignment, and after alignment. ShapeIt2 says there are no aligned sites between my data and the reference whether I do strand alignment or not. For strand alignment I used Genotype Harmonizer which gives the following output for chr22 of my data:

"Input data loaded

Reference data loaded

Beginning alignment

Iteration 1 - 10,000 variants processed
Iteration 1 - 20,000 variants processed
Iteration 1 - 30,000 variants processed
Iteration 1 - 40,000 variants processed
Iteration 1 - 50,000 variants processed
Iteration 1 - 60,000 variants processed
Iteration 1 - Completed, non A/T and non G/C SNPs are aligned 36,749 found and 0 swapped
Iteration 2 - Skipped, non A/T and non G/C SNPs are not LD checked 
Iteration 3 - 10,000 variants processed (1,468 G/C or A/T SNPs checked)
Iteration 3 - 20,000 variants processed (2,953 G/C or A/T SNPs checked)
Iteration 3 - 30,000 variants processed (4,407 G/C or A/T SNPs checked)
Iteration 3 - 40,000 variants processed (5,897 G/C or A/T SNPs checked)
Iteration 3 - Completed, non A/T and non G/C SNPs are aligned. Extra LD check skipped

Swapped 32 A/T or G/C variants based on LD patterns

Alignment complete

Excluded in total 22,173 variants during alignment phase

Writing results

Output data written

Program complete"

So during alignment, 22,172 of my 60,000+ variants in chr22 were not found in the reference. Then when I run ShapeIt2 I get:

"Segmented HAPlotype Estimation & Imputation Tool
  * Authors : Olivier Delaneau, Jared O'Connell, Jean-François Zagury, Jonathan Marchini
  * Contact : send an email to the OXSTATGEN mail list https://www.jiscmail.ac.uk/cgi-bin/webadmin?A0=OXSTATGEN
  * Webpage : https://mathgen.stats.ox.ac.uk/shapeit
  * Version : v2.r904
  * Date    : 27/10/2019 21:41:47
  * LOGfile : [trio_chr22_harmonized.checks.log]

MODE -summarise : GENERATING SUMMARY STATISTICS OF THE INPUT DATA
  * Autosome (chr1 ... chr22)
  * Reference panel of haplotypes used


Parameters :
  * Seed : 1572212507
  * Parallelisation: 1 threads
  * Ref allele is NOT aligned on the reference genome

Reading site list in [trio_chr22_harmonized.bim]
  * 40639 sites included

Reading sample list in [trio_chr22_harmonized.fam]
  * 3 samples included
  * 0 unrelateds / 0 duos / 1 trios in 1 different families

Reading genotypes in [trio_chr22_harmonized.bed]
  * Plink binary file SNP-major mode

Reading sample list [/references/1000GP_Phase3/1000GP_Phase3_chr22.legend]
  * 2220434 reference haplotypes included

Reading SNPs in [/references/1000GP_Phase3/1000GP_Phase3_chr22.haps]
  * 0 reference panel sites included
  * 1110216 reference panel sites excluded

ERROR: Reference and Main panels are not well aligned:
  * #Missing sites in reference panel = 40639
  * #Misaligned sites between panels = 0
  * #Multiple alignments between panels = 0"

So this tells me that all of my 40639 sites in my trio are not found in the reference panel. However, according to Genotype Harmonizer all of my sites (after it removed 22,173) should be in the reference panel and according to my own scripts my sites are in the reference panel. Also worth mentioning, my original VCF files were aligned to GRCh38 but I used liftover to convert the sites to GRCh37. So I know that my positions are in GRCh37 coordinates.

[Kevin Blighe]

I see... can you show the code that you are using? I mean, I literally run SHAPEIT at least 3 times.

With this first code chunk, an error is thrown but it generates the 'exclude' file:

shapeit \
      -check \
      -B temp \
      -M library/1000GP_Phase3/genetic_map_chr"${chr}"_combined_b37.txt \
      --input-ref library/1000GP_Phase3/1000GP_Phase3_chr"${chr}".hap.gz library/1000GP_Phase3/1000GP_Phase3_chr"${chr}".legend.gz library/1000GP_Phase3/1000GP_Phase3.sample \
      --output-log Prephased/GSA_QCd_chr"${chr}"_alignments \
      -T 8 ;

With this second code chunk, the check is re-performed, but excluding the problematic variants (this runs to completion, without errors):

shapeit \
    -check \
    -B temp \
    -M library/1000GP_Phase3/genetic_map_chr"${chr}"_combined_b37.txt \
    --input-ref library/1000GP_Phase3/1000GP_Phase3_chr"${chr}".hap.gz library/1000GP_Phase3/1000GP_Phase3_chr"${chr}".legend.gz library/1000GP_Phase3/1000GP_Phase3.sample \
    --exclude-snp Prephased/GSA_QCd_chr"${chr}"_alignments.snp.strand.exclude \
    -T 8 ;

Now I pre-phase:

shapeit \
    -B temp \
    -M library/1000GP_Phase3/genetic_map_chr"${chr}"_combined_b37.txt \
    --input-ref library/1000GP_Phase3/1000GP_Phase3_chr"${chr}".hap.gz library/1000GP_Phase3/1000GP_Phase3_chr"${chr}".legend.gz library/1000GP_Phase3/1000GP_Phase3.sample \
    --exclude-snp Prephased/GSA_QCd_chr"${chr}"_alignments.snp.strand.exclude \
    -O Prephased/GSA_QCd_chr"${chr}"_1KGphased \
    -T 8 ;

[Dustin]

I use the same ShapeIt code as you, with one difference. I was putting "--input-ref file.sample file.hap.gz file.legend.gz" instead of "--input-ref file.hap.gz file.legend.gz file.sample". The order made the difference and now ShapeIt is recognizing that my sample sites are in the reference. Seems crazy that the order of input matters that much! Thanks for your help Kevin.