renluyao
/
high_confidence_calls_intergration
關注 1
收藏 0
複製 0
程式碼 問題 0 版本發佈 0 Wiki 活動
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.
13 提交
1 分支
目錄樹: f9bbc8520e
分支列表 標籤
${ item.name }
建立分支 ${ searchTerm }
從 'f9bbc8520e'
${ noResults }
high_confidence_calls_inter.../codescripts/high_confidence_call_vote.py

198 line
8.2KB
原始文件 Blame 歷史記錄 

  1. # import modules

  2. import sys, argparse, os

  3. import fileinput

  4. import re

  5. import pandas as pd

  6. from operator import itemgetter

  7. from collections import Counter

  8. from itertools import islice  

  9. 

  10. # input arguments

  11. parser = argparse.ArgumentParser(description="this script is to count voting number")

  12. 

  13. parser.add_argument('-vcf', '--multi_sample_vcf', type=str, help='The VCF file you want to count the voting number',  required=True)

  14. parser.add_argument('-dup', '--dup_list', type=str, help='Duplication list',  required=True)

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

  16. parser.add_argument('-prefix', '--prefix', type=str, help='Prefix of output file name',  required=True)

  17. 

  18. args = parser.parse_args()

  19. multi_sample_vcf = args.multi_sample_vcf

  20. dup_list = args.dup_list

  21. prefix = args.prefix

  22. sample_name = args.sample_name

  23. 

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

  25. ##fileDate=20191224

  26. ##source=high_confidence_calls_intergration(choppy app)

  27. ##reference=GRCh38.d1.vd1

  28. ##INFO=<ID=PCT,Number=1,Type=Float,Description="Percentage of votes">

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

  30. ##FORMAT=<ID=TWINS,Number=0,Type=Flag,Description="0 for sister consistent, 1 for sister inconsistent">

  31. ##FORMAT=<ID=TRIO5,Number=0,Type=Flag,Description="0 for trio consistent, 1 for trio inconsistent">

  32. ##FORMAT=<ID=TRIO6,Number=0,Type=Flag,Description="0 for trio consistent, 1 for trio inconsistent">

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  56. '''

  57. 

  58. # read in duplication list

  59. dup = pd.read_table(dup_list,header=None)

  60. var_dup = dup[0].tolist()

  61. 

  62. # output file

  63. file_name = prefix + '_annotated.vcf'

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

  65. 

  66. # write VCF

  67. outputcolumn = '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tQuartet_DNA_BGI_SEQ2000_BGI_1_20180518_LCL5\tQuartet_DNA_BGI_SEQ2000_BGI_2_20180530_LCL5\tQuartet_DNA_BGI_SEQ2000_BGI_3_20180530_LCL5\tQuartet_DNA_BGI_SEQ2000_WGE_1_20190402_LCL5\tQuartet_DNA_BGI_SEQ2000_WGE_2_20190402_LCL5\tQuartet_DNA_BGI_SEQ500_BGI_1_20180328_LCL5 \tQuartet_DNA_BGI_SEQ500_BGI_2_20180328_LCL5\tQuartet_DNA_BGI_SEQ500_BGI_3_20180328_LCL5\tQuartet_DNA_ILM_Nova_ARD_1_20181108_LCL5\tQuartet_DNA_ILM_Nova_ARD_2_20181108_LCL5\tQuartet_DNA_ILM_Nova_ARD_3_20181108_LCL5\tQuartet_DNA_ILM_Nova_ARD_4_20190111_LCL5\tQuartet_DNA_ILM_Nova_ARD_5_20190111_LCL5\tQuartet_DNA_ILM_Nova_ARD_6_20190111_LCL5\tQuartet_DNA_ILM_Nova_BRG_1_20171024_LCL5\tQuartet_DNA_ILM_Nova_BRG_1_20180930_LCL5\tQuartet_DNA_ILM_Nova_BRG_2_20180930_LCL5\tQuartet_DNA_ILM_Nova_BRG_3_20180930_LCL5\tQuartet_DNA_ILM_Nova_GAC_1_20171025_LCL5\tQuartet_DNA_ILM_Nova_NVG_1_20171024_LCL5\tQuartet_DNA_ILM_Nova_WUX_1_20171024_LCL5\tQuartet_DNA_ILM_XTen_ARD_1_20170403_LCL5\tQuartet_DNA_ILM_XTen_ARD_2_20170403_LCL5\tQuartet_DNA_ILM_XTen_ARD_3_20170403_LCL5\tQuartet_DNA_ILM_XTen_NVG_1_20170329_LCL5\tQuartet_DNA_ILM_XTen_NVG_2_20170329_LCL5\tQuartet_DNA_ILM_XTen_NVG_3_20170329_LCL5\tQuartet_DNA_ILM_XTen_WUX_1_20170216_LCL5\tQuartet_DNA_ILM_XTen_WUX_2_20170216_LCL5\tQuartet_DNA_ILM_XTen_WUX_3_20170216_LCL5\tQuartet_DNA_ILM_XTen_WUX_4_20180703_LCL5\tQuartet_DNA_ILM_XTen_WUX_5_20180703_LCL5\tQuartet_DNA_ILM_XTen_WUX_6_20180703_LCL5' +'\t'+ sample_name+'_pcr'+'\t' + sample_name+'_pcr-free'+ '\t'+ sample_name +'consensus' + '\n'

  68. outfile.write(vcf_header)

  69. outfile.write(outputcolumn)

  70. 

  71. #function

  72. def vote_percentage(strings):

  73. 	strings = [x.replace('0/0','.') for x in strings]

  74. 	gt = [x.split(':')[0] for x in strings]

  75. 	percentage = round((33 - gt.count('.'))/33,2)

  76. 	return(str(percentage))

  77. 

  78. def decide_by_rep(strings):

  79. 	consensus_rep = ''

  80. 	mendelian = [x[-5:] for x in strings]

  81. 	strings = [x.replace('.','0/0') for x in strings]

  82. 	gt = [x.split(':')[0] for x in strings]

  83. 	# mendelian consistent?

  84. 	mendelian_dict = Counter(mendelian)

  85. 	highest_mendelian = mendelian_dict.most_common(1)

  86. 	candidate_mendelian = highest_mendelian[0][0]

  87. 	freq_mendelian = highest_mendelian[0][1]

  88. 	if (candidate_mendelian == '1:1:1') and (freq_mendelian >= 2):

  89. 		gt_num_dict = Counter(gt)

  90. 		highest_gt = gt_num_dict.most_common(1)

  91. 		candidate_gt = highest_gt[0][0]

  92. 		freq_gt = highest_gt[0][1]

  93. 		if (candidate_gt != '0/0') and (freq_gt >= 2):

  94. 			consensus_rep = candidate_gt

  95. 		elif (candidate_gt == '0/0') and (freq_gt >= 2):

  96. 			consensus_rep = '0/0'

  97. 		else:

  98. 			consensus_rep = 'inconGT'

  99. 	elif (candidate_mendelian == '.') and (freq_mendelian >= 2):

  100. 		consensus_rep = 'noInfo'

  101. 	else:

  102. 		consensus_rep = 'inconMen'

  103. 	return consensus_rep

  104. 

  105. 

  106. def main():

  107. 	for line in fileinput.input(multi_sample_vcf):

  108. 		headline = re.match('^\#',line)

  109. 		if headline is not None:

  110. 			pass

  111. 		else:

  112. 			line = line.strip()

  113. 			strings = line.split('\t')

  114. 			variant_id = '_'.join([strings[0],strings[1]])

  115. 			# check if the variants location is duplicated

  116. 			if variant_id in var_dup:

  117. 				outLine = '\t'.join(strings) + '\t' + '.' +'\t' + '.' + '\t' + 'dupVar' + '\n'

  118. 				outfile.write(outLine)

  119. 			else:

  120. 				# pre-define

  121. 				pcr_consensus = ''

  122. 				pcr_free_consensus = ''

  123. 				consensus_call = ''

  124. 				# pcr 

  125. 				pcr = itemgetter(*[9,10,11,12,14,15,16,23,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41])(strings)

  126. 				SEQ2000 = decide_by_rep(pcr[0:4])

  127. 				SEQ500 = decide_by_rep(pcr[4:7])

  128. 				Nova = decide_by_rep(pcr[7:11])

  129. 				XTen_ARD = decide_by_rep(pcr[11:14])

  130. 				XTen_NVG = decide_by_rep(pcr[14:17])

  131. 				XTen_WUX_1 = decide_by_rep(pcr[17:20])

  132. 				XTen_WUX_2 = decide_by_rep(pcr[20:23])

  133. 				sequence_site = [SEQ2000,SEQ500,Nova,XTen_ARD,XTen_NVG,XTen_WUX_1,XTen_WUX_2]

  134. 				sequence_dict = Counter(sequence_site)

  135. 				highest_sequence = sequence_dict.most_common(1)

  136. 				candidate_sequence = highest_sequence[0][0]

  137. 				freq_sequence = highest_sequence[0][1]

  138. 				if freq_sequence > 4:

  139. 					pcr_consensus = candidate_sequence

  140. 				else:

  141. 					pcr_consensus = 'inconSequenceSite'

  142. 				# pcr-free

  143. 				pcr_free = itemgetter(*[13,17,18,19,20,21,22,24,25,26])(strings)

  144. 				SEQ2000 = decide_by_rep(pcr_free[0])

  145. 				Nova_ARD_1 = decide_by_rep(pcr_free[1:4])

  146. 				Nova_ARD_2 = decide_by_rep(pcr_free[4:7])

  147. 				Nova_BRG = decide_by_rep(pcr_free[7:10])

  148. 				sequence_site = [SEQ2000,Nova_ARD_1,Nova_ARD_2,Nova_BRG]

  149. 				highest_sequence = sequence_dict.most_common(1)

  150. 				candidate_sequence = highest_sequence[0][0]

  151. 				freq_sequence = highest_sequence[0][1]

  152. 				if freq_sequence > 2:

  153. 					pcr_free_consensus = candidate_sequence

  154. 				else:

  155. 					pcr_free_consensus = 'inconSequenceSite'

  156. 				# pcr and pcr-free

  157. 				tag = ['inconGT','noInfo','inconMen','inconSequenceSite']

  158. 				if (pcr_consensus == pcr_free_consensus) and (pcr_consensus not in tag) and (pcr_consensus != '0/0'):

  159. 					consensus_call = pcr_consensus

  160. 					strings[6] = 'reproducible'

  161. 				elif (pcr_consensus in tag) or (pcr_free_consensus in tag):

  162. 					consensus_call = 'filtered'

  163. 					strings[6] = '.'

  164. 				elif (pcr_consensus == '0/0') and (pcr_free_consensus not in tag) and (pcr_free_consensus != '0/0'):

  165. 					consensus_call = 'pcr-free-speicifc'

  166. 					strings[6] = '.'

  167. 				elif (pcr_consensus != '0/0') and (pcr_consensus not in tag) and (pcr_free_consensus == '0/0'):

  168. 					consensus_call = 'pcr-speicifc'

  169. 					strings[6] = '.'

  170. 				elif (pcr_consensus == '0/0') and (pcr_free_consensus == '0/0'):

  171. 					consensus_call = 'confirm for parents'

  172. 					strings[6] = '.'					

  173. 				else:

  174. 					consensus_call = 'filtered'

  175. 					strings[6] = '.'

  176. 				# percentage

  177. 				percentage = vote_percentage(strings[9:])

  178. 				strings[7] = 'PCT=' + percentage

  179. 				# output 

  180. 				outLine = '\t'.join(strings) + '\t' + pcr_consensus +'\t' + pcr_free_consensus + '\t' + consensus_call + '\n'

  181. 				outfile.write(outLine)

  182. 

  183. 

  184. if __name__ == '__main__':

  185. 	main()

  186. 

  187. 

  188. 

  189. 

  190. 

  191. 

  192. 

  193. 

  194. 

  195. 

  196.