<commit_before><commit_msg>add script to add HGVS info fields to a VCF file<commit_after>"""
this script add three more INFO fields to a VCF file
these three INFO fields are:
HGVS_G: HGVS genomic coordinate
HGVS_C: HGVS cDNA coordinate(s)
HGVS_P: HGVS protein coordinate(s)

Note: this script can only handle SNPs, cann't do indels
"""

import hgvs.dataproviders.uta
import hgvs.parser
import hgvs.variantmapper


def main():
	old_vcf = open("data/brca.clinvar.10line.vcf", "r")
	new_vcf = open("data/brca.clinvar.10line.addHGVS.vcf", "w")

	for line in old_vcf:
		if "#" in line:
			new_vcf.write(line)
		else:
			line = add_HGVS(line)
			new_vcf.write(line)

def add_HGVS(line):
	items = line.strip().split("\t")
	chrom = items[0]
	pos = items[1]
	ref = items[3]
	alt = items[4]

	hgvs_g = "NC_0000" + str(chrom) + ".10:g." + str(pos) + ref + ">" + alt
	print hgvs_g
	hp = hgvs.parser.Parser()
	hgvs_g = hp.parse_hgvs_variant(hgvs_g)
	hdp = hgvs.dataproviders.uta.connect()
	evm = hgvs.variantmapper.EasyVariantMapper(hdp,
        primary_assembly='GRCh37', alt_aln_method='splign')
	transcripts = evm.relevant_transcripts(hgvs_g)
	for transcript in transcripts:
		hgvs_c = evm.g_to_c(hgvs_g, transcript)
		hgvs_p = evm.c_to_p(hgvs_c)
	
	items[-1] += ";HGVS_G=" + str(hgvs_g) + ";HGVS_C=" + str(hgvs_c) + ";HGVS_P=s" + str(hgvs_p)
	new_line = "\t".join(items) + "\n"
	return new_line

if __name__ == "__main__":
	main()








{'fraction_non_alphanumeric': 0.14275717284814557, 'fraction_numerical': 0.011896431070678797, 'mean_word_length': 4.39622641509434, 'pattern_counts': {'":': 1, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2523085', 'n_tokens_mistral': 671, 'n_tokens_neox': 598, 'n_words': 133}