1
$\begingroup$

Here are the first few columns of my VCF v4.2 file. As you can see, the alt allele is missing. The VCF's refs and alts don't correspond to the reference genome.

#CHROM POS ID REF ALT
9 77051 rs9408670 T C
9 114350 rs12685480 A .
9 116581 rs558835286 G .
9 117550 rs267602077 A .
9 117720 rs10959293 T .
9 126903 rs62533602 G A
9 133828 rs78804900 C T
9 154736 rs200683662 T .
9 154933 rs17795008 T G
9 174425 rs2998944 T C
9 174500 rs142581442 G .
9 175050 rs62532020 G A
9 185632 rs2992854 C A
9 188894 rs3951559 G A
9 191843 rs3008128 G A
9 202061 rs117675283 T .
9 204201 rs10964134 G A

The missing alt allele happens when the genotype is 0/0. The VCF file isn't saying anything incorrect. However, there is a tool I'm trying to use that filters out "0/."genotypes.

A better way of formatting them (which I'm trying to achieve) would have the ref/alt alleles correspond to the reference genome, and just have a redundant letter in the case of a homozygous genotype.

I'd like my VCF to display the correct alt/ref for that position even if there are no samples with a genotype corresponding to the alt allele. How do I do this? I'm not sure what this process is called.

Ideal output:

#CHROM POS ID REF ALT
9 77051 rs9408670 T C
9 114350 rs12685480 A G
9 116581 rs558835286 G A
$\endgroup$
4
  • $\begingroup$ Are you sure you need to update the alt? The problem with 0/. genotypes is that period usually means “missing” so that genotype is kinda bizarre. There’s one allele we know, the reference, but also one we couldn’t figure out? If you think these should actually be homozygous reference genotypes, convert them to “0/0”. $\endgroup$ Jul 1, 2023 at 22:07
  • $\begingroup$ If you mean the tool can’t handle “.” in the ALT column, maybe just use awk to update all the dots to “N”? That’s a valid VCF allele but won’t be a ref. $\endgroup$ Jul 1, 2023 at 22:13
  • $\begingroup$ It would still filter that out, unfortunately. My solution for now is to create valid letters (ATGC) for the ref, replacing dots. $\endgroup$
    – BigMistake
    Jul 2, 2023 at 4:55
  • $\begingroup$ You should include the genotype column of the VCF, so that we have a clearer picture. As has been mentioned, it is not the same a missing (./.) genotype than a homozigous reference. The former representation usually means the actual genotype could not be establish by the variant caller and assuming hom ref might cause mistakes in dowstream analyses a $\endgroup$
    – JRodrigoF
    Jul 5, 2023 at 13:30

1 Answer 1

1
$\begingroup$

Let me know if you see any potential issues with this solution. It just replaces the dots with random letters (easier/faster than aligning to reference genome).

def add_ref():
    # Directory path
    dir_path = '/path/'

    # Function to process genotypes in each row
    def process_genotypes(row):
        genotypes = np.array([x.split(":")[0] for x in row[9:]])
        if row[4] == "." and not any(val in ["1/1", "0/1", "1/0"] for val in genotypes): # to make sure I'm not introducing an allele that isn't present
            # also can do if row[3] in ["A", "T"]: then random.choice(["C", "G"])
            if row[3] == "A":
                row[4] = "C"
            elif row[3] == "T":
                row[4] = "G"
            elif row[3] == "C":
                row[4] = random.choice(["T"]) # I know
            elif row[3] == "G":
                row[4] = random.choice(["A"])
        return row

    # Get all .vcf.gz files in the directory
    files = glob.glob(dir_path + '*.vcf.gz')

    for file in files:
        print(f"Processing {file}...")
        # Name of the new file
        new_file = os.path.join(dir_path, 'new_' + os.path.basename(file))
        with gzip.open(file, 'rt') as in_f, gzip.open(new_file, 'wt') as out_f:
            for line in in_f:
                if line.startswith('#'):
                    out_f.write(line)
                else:
                    row = line.strip().split('\t')
                    new_row = process_genotypes(row)
                    out_f.write('\t'.join(new_row) + '\n')
$\endgroup$
3
  • 1
    $\begingroup$ Okay I see. Python pandas is an alternative. If this approach runs out of memory its worth thinking about. $\endgroup$
    – M__
    Jul 2, 2023 at 10:02
  • $\begingroup$ You best make sure to that the downstream program doesn’t use the alt allele field for something, or else it could cause problems $\endgroup$
    – user438383
    Jul 2, 2023 at 11:25
  • $\begingroup$ One of the downstream programs is imputation. Not sure if it uses the alt, but I don't see why it would. $\endgroup$
    – BigMistake
    Jul 3, 2023 at 19:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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