Answers.Body: <p>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.</p>
  <p>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, <strong>bcftools norm</strong> 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, ... .</p>
  <pre><code>bcftools norm --check-ref s --fasta-ref &lt;reference_file&gt; &lt;my_VCF_file&gt;
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-11-11T09:31:26.577Z
Answers.LastActivityDate: 2022-11-11T09:31:26.577Z
Answers.UserDisplayName: JRodrigoF
Answers.Body: <p>Thanks for @JRodrigoF for pointing towards <code>bcftools norm --check-ref s --fasta-ref</code>. This will indeed fix the initial REF prefix diff issue of my question.</p>
  <p>If anyone else stumbles over this question while facing a similar error to :<code>[E::faidx_adjust_position] The sequence &quot;1&quot; was not found faidx_fetch_seq failed at 1:69869</code></p>
  <p>Please try to re-index your FASTA-file using samtools
  <code>samtools faidx reference_file.fa</code> Most of the times this seems to solve this kind of error!</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-11-10T10:38:34.793Z
Answers.LastActivityDate: 2022-11-10T10:50:50.360Z
Answers.LastEditDate: 2022-11-10T10:50:50.360Z
Answers.UserDisplayName: mugdi
Body: <p>I lifted some unmerged trio files from hg19 to hg38 using <a href="https://github.com/sritchie73/liftOverPlink" rel="nofollow noreferrer">liftOverPlink</a>. 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.</p>
  <p>My command to merge is</p>
  <pre><code>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
  </code></pre>
  <p>First, I tried to remove those alleles in question but I realized that these are just to many to do this by hand.</p>
  <p>I used R to edit the file</p>
  <pre><code>file &lt;- read_delim(
    gzfile(
      str_glue(&quot;{some_path}file.vcf.gz&quot;)),
    delim=&quot;|&quot;,  quote='', escape_backslash = T, col_names = F) %&gt;%
  #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)) %&gt;%
    write_delim(file = str_glue(&quot;{some_path}some_file_v1.vcf&quot;),
                                   k                        col_names = F, quote = 'none', escape = 'none')
  
  #zip the file
  call = str_glue(&quot;bcftools view {some_path}some_file_v1.vcf -o {some_path}some_file_v1.vcf.gz&quot;)
  system(call)
  #index the new file
  call = str_glue(&quot;bcftools index -f {f}.gz&quot;)
  system(call)
  </code></pre>
  <p><strong>Q1: How can one fix the large amount of REF prefix differences in the chromosome files?</strong></p>
  <p><strong>Q2: If one uses the fixref plugin of bcftools which is the correct reference file to use after lifting was performed?</strong></p>
  <p><strong>EDIT:</strong> JRodrigoF recommended to use the <code>bcftools norm --check-ref s --fasta-ref GRCh38_full_analysis_set_plus_decoy_hla.fa </code> to fix the issue. It seems this produces the following error on my end:
  <code> [E::faidx_adjust_position] The sequence &quot;1&quot; was not found faidx_fetch_seq failed at 1:69869</code></p>
  <p>To make sure this is not a problem created during liftover I will post my liftover procedure here:</p>
  <pre><code>library(tidyverse)
  
  
  #set working directory to data directory
  trio_wd &lt;- str_glue(here::here(),'/trio/SomeName/')
  
  #create file list for raw data
  file_list &lt;- str_c(trio_wd,dir(trio_wd)) %&gt;% str_extract('.+txt') %&gt;%
    {.[!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,&quot;.txt&quot;)}
                        ')
        system(call)
      }
      
      #create file list for map data
      file_list &lt;- str_c(trio_wd,dir(trio_wd)) %&gt;% str_extract('.+map') %&gt;%
        {.[!is.na(.)]}
      
      for (f in file_list){
      call &lt;- str_glue('python3 reference/liftOverPlink.py --map {f} --out {str_replace(f,&quot;.map&quot;,&quot;_lifted&quot;)} --chain reference/hg19ToHg38.over.chain.gz  ')
      system(call)
      
      #create valid mapping file for hg19
      call &lt;- str_glue('python3 reference/rmBadLifts.py -m {str_replace(f,&quot;.map&quot;,&quot;_lifted&quot;)}.map --out {str_replace(f,&quot;.map&quot;,&quot;_good_lifted&quot;)}.map -l {str_replace(f,&quot;.map&quot;,&quot;_bad_lifted.dat&quot;)}')
      system(call)
      
      call &lt;- str_glue('cut -f 2 {str_replace(f,&quot;.map&quot;,&quot;_bad_lifted.dat&quot;)} &gt; {str_replace(f,&quot;.map&quot;,&quot;_exclude_lifted.dat&quot;)}')
      system(call)
      call &lt;- str_glue('cut -f 4 {str_replace(f,&quot;.map&quot;,&quot;_lifted&quot;)}.bed.unlifted | sed &quot;/^#/d&quot; &gt;&gt; {str_replace(f,&quot;.map&quot;,&quot;_exclude_lifted.dat&quot;)} ')
      system(call)
      # Note: this will clobber the lifted MAP file generated by `liftOverPlink`:
      call &lt;- str_glue('plink --file {str_remove(f,&quot;.map&quot;)} --recode --out {str_replace(f,&quot;.map&quot;,&quot;_lifted&quot;)} --exclude {str_replace(f,&quot;.map&quot;,&quot;_exclude_lifted.dat&quot;)} ')
      system(call)
      call &lt;- str_glue('plink --ped {str_replace(f,&quot;.map&quot;,&quot;_lifted.ped&quot;)} --map {str_replace(f,&quot;.map&quot;,&quot;_good_lifted&quot;)}.map --recode --out {str_replace(f,&quot;.map&quot;,&quot;_final&quot;)}')
      system(call)
      }
  </code></pre>
  <p><strong>EDIT II:</strong>
  Sample of the file including header</p>
  <pre><code>##fileformat=VCFv4.2
  ##FILTER=&lt;ID=PASS,Description=&quot;All filters passed&quot;&gt;
  ##fileDate=20221108
  ##source=PLINKv1.90
  ##contig=&lt;ID=1,length=248928329&gt;
  ##contig=&lt;ID=2,length=242157857&gt;
  ##contig=&lt;ID=3,length=198157342&gt;
  ##contig=&lt;ID=4,length=190106194&gt;
  ##contig=&lt;ID=5,length=181288140&gt;
  ##contig=&lt;ID=6,length=170595721&gt;
  ##contig=&lt;ID=7,length=159331484&gt;
  ##contig=&lt;ID=8,length=145075470&gt;
  ##contig=&lt;ID=9,length=138211490&gt;
  ##contig=&lt;ID=10,length=133659511&gt;
  ##contig=&lt;ID=11,length=135075227&gt;
  ##contig=&lt;ID=12,length=133261768&gt;
  ##contig=&lt;ID=13,length=114341522&gt;
  ##contig=&lt;ID=14,length=106873816&gt;
  ##contig=&lt;ID=15,length=101922277&gt;
  ##contig=&lt;ID=16,length=90172890&gt;
  ##contig=&lt;ID=17,length=83203771&gt;
  ##contig=&lt;ID=18,length=80252738&gt;
  ##contig=&lt;ID=19,length=58586567&gt;
  ##contig=&lt;ID=20,length=64323573&gt;
  ##contig=&lt;ID=21,length=46679699&gt;
  ##contig=&lt;ID=22,length=50776369&gt;
  ##INFO=&lt;ID=PR,Number=0,Type=Flag,Description=&quot;Provisional reference allele, may not be based on real reference genome&quot;&gt;
  ##FORMAT=&lt;ID=GT,Number=1,Type=String,Description=&quot;Genotype&quot;&gt;
  ##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
  </code></pre>
  <p><strong>EDIT III:</strong>
  Output of <code>grep '^&gt;' reference/GRCh38_full_analysis_set_plus_decoy_hla.fa</code></p>
  <pre><code>&gt;chr1  AC:CM000663.2  gi:568336023  LN:248956422  rl:Chromosome  M5:6aef897c3d6ff0c78aff06ac189178dd  AS:GRCh38
  &gt;chr2  AC:CM000664.2  gi:568336022  LN:242193529  rl:Chromosome  M5:f98db672eb0993dcfdabafe2a882905c  AS:GRCh38
  &gt;chr3  AC:CM000665.2  gi:568336021  LN:198295559  rl:Chromosome  M5:76635a41ea913a405ded820447d067b0  AS:GRCh38
  &gt;chr4  AC:CM000666.2  gi:568336020  LN:190214555  rl:Chromosome  M5:3210fecf1eb92d5489da4346b3fddc6e  AS:GRCh38
  &gt;chr5  AC:CM000667.2  gi:568336019  LN:181538259  rl:Chromosome  M5:a811b3dc9fe66af729dc0dddf7fa4f13  AS:GRCh38  hm:47309185-49591369
  </code></pre>
  <p><strong>EDIT IV:</strong> Edited the fasta-file and the corresponding index file as mentioned in the comments. Still, the error from EDIT I sustains.</p>
  <p>First five chromosomes of the modified FASTA-file:</p>
  <pre><code>grep '^&gt;' reference/GRCh38_full_analysis_set_plus_decoy_hla_modified.fa
  &gt;1  AC:CM000663.2  gi:568336023  LN:248956422  rl:Chromosome  M5:6aef897c3d6ff0c78aff06ac189178dd  AS:GRCh38
  &gt;2  AC:CM000664.2  gi:568336022  LN:242193529  rl:Chromosome  M5:f98db672eb0993dcfdabafe2a882905c  AS:GRCh38
  &gt;3  AC:CM000665.2  gi:568336021  LN:198295559  rl:Chromosome  M5:76635a41ea913a405ded820447d067b0  AS:GRCh38
  &gt;4  AC:CM000666.2  gi:568336020  LN:190214555  rl:Chromosome  M5:3210fecf1eb92d5489da4346b3fddc6e  AS:GRCh38
  &gt;5  AC:CM000667.2  gi:568336019  LN:181538259  rl:Chromosome  M5:a811b3dc9fe66af729dc0dddf7fa4f13  AS:GRCh38  hm:47309185-49591369
  </code></pre>
  <p>Addjusted the <code>bcftools</code>-call of EDIT I afterwards so the modified FASTA-file will be used</p>
  <pre><code>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 &quot;1&quot; was not found
  faidx_fetch_seq failed at 1:69869
  </code></pre>
  <p>View region in question</p>
  <pre><code>##fileformat=VCFv4.2
  ##FILTER=&lt;ID=PASS,Description=&quot;All filters passed&quot;&gt;
  ##fileDate=20221108
  ##source=PLINKv1.90
  ##contig=&lt;ID=1,length=248928329&gt;
  ##contig=&lt;ID=2,length=242157857&gt;
  ##contig=&lt;ID=3,length=198157342&gt;
  ##contig=&lt;ID=4,length=190106194&gt;
  ##contig=&lt;ID=5,length=181288140&gt;
  ##contig=&lt;ID=6,length=170595721&gt;
  ##contig=&lt;ID=7,length=159331484&gt;
  ##contig=&lt;ID=8,length=145075470&gt;
  ##contig=&lt;ID=9,length=138211490&gt;
  ##contig=&lt;ID=10,length=133659511&gt;
  ##contig=&lt;ID=11,length=135075227&gt;
  ##contig=&lt;ID=12,length=133261768&gt;
  ##contig=&lt;ID=13,length=114341522&gt;
  ##contig=&lt;ID=14,length=106873816&gt;
  ##contig=&lt;ID=15,length=101922277&gt;
  ##contig=&lt;ID=16,length=90172890&gt;
  ##contig=&lt;ID=17,length=83203771&gt;
  ##contig=&lt;ID=18,length=80252738&gt;
  ##contig=&lt;ID=19,length=58586567&gt;
  ##contig=&lt;ID=20,length=64323573&gt;
  ##contig=&lt;ID=21,length=46679699&gt;
  ##contig=&lt;ID=22,length=50776369&gt;
  ##INFO=&lt;ID=PR,Number=0,Type=Flag,Description=&quot;Provisional reference allele, may not be based on real reference genome&quot;&gt;
  ##FORMAT=&lt;ID=GT,Number=1,Type=String,Description=&quot;Genotype&quot;&gt;
  ##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
  </code></pre>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-11-09T09:44:57.110Z
Comments.Text: 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
Comments.UserDisplayName: JRodrigoF
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-11-09T10:12:53.007Z
Comments.Text: 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.
Comments.UserDisplayName: mugdi
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-11-09T15:43:14.470Z
Comments.Text: 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,
Comments.UserDisplayName: JRodrigoF
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-11-09T15:58:52.367Z
Comments.Text: Added a sample of the file including the header. Does not seem that this is the case for me.
Comments.UserDisplayName: mugdi
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-11-09T18:17:38.460Z
Comments.Text: That error just means you have `1` in your VCF and probably `chr1` in your fasta sequence.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-11-09T19:50:27.660Z
Comments.Text: How can I solve this error?
Comments.UserDisplayName: mugdi
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-11-09T21:17:32.283Z
Comments.Text: can you write here or include in the post the names of the chromosomes in your fasta file (reference). You can fetch the first 5 with grep '^>' <my_reference_file>
Comments.UserDisplayName: JRodrigoF
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-11-09T21:48:27.527Z
Comments.Text: Done. So I guess I have to edit the fasta file and remove the "chr" string infront of each position?
Comments.UserDisplayName: mugdi
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-11-09T22:02:59.760Z
Comments.Text: 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,
Comments.UserDisplayName: JRodrigoF
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-11-10T10:00:04.030Z
Comments.Text: Done. Also adjusted the index file of the fast file in the same way. Still, I receive the same error as before. I added the the output of the region in question to the section EDIT IV of the question.
Comments.UserDisplayName: mugdi
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-11-10T20:36:33.440Z
Comments.Text: Good. I'll wrap it up as an answer too and maybe you could mark it as such so that I earn a few points,
Comments.UserDisplayName: JRodrigoF
Tags: <bcftools><liftover><merge>
Title: Large amount of REF prefix differences after lifting
ThreadText: # Large amount of REF prefix differences after lifting
  
  (Asked by: username_0 on 2022-11-08)
  
  I lifted some unmerged trio files from hg19 to hg38 using [liftOverPlink](https://github.com/sritchie73/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:** username_1 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
  
  ```
  
  **Tags:** bcftools, liftover, merge
  
  > **Comments:**
  > **username_1**: 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  
  > **username_0**: 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.  
  > **username_1**: 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,  
  > **username_0**: Added a sample of the file including the header. Does not seem that this is the case for me.  
  > **username_2**: That error just means you have `1` in your VCF and probably `chr1` in your fasta sequence.  
  > **username_0**: How can I solve this error?  
  > **username_1**: can you write here or include in the post the names of the chromosomes in your fasta file (reference). You can fetch the first 5 with grep '^>' <my_reference_file>  
  > **username_0**: Done. So I guess I have to edit the fasta file and remove the "chr" string infront of each position?  
  > **username_1**: 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,  
  > **username_0**: Done. Also adjusted the index file of the fast file in the same way. Still, I receive the same error as before. I added the the output of the region in question to the section EDIT IV of the question.  
  > **username_1**: Good. I'll wrap it up as an answer too and maybe you could mark it as such so that I earn a few points,  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_1 on 2022-11-11. Score: 1 | Accepted answer
  
  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>
  
  ```
  
  
  ---
  
  ### Answer by username_0 on 2022-11-10. Score: 0
  
  Thanks for @username_1 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!
  
  
  ---
  
  