renluyao
/
high_confidence_calls_manuscript
Watch 1
Star 0
Fork 0
Code Issues 0 Pull Requests 0 Releases 0 Wiki Activity
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
10 Commits
1 Branch
Tree: 3ccd66f2ea
Branches Tags
${ item.name }
Create branch ${ searchTerm }
from '3ccd66f2ea'
${ noResults }
high_confidence_calls_manus.../codescripts/merge_mendelian_vcfinfo.py

128 lines
5.0KB
Raw Blame History 

  1. from __future__ import division

  2. import pandas as pd

  3. import sys, argparse, os

  4. import fileinput

  5. import re

  6. 

  7. # input arguments

  8. parser = argparse.ArgumentParser(description="this script is to get final high confidence calls and information of all replicates")

  9. 

  10. parser.add_argument('-vcfInfo', '--vcfInfo', type=str, help='The txt file of variants information, this file is named as prefix__variant_quality_location.txt',  required=True)

  11. parser.add_argument('-mendelianInfo', '--mendelianInfo', type=str, help='The merged mendelian information of all samples',  required=True)

  12. parser.add_argument('-sample', '--sample_name', type=str, help='which sample of quartet',  required=True)

  13. 

  14. 

  15. args = parser.parse_args()

  16. vcfInfo = args.vcfInfo

  17. mendelianInfo = args.mendelianInfo

  18. sample_name = args.sample_name

  19. 

  20. 

  21. #GT:TWINS:TRIO5:TRIO6:DP:AF:GQ:QD:MQ:FS:QUAL

  22. 

  23. vcf_header = '''##fileformat=VCFv4.2

  24. ##fileDate=20200331

  25. ##source=high_confidence_calls_intergration(choppy app)

  26. ##reference=GRCh38.d1.vd1

  27. #FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">

  28. #FORMAT=<ID=TWINS,Number=1,Type=String,Description="1 for sister consistent, 0 for sister different">

  29. #FORMAT=<ID=TRIO5,Number=1,Type=String,Description="1 for LCL7, LCL8 and LCL5 mendelian consistent, 0 for family violation">

  30. #FORMAT=<ID=TRIO6,Number=1,Type=String,Description="1 for LCL7, LCL8 and LCL6 mendelian consistent, 0 for family violation">

  31. ##FORMAT=<ID=DP,Number=1,Type=Int,Description="Depth">

  32. ##FORMAT=<ID=AF,Number=1,Type=Float,Description="Allele frequency">

  33. ##FORMAT=<ID=GQ,Number=1,Type=Float,Description="Genotype quality">

  34. ##FORMAT=<ID=QD,Number=1,Type=Float,Description="Variant Confidence/Quality by Depth">

  35. ##FORMAT=<ID=MQ,Number=1,Type=Float,Description="Mapping quality">

  36. ##FORMAT=<ID=FS,Number=1,Type=Float,Description="Phred-scaled p-value using Fisher's exact test to detect strand bia">

  37. ##FORMAT=<ID=QUAL,Number=1,Type=Float,Description="variant quality">

  38. ##contig=<ID=chr1,length=248956422>

  39. ##contig=<ID=chr2,length=242193529>

  40. ##contig=<ID=chr3,length=198295559>

  41. ##contig=<ID=chr4,length=190214555>

  42. ##contig=<ID=chr5,length=181538259>

  43. ##contig=<ID=chr6,length=170805979>

  44. ##contig=<ID=chr7,length=159345973>

  45. ##contig=<ID=chr8,length=145138636>

  46. ##contig=<ID=chr9,length=138394717>

  47. ##contig=<ID=chr10,length=133797422>

  48. ##contig=<ID=chr11,length=135086622>

  49. ##contig=<ID=chr12,length=133275309>

  50. ##contig=<ID=chr13,length=114364328>

  51. ##contig=<ID=chr14,length=107043718>

  52. ##contig=<ID=chr15,length=101991189>

  53. ##contig=<ID=chr16,length=90338345>

  54. ##contig=<ID=chr17,length=83257441>

  55. ##contig=<ID=chr18,length=80373285>

  56. ##contig=<ID=chr19,length=58617616>

  57. ##contig=<ID=chr20,length=64444167>

  58. ##contig=<ID=chr21,length=46709983>

  59. ##contig=<ID=chr22,length=50818468>

  60. ##contig=<ID=chrX,length=156040895>

  61. '''

  62. 

  63. # output file

  64. file_name = sample_name + '_mendelian_vcfInfo.vcf'

  65. outfile = open(file_name,'w')

  66. 

  67. outputcolumn = '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t' +  sample_name + '\n'

  68. outfile.write(vcf_header)

  69. outfile.write(outputcolumn)

  70. 

  71. # input files

  72. vcf_info = pd.read_table(vcfInfo)

  73. mendelian_info = pd.read_table(mendelianInfo)

  74. 

  75. merged_df = pd.merge(vcf_info, mendelian_info,  how='outer', left_on=['#CHROM','POS'], right_on = ['#CHROM','POS'])

  76. merged_df = merged_df.fillna('.')

  77. 

  78. #

  79. def parse_INFO(info):

  80. 	strings = info.strip().split(';')

  81. 	keys = []

  82. 	values = []

  83. 	for i in strings:

  84. 		kv = i.split('=')

  85. 		if kv[0] == 'DB':

  86. 			keys.append('DB')

  87. 			values.append('1')

  88. 		else:

  89. 			keys.append(kv[0])

  90. 			values.append(kv[1])

  91. 	infoDict = dict(zip(keys, values))

  92. 	return infoDict

  93. #

  94. for row in merged_df.itertuples():

  95. 	if row.Quartet_DNA_BGI_SEQ2000_BGI_1_20180518_LCL5 != '.':

  96. 		# format

  97. 		# GT:TWINS:TRIO5:TRIO6:DP:AF:GQ:QD:MQ:FS:QUAL

  98. 		FORMAT_x = row.Quartet_DNA_BGI_SEQ2000_BGI_LCL5_1_20180518.split(':')

  99. 		ALT = int(FORMAT_x[1].split(',')[1])

  100. 		if int(FORMAT_x[2]) != 0:

  101. 			AF = round(ALT/int(FORMAT_x[2]),2)

  102. 		else:

  103. 			AF = '.'

  104. 		INFO_x = parse_INFO(row.INFO_x)

  105. 		if FORMAT_x[2] == '0':

  106. 			INFO_x['QD'] = '.'

  107. 		else:

  108. 			pass

  109. 		FORMAT = row.Quartet_DNA_BGI_SEQ2000_BGI_1_20180518_LCL5 + ':' + FORMAT_x[2] + ':' + str(AF) + ':' + FORMAT_x[3] + ':' + INFO_x['QD'] + ':' + INFO_x['MQ'] + ':' + INFO_x['FS'] + ':' + str(row.QUAL_x)

  110. 		# outline

  111. 		outline = row._1 + '\t' + str(row.POS) + '\t' + row.ID_x + '\t' + row.REF_y + '\t' + row.ALT_y + '\t' + '.' + '\t' + '.' + '\t' + '.' + '\t' + 'GT:TWINS:TRIO5:TRIO6:DP:AF:GQ:QD:MQ:FS:QUAL' + '\t' + FORMAT + '\n'

  112. 	else:

  113. 		FORMAT_x = row.Quartet_DNA_BGI_SEQ2000_BGI_LCL5_1_20180518.split(':')

  114. 		ALT = int(FORMAT_x[1].split(',')[1])

  115. 		if int(FORMAT_x[2]) != 0:

  116. 			AF = round(ALT/int(FORMAT_x[2]),2)

  117. 		else:

  118. 			AF = '.'

  119. 		INFO_x = parse_INFO(row.INFO_x)

  120. 		if FORMAT_x[2] == '0':

  121. 			INFO_x['QD'] = '.'

  122. 		else:

  123. 			pass

  124. 		FORMAT = '.:.:.:.' + ':' + FORMAT_x[2] + ':' + str(AF) + ':' + FORMAT_x[3] + ':' + INFO_x['QD'] + ':' + INFO_x['MQ'] + ':' + INFO_x['FS'] + ':' + str(row.QUAL_x)

  125. 		# outline

  126. 		outline = row._1 + '\t' + str(row.POS) + '\t' + row.ID_x + '\t' + row.REF_y + '\t' + row.ALT_y + '\t' + '.' + '\t' + '.' + '\t' + '.' + '\t' + 'GT:TWINS:TRIO5:TRIO6:DP:AF:GQ:QD:MQ:FS:QUAL' + '\t' + FORMAT + '\n'

  127. 	outfile.write(outline)