renluyao
/
high_confidence_calls_manuscript
Observar 1
Favorito 0
Fork 0
Código Issues 0 Pull requests 0 Versões 0 Wiki Atividade
Você não pode selecionar mais de 25 tópicos Os tópicos devem começar com uma letra ou um número, podem incluir traços ('-') e podem ter até 35 caracteres.
13 Commits
1 Branch
Tag: b3dfa2e638
Branches Tags
${ item.name }
Criar branch ${ searchTerm }
de b3dfa2e638
${ noResults }
high_confidence_calls_manus.../codescripts/merge_mendelian_vcfinfo.py

130 linhas
5.0KB
Original Anotar Histórico 

  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=Flag,Description="1 for sister consistent, 0 for sister different">

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

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

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

  32. ##FORMAT=<ID=ALT,Number=1,Type=Integer,Description="Alternative Depth">

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  62. '''

  63. 

  64. # output file

  65. file_name = sample_name + '_mendelian_vcfInfo.vcf'

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

  67. 

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

  69. outfile.write(vcf_header)

  70. outfile.write(outputcolumn)

  71. 

  72. # input files

  73. vcf_info = pd.read_table(vcfInfo)

  74. mendelian_info = pd.read_table(mendelianInfo)

  75. 

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

  77. merged_df = merged_df.fillna('.')

  78. 

  79. #

  80. def parse_INFO(info):

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

  82. 	keys = []

  83. 	values = []

  84. 	for i in strings:

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

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

  87. 			keys.append('DB')

  88. 			values.append('1')

  89. 		else:

  90. 			keys.append(kv[0])

  91. 			values.append(kv[1])

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

  93. 	return infoDict

  94. #

  95. for row in merged_df.itertuples():

  96. 	if row[18] != '.':

  97. 		# format

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

  99. 		FORMAT_x = row[10].split(':')

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

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

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

  103. 		else:

  104. 			AF = '.'

  105. 		INFO_x = parse_INFO(row.INFO_x)

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

  107. 			INFO_x['QD'] = '.'

  108. 		else:

  109. 			pass

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

  111. 		# outline

  112. 		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:ALT:AF:GQ:QD:MQ:FS:QUAL' + '\t' + FORMAT + '\n'

  113. 	else:

  114. 		rawGT = row[10].split(':')

  115. 		FORMAT_x = row[10].split(':')

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

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

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

  119. 		else:

  120. 			AF = '.'

  121. 		INFO_x = parse_INFO(row.INFO_x)

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

  123. 			INFO_x['QD'] = '.'

  124. 		else:

  125. 			pass

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

  127. 		# outline

  128. 		outline = row._1 + '\t' + str(row.POS) + '\t' + row.ID_x + '\t' + row.REF_x + '\t' + row.ALT_x + '\t' + '.' + '\t' + '.' + '\t' + '.' + '\t' + 'GT:TWINS:TRIO5:TRIO6:DP:ALT:AF:GQ:QD:MQ:FS:QUAL:rawGT' + '\t' + FORMAT + '\n'

  129. 	outfile.write(outline)