renluyao
/
high_confidence_calls_manuscript
Sledovat 1
Oblíbit 0
Rozštěpit 0
Zdrojový kód Úkoly 0 Požadavky na natažení 0 Vydání 0 Wiki Aktivita
Nelze vybrat více než 25 témat Téma musí začínat písmenem nebo číslem, může obsahovat pomlčky („-“) a může být dlouhé až 35 znaků.
2 Revize
1 Větev
Strom: 356448c941
Větve Značky
${ item.name }
Vytvořit větev ${ searchTerm }
z „356448c941“
${ noResults }
high_confidence_calls_manus.../codescripts/high_confidence_call_vote.py

256 lines
10KB
Surový Blame Historie 

  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=DPCT,Number=1,Type=Float,Description="Percentage of detected votes">

  29. ##INFO=<ID=VPCT,Number=1,Type=Float,Description="Percentage of consnesus votes">

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  54. '''

  55. 

  56. # read in duplication list

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

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

  59. 

  60. # output file

  61. file_name = prefix + '_annotated.vcf'

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

  63. 

  64. # write VCF

  65. 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' + '\t' + sample_name + '_consensus_alt_seq' +'\n'

  66. outfile.write(vcf_header)

  67. outfile.write(outputcolumn)

  68. 

  69. #function

  70. 

  71. def detected_percentage(strings):

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

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

  74. 	percentage = round((33 - gt.count('.'))/33,4)

  75. 	return(str(percentage))

  76. 

  77. def vote_percentage(strings,consensus_call):

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

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

  80. 	gt = list(map(gt_uniform,[i for i in gt]))

  81. 	percentage = round(gt.count(consensus_call)/33,4)

  82. 	return(str(percentage))

  83. 

  84. def family_vote(strings,consensus_call):

  85. 	pass

  86. 

  87. def gt_uniform(strings):

  88. 	uniformed_gt = ''

  89. 	allele1 = strings.split('/')[0]

  90. 	allele2 = strings.split('/')[1]

  91. 	if int(allele1) > int(allele2):

  92. 		uniformed_gt = allele2 + '/' + allele1

  93. 	else:

  94. 		uniformed_gt = allele1 + '/' + allele2

  95. 	return uniformed_gt

  96. 

  97. def decide_by_rep(strings):

  98. 	consensus_rep = ''

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

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

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

  102. 	# modified gt turn 2/1 to 1/2

  103. 	gt = list(map(gt_uniform,[i for i in gt]))

  104. 	# mendelian consistent?

  105. 	mendelian_dict = Counter(mendelian)

  106. 	highest_mendelian = mendelian_dict.most_common(1)

  107. 	candidate_mendelian = highest_mendelian[0][0]

  108. 	freq_mendelian = highest_mendelian[0][1]

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

  110. 		gt_num_dict = Counter(gt)

  111. 		highest_gt = gt_num_dict.most_common(1)

  112. 		candidate_gt = highest_gt[0][0]

  113. 		freq_gt = highest_gt[0][1]

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

  115. 			consensus_rep = candidate_gt

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

  117. 			consensus_rep = '0/0'

  118. 		else:

  119. 			consensus_rep = 'inconGT'

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

  121. 		consensus_rep = 'noInfo'

  122. 	else:

  123. 		consensus_rep = 'inconMen'

  124. 	return consensus_rep

  125. 

  126. 

  127. def main():

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

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

  130. 		if headline is not None:

  131. 			pass

  132. 		else:

  133. 			line = line.strip()

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

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

  136. 			# check if the variants location is duplicated

  137. 			if variant_id in var_dup:

  138. 				strings[6] = 'dupVar'

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

  140. 				outfile.write(outLine)

  141. 			else:

  142. 				# pre-define

  143. 				pcr_consensus = ''

  144. 				pcr_free_consensus = ''

  145. 				consensus_call = ''

  146. 				consensus_alt_seq = '.'

  147. 				# pcr 

  148. 				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)

  149. 				SEQ2000 = decide_by_rep(pcr[0:3])

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

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

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

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

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

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

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

  157. 				sequence_dict = Counter(sequence_site)

  158. 				highest_sequence = sequence_dict.most_common(1)

  159. 				candidate_sequence = highest_sequence[0][0]

  160. 				freq_sequence = highest_sequence[0][1]

  161. 				if freq_sequence > 4:

  162. 					pcr_consensus = candidate_sequence

  163. 				else:

  164. 					pcr_consensus = 'inconSequenceSite'

  165. 				# pcr-free

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

  167. 				#SEQ2000 = decide_by_rep(pcr_free[0])

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

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

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

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

  172. 				highest_sequence = sequence_dict.most_common(1)

  173. 				candidate_sequence = highest_sequence[0][0]

  174. 				freq_sequence = highest_sequence[0][1]

  175. 				if freq_sequence > 2:

  176. 					pcr_free_consensus = candidate_sequence

  177. 				else:

  178. 					pcr_free_consensus = 'inconSequenceSite'

  179. 				# pcr and pcr-free

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

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

  182. 					consensus_call = pcr_consensus

  183. 					VPCT = vote_percentage(strings[9:],consensus_call)

  184. 					strings[7] = 'VPCT=' + VPCT

  185. 					DPCT = detected_percentage(strings[9:])

  186. 					strings[7] = strings[7] + ';DPCT=' + DPCT

  187. 					# Delete multiple alternative genotype to necessary expression

  188. 					strings[6] = 'reproducible'

  189. 					alt = strings[4]

  190. 					alt_gt = alt.split(',')

  191. 					if len(alt_gt) > 1:

  192. 						allele1 = consensus_call.split('/')[0]

  193. 						allele2 = consensus_call.split('/')[1]

  194. 						if allele1 == '0':

  195. 							allele2_seq = alt_gt[int(allele2) - 1]

  196. 							consensus_alt_seq = allele2_seq

  197. 							consensus_call = '0/1'

  198. 						else:

  199. 							allele1_seq = alt_gt[int(allele1) - 1]

  200. 							allele2_seq = alt_gt[int(allele2) - 1]

  201. 							if int(allele1) > int(allele2):

  202. 								consensus_alt_seq = allele2_seq + ',' + allele1_seq

  203. 								consensus_call = '1/2'

  204. 							elif int(allele1) < int(allele2):

  205. 								consensus_alt_seq = allele1_seq + ',' + allele2_seq

  206. 								consensus_call = '1/2'

  207. 							else:

  208. 								consensus_alt_seq = allele1_seq 

  209. 								consensus_call = '1/1'

  210. 					else:

  211. 						consensus_alt_seq = alt

  212. 				elif (pcr_consensus in tag) and (pcr_free_consensus in tag):

  213. 					consensus_call = 'filtered'

  214. 					strings[6] = 'filtered'

  215. 					DPCT = detected_percentage(strings[9:])

  216. 					strings[7] = 'DPCT=' + DPCT

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

  218. 					consensus_call = 'pcr-free-speicifc'

  219. 					strings[6] = 'pcr-free-speicifc'

  220. 					DPCT = detected_percentage(strings[9:])

  221. 					strings[7] = 'DPCT=' + DPCT

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

  223. 					consensus_call = 'pcr-speicifc'

  224. 					strings[6] = 'pcr-speicifc'

  225. 					DPCT = detected_percentage(strings[9:])

  226. 					strings[7] = 'DPCT=' + DPCT

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

  228. 					consensus_call = 'confirm for parents'

  229. 					strings[6] = 'confirm for parents'					

  230. 					DPCT = detected_percentage(strings[9:])

  231. 					strings[7] = 'DPCT=' + DPCT

  232. 				else:

  233. 					consensus_call = 'filtered'

  234. 					strings[6] = 'filtered'

  235. 					DPCT = detected_percentage(strings[9:])

  236. 					strings[7] = 'DPCT=' + DPCT

  237. 				# output 

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

  239. 				outfile.write(outLine)

  240. 

  241. 

  242. if __name__ == '__main__':

  243. 	main()

  244. 

  245. 

  246. 

  247. 

  248. 

  249. 

  250. 

  251. 

  252. 

  253. 

  254.