Large amount of REF prefix differences after lifting

Question

I lifted some unmerged trio files from hg19 to hg38 using liftOverPlink. It looks like that several houndred allels on each chromosome of one file shifted/switched to non-standard alleles. I recorded some of the differences. And will add them to the post. I will also add a sample of the file to this post. My goal is to merge the lifted trio files to one singe file.

My command to merge is

bcftools  merge -l merge.txt --force-samples -Ob -o file.bcf


The REF prefixes differ: C vs N (1,1)
Failed to merge alleles at 1:3940550 in file_1.vcf.gz

First, I tried to remove those alleles in question but I realized that these are just to many to do this by hand.

I used R to edit the file

file <- read_delim(
  gzfile(
    str_glue("{some_path}file.vcf.gz")),
  delim="|",  quote='', escape_backslash = T, col_names = F) %>%
#A list of all the REF prefix differences: number equals position in chromsome 1
  #630128 N vs A, 791532, N vs A, 841166 G vs A, 937572 C vs T, 1110226 N vs A, 1486611 N vs G, 1554290 C vs G,
  #1654663 N vs T, 1671724 C vs T, 1755523 C vs T, 1774697 C vs T, 1916092 C vs T, 1949922 C vs T, 1959595 C vs n
  #1959717 C vs T, 1983048 A vs G, 1991690 A vs G, 1999315 C vs T 2096621 C vs T 2107385 A vs G, 2124676 A vs G
  #2274142 A vs C, 2284952 A vs G, 2293812 C vs T, 2303949 N vs C, 2321888 T vs N, 2351827 C vs T,2364905 N vs G
  #2364905 T vs C, 2420707 T vs C, 2431388 T vs C, 2441501 T vs C, 2461209 C vs G, 2614790 A vs G, 2620461 G vs A
  #2622185 C vs T, 2792599 A vs C, 2861494 C vs T 2864702 G vs N, 2884008 T vs C , 2894837 A vs C, 2925918 T vs C,
  #2944366 T vs C, 3097126 G vs A,, 3107701 C vs A, 3109422  C vs A, 3112423 G vs A, , 3119756 C vs T, 3146311 G vs T
  #3282557 G vs A,, 3294386 A vs G, 3294894 A vs C, 3303337 A vs G, 3320610 T vs C, 3335540 T vs C, 3344952 T vs C
  #3345118 T vs C 3355977 G vs A, 3401871 G vs A, 3443779 G vs T,, 3534441 T vs N, 3591297 A vs G,
  #3595370  A vs G, 3601923 T vs C, 3624792 G vs N, 3689917 T vs C, 3689917 G vs A, 3695328 G vs A, 3696503 G vs A
  #3762340 T vs C, 3775163 A vs C, 3870048 T vs Cl 3924314 G vs A3935847 T vs C, 3940550 G vs A
   filter(!grepl('.+630128.+|.+791532.+|.+841166.+|.+937572.+|.+1110226.+|.+1486611.+|.+1554290.+|.+1654663.+|.+1671724.+|.+1755523.+|.+1774697.+|.+1916092.+|.+1949922.+|.+1959595.+|.+1959717.+|.+1983048.+|.+1991690.+|.+1999315.+|.+2096621.+|.+2107385.+|.+2124676.+|.+2274142.+|.+2284952.+|.+2293812.+|.+2303949.+|.+2321888.+|.+2324778.+|.+2351678.+|.+2351827.+|.+2364905.+|.+2364905.+|.+2420707.+|.+2431388.+|.+2441501.+|.+2461209.+|.+2614790.+|.+2620461.+|.+2622185.+|.+2792599.+|.+2861494.+|.+2864702.+|.+2884008.+|.+2894837.+|.+2925918.+|.+2944366.+|.+3097126.+|.+3107701.+|.+3109422.+|.+3112423.+|.+3119756.+|.+3146311.+|.+3282557.+|.+3294386.+|.+3294894.+|.+3303337.+|.+3320610.+|.+3335540.+|.+3344952.+|.+3345118.+|.+3355977.+|.+3401871.+|.+3443779.+|.+3534441.+|.+3591297.+|.+3595370.+|.+3601923.+|.+3624792.+|.+3689917.+|.+3689917.+|.+3695328.+|.+3696503.+|.+3762340.+|.+3775163.+|.+3870048.+|.+3924314.+|.+3935847.+|.+3940550.+',X1)) %>%
  write_delim(file = str_glue("{some_path}some_file_v1.vcf"),
                                 k                        col_names = F, quote = 'none', escape = 'none')

#zip the file
call = str_glue("bcftools view {some_path}some_file_v1.vcf -o {some_path}some_file_v1.vcf.gz")
system(call)
#index the new file
call = str_glue("bcftools index -f {f}.gz")
system(call)

Q1: How can one fix the large amount of REF prefix differences in the chromosome files?

Q2: If one uses the fixref plugin of bcftools which is the correct reference file to use after lifting was performed?

EDIT: JRodrigoF recommended to use the bcftools norm --check-ref s --fasta-ref GRCh38_full_analysis_set_plus_decoy_hla.fa to fix the issue. It seems this produces the following error on my end: [E::faidx_adjust_position] The sequence "1" was not found faidx_fetch_seq failed at 1:69869

To make sure this is not a problem created during liftover I will post my liftover procedure here:

library(tidyverse)


#set working directory to data directory
trio_wd <- str_glue(here::here(),'/trio/SomeName/')

#create file list for raw data
file_list <- str_c(trio_wd,dir(trio_wd)) %>% str_extract('.+txt') %>%
  {.[!is.na(.)]}


    #generate map files from 23andme text file
    for(f in file_list){
      call = str_glue('plink --23file {f} --recode --out {str_remove(f,".txt")}
                      ')
      system(call)
    }
    
    #create file list for map data
    file_list <- str_c(trio_wd,dir(trio_wd)) %>% str_extract('.+map') %>%
      {.[!is.na(.)]}
    
    for (f in file_list){
    call <- str_glue('python3 reference/liftOverPlink.py --map {f} --out {str_replace(f,".map","_lifted")} --chain reference/hg19ToHg38.over.chain.gz  ')
    system(call)
    
    #create valid mapping file for hg19
    call <- str_glue('python3 reference/rmBadLifts.py -m {str_replace(f,".map","_lifted")}.map --out {str_replace(f,".map","_good_lifted")}.map -l {str_replace(f,".map","_bad_lifted.dat")}')
    system(call)
    
    call <- str_glue('cut -f 2 {str_replace(f,".map","_bad_lifted.dat")} > {str_replace(f,".map","_exclude_lifted.dat")}')
    system(call)
    call <- str_glue('cut -f 4 {str_replace(f,".map","_lifted")}.bed.unlifted | sed "/^#/d" >> {str_replace(f,".map","_exclude_lifted.dat")} ')
    system(call)
    # Note: this will clobber the lifted MAP file generated by `liftOverPlink`:
    call <- str_glue('plink --file {str_remove(f,".map")} --recode --out {str_replace(f,".map","_lifted")} --exclude {str_replace(f,".map","_exclude_lifted.dat")} ')
    system(call)
    call <- str_glue('plink --ped {str_replace(f,".map","_lifted.ped")} --map {str_replace(f,".map","_good_lifted")}.map --recode --out {str_replace(f,".map","_final")}')
    system(call)
    }

EDIT II: Sample of the file including header

##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##fileDate=20221108
##source=PLINKv1.90
##contig=<ID=1,length=248928329>
##contig=<ID=2,length=242157857>
##contig=<ID=3,length=198157342>
##contig=<ID=4,length=190106194>
##contig=<ID=5,length=181288140>
##contig=<ID=6,length=170595721>
##contig=<ID=7,length=159331484>
##contig=<ID=8,length=145075470>
##contig=<ID=9,length=138211490>
##contig=<ID=10,length=133659511>
##contig=<ID=11,length=135075227>
##contig=<ID=12,length=133261768>
##contig=<ID=13,length=114341522>
##contig=<ID=14,length=106873816>
##contig=<ID=15,length=101922277>
##contig=<ID=16,length=90172890>
##contig=<ID=17,length=83203771>
##contig=<ID=18,length=80252738>
##contig=<ID=19,length=58586567>
##contig=<ID=20,length=64323573>
##contig=<ID=21,length=46679699>
##contig=<ID=22,length=50776369>
##INFO=<ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on real reference genome">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##bcftools_viewVersion=1.13+htslib-1.13+ds
##bcftools_viewCommand=view -r 1 -o /home/oem/Documents/GitHub/Genomics_prac_guide/trio/some_name/some_name//genome_some_name_v5_Full_20191220113716_final_1.vcf.gz /home/oem/Documents/GitHub/Genomics_prac_guide/trio/some_name/some_name_v5_Full_20191220113716_final.vcf.gz; Date=Tue Nov  8 13:02:01 2022
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  FAM001_ID001
1   69869   rs548049170 T   .   .   .   PR  GT  0/0
1   74792   rs13328684  N   .   .   .   PR  GT  ./.
1   792461  rs116587930 G   .   .   .   PR  GT  0/0
1   817341  rs3131972   G   .   .   .   PR  GT  0/0
1   818725  rs12184325  C   .   .   .   PR  GT  0/0
1   820888  rs12567639  A   .   .   .   PR  GT  0/0
1   823656  rs114525117 G   .   .   .   PR  GT  0/0
1   858952  rs12127425  G   .   .   .   PR  GT  0/0
1   866156  rs79373928  T   .   .   .   PR  GT  0/0
1   880041  rs72888853  N   .   .   .   PR  GT  ./.
1   889018  rs7538305   A   .   .   .   PR  GT  0/0
1   894801  rs28444699  A   .   .   .   PR  GT  0/0
1   895351  i713449 N   .   .   .   PR  GT  ./.
1   899450  rs116452738 G   .   .   .   PR  GT  0/0
1   899712  rs72631887  T   .   .   .   PR  GT  0/0
1   903285  rs28678693  T   C   .   .   PR  GT  0/1
1   905373  rs4970382   T   C   .   .   PR  GT  0/1
1   911428  rs4475691   C   .   .   .   PR  GT  0/0
1   916010  rs72631889  G   .   .   .   PR  GT  0/0

EDIT III: Output of grep '^>' reference/GRCh38_full_analysis_set_plus_decoy_hla.fa

>chr1  AC:CM000663.2  gi:568336023  LN:248956422  rl:Chromosome  M5:6aef897c3d6ff0c78aff06ac189178dd  AS:GRCh38
>chr2  AC:CM000664.2  gi:568336022  LN:242193529  rl:Chromosome  M5:f98db672eb0993dcfdabafe2a882905c  AS:GRCh38
>chr3  AC:CM000665.2  gi:568336021  LN:198295559  rl:Chromosome  M5:76635a41ea913a405ded820447d067b0  AS:GRCh38
>chr4  AC:CM000666.2  gi:568336020  LN:190214555  rl:Chromosome  M5:3210fecf1eb92d5489da4346b3fddc6e  AS:GRCh38
>chr5  AC:CM000667.2  gi:568336019  LN:181538259  rl:Chromosome  M5:a811b3dc9fe66af729dc0dddf7fa4f13  AS:GRCh38  hm:47309185-49591369

EDIT IV: Edited the fasta-file and the corresponding index file as mentioned in the comments. Still, the error from EDIT I sustains.

First five chromosomes of the modified FASTA-file:

grep '^>' reference/GRCh38_full_analysis_set_plus_decoy_hla_modified.fa
>1  AC:CM000663.2  gi:568336023  LN:248956422  rl:Chromosome  M5:6aef897c3d6ff0c78aff06ac189178dd  AS:GRCh38
>2  AC:CM000664.2  gi:568336022  LN:242193529  rl:Chromosome  M5:f98db672eb0993dcfdabafe2a882905c  AS:GRCh38
>3  AC:CM000665.2  gi:568336021  LN:198295559  rl:Chromosome  M5:76635a41ea913a405ded820447d067b0  AS:GRCh38
>4  AC:CM000666.2  gi:568336020  LN:190214555  rl:Chromosome  M5:3210fecf1eb92d5489da4346b3fddc6e  AS:GRCh38
>5  AC:CM000667.2  gi:568336019  LN:181538259  rl:Chromosome  M5:a811b3dc9fe66af729dc0dddf7fa4f13  AS:GRCh38  hm:47309185-49591369

Addjusted the bcftools-call of EDIT I afterwards so the modified FASTA-file will be used

call = str_glue('bcftools  norm --check-ref s {f} --fasta-ref reference/GRCh38_full_analysis_set_plus_decoy_hla_modified.fa -o {out} ')
[E::faidx_adjust_position] The sequence "1" was not found
faidx_fetch_seq failed at 1:69869

View region in question

##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##fileDate=20221108
##source=PLINKv1.90
##contig=<ID=1,length=248928329>
##contig=<ID=2,length=242157857>
##contig=<ID=3,length=198157342>
##contig=<ID=4,length=190106194>
##contig=<ID=5,length=181288140>
##contig=<ID=6,length=170595721>
##contig=<ID=7,length=159331484>
##contig=<ID=8,length=145075470>
##contig=<ID=9,length=138211490>
##contig=<ID=10,length=133659511>
##contig=<ID=11,length=135075227>
##contig=<ID=12,length=133261768>
##contig=<ID=13,length=114341522>
##contig=<ID=14,length=106873816>
##contig=<ID=15,length=101922277>
##contig=<ID=16,length=90172890>
##contig=<ID=17,length=83203771>
##contig=<ID=18,length=80252738>
##contig=<ID=19,length=58586567>
##contig=<ID=20,length=64323573>
##contig=<ID=21,length=46679699>
##contig=<ID=22,length=50776369>
##INFO=<ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on real reference genome">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##bcftools_viewVersion=1.13+htslib-1.13+ds
##bcftools_viewCommand=view -r 1 -o /home/oem/Documents/GitHub/Genomics_prac_guide/trio/Some_name/Some_name//genome_Some_name_v5_Full_20191220113716_final_1.vcf.gz /home/oem/Documents/GitHub/Genomics_prac_guide/trio/Some_name/genome_Some_name_v5_Full_20191220113716_final.vcf.gz; Date=Tue Nov  8 13:02:01 2022
##bcftools_viewCommand=view -o /home/oem/Documents/GitHub/Genomics_prac_guide/trio/Some_name/Some_name/genome_Some_name_v5_Full_20191220113716_final_1_1.vcf.gz /home/oem/Documents/GitHub/Genomics_prac_guide/trio/Some_name/Some_name/genome_Some_name_v5_Full_20191220113716_final_1_1.vcf; Date=Tue Nov  8 18:16:55 2022
##bcftools_viewCommand=view -r 1:69869 /home/oem/Documents/GitHub/Genomics_prac_guide/trio/Some_name/Some_name/genome_Some_name_v5_Full_20191220113716_final_1_1.vcf.gz; Date=Thu Nov 10 10:56:54 2022
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  FAM001_ID001
1   69869   rs548049170 T   .   .   .   PR  GT  0/0

What about using bctools norm instead to fix the REF allele whenever needed? bcftools norm --check-ref s --fasta-ref GRCh38_full_analysis_set_plus_decoy_hla.fa Or any other flavour of GRCh38 you are using, since the issue to solve is after the lifting over, so v38 has to be used in this case — JRodrigoF, Commented Nov 9, 2022 at 9:44
If I do this I will recieve the following error msg: [E::faidx_adjust_position] The sequence "1" was not found faidx_fetch_seq failed at 1:69869 Maybe something went wrong during the lifting? I will add my lifting procedure to the OP. — mugdi, Commented Nov 9, 2022 at 10:12
looks to me that the header of your lifted VCF now in hg38 is not compatible with the sequence names contained in your reference fasta file; for example if in the hg38 version you have names like "chr1", ... and in the new VCF you have "1", an error like this would come up, — JRodrigoF, Commented Nov 9, 2022 at 15:43
That error just means you have 1 in your VCF and probably chr1 in your fasta sequence. — terdon, Commented Nov 9, 2022 at 18:17
Yes. You can do things like this with sed command. sed 's/>chr/>/' reference/GRCh38_full_analysis_set_plus_decoy_hla.fa > reference/GRCh38_full_analysis_set_plus_decoy_hla_modified.fa Also I would check if the plink lifter does not have any relevant option regarding this issue. My expectation would be that the lifter produces a new VCF that contains the correct contig names as taken from the new reference, ... but here is not being the case. I would usually go for the lifter tool from GATK for example, — JRodrigoF, Commented Nov 9, 2022 at 22:02

JRodrigoF · Accepted Answer · 2022-11-11 09:31:26Z

It is important that after lifting over a VCF/BCF file from one coordinate system to another, the header of the resulting file contains the chromosome/contig names as present in the coordinate system of the new reference genome. it works as a double check that the assumptions of the coordinate system (particular reference) are consistent with the actual data (VCF lines). PLINK files for example do not have a header section and this doesn't not apply.

Looks like BCFTools (merge) complains whenever if finds positions that have one base ('C') as reference in one of the input files, and another base ('N') in another file. When working with VCF/BCF files, bcftools norm can be used to fix 'N' reference position into their actual base (if that is the case) and/or fix bases listed as reference (e.g. coming from a PLINK major-allele based BED file), but which actually correspond to an alternate allele based on either GRCh37/GRCh38, ... .

bcftools norm --check-ref s --fasta-ref <reference_file> <my_VCF_file>

mugdi · Accepted Answer · 2022-11-10 10:50:50Z

0

Thanks for @JRodrigoF for pointing towards bcftools norm --check-ref s --fasta-ref. This will indeed fix the initial REF prefix diff issue of my question.

If anyone else stumbles over this question while facing a similar error to :[E::faidx_adjust_position] The sequence "1" was not found faidx_fetch_seq failed at 1:69869

Please try to re-index your FASTA-file using samtools samtools faidx reference_file.fa Most of the times this seems to solve this kind of error!

edited Nov 10, 2022 at 10:50

answered Nov 10, 2022 at 10:38

mugdi

1455 bronze badges

Add a comment |

Stack Exchange Network

Large amount of REF prefix differences after lifting

2 Answers 2

Your Answer

Not the answer you're looking for? Browse other questions tagged
bcftools
liftover
merge
or ask your own question.

Linked

Hot Network Questions

Large amount of REF prefix differences after lifting

2 Answers 2

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged bcftoolsliftovermerge or ask your own question.

Linked

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
bcftools
liftover
merge
or ask your own question.