renluyao
/
quartet-wgs-germline-data-generation-qc
Segui 1
Vota 0
Forka 0
Codice Problemi 0 Pull Requests 0 Rilasci 2 Wiki Attività
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.
15 Commit
1 Ramo (Branch)
Albero (Tree): 211eb24ff3
Rami (Branch) Tag
${ item.name }
Crea branch ${ searchTerm }
da '211eb24ff3'
${ noResults }
quartet-wgs-germline-data-g.../codescripts/voted_by_vcfinfo_mendeliani...

296 lines
18KB
Originale Blame Cronologia 

  1. from __future__ import division

  2. from glob import glob

  3. import sys, argparse, os

  4. import fileinput

  5. import re

  6. import pandas as pd

  7. from operator import itemgetter

  8. from collections import Counter

  9. from itertools import islice

  10. from numpy import *

  11. import statistics

  12. 

  13. # input arguments

  14. parser = argparse.ArgumentParser(description="this script is to merge mendelian and vcfinfo, and extract high_confidence_calls")

  15. 

  16. parser.add_argument('-folder', '--folder', type=str, help='directory that holds all the mendelian info',  required=True)

  17. parser.add_argument('-vcf', '--vcf', type=str, help='merged multiple sample vcf',  required=True)

  18. 

  19. 

  20. args = parser.parse_args()

  21. folder = args.folder

  22. vcf = args.vcf

  23. 

  24. # input files

  25. folder = folder + '/*.txt'

  26. filenames = glob(folder)

  27. dataframes = []

  28. for filename in filenames:

  29. 	dataframes.append(pd.read_table(filename,header=None))

  30. 

  31. dfs = [df.set_index([0, 1]) for df in dataframes]

  32. merged_mendelian = pd.concat(dfs, axis=1).reset_index()

  33. family_name = [i.split('/')[-1].replace('.txt','') for i in filenames]

  34. columns = ['CHROM','POS'] + family_name

  35. merged_mendelian.columns = columns

  36. 

  37. vcf_dat = pd.read_table(vcf)

  38. 

  39. merged_df = pd.merge(merged_mendelian, vcf_dat,  how='outer', left_on=['CHROM','POS'], right_on = ['#CHROM','POS'])

  40. merged_df = merged_df.fillna('nan')

  41. 

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

  43. ##fileDate=20200501

  44. ##source=high_confidence_calls_intergration(choppy app)

  45. ##reference=GRCh38.d1.vd1

  46. ##INFO=<ID=VOTED,Number=1,Type=Integer,Description="Number mendelian consisitent votes">

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

  48. ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Sum depth of all samples">

  49. ##FORMAT=<ID=ALT,Number=1,Type=Integer,Description="Sum alternative depth of all samples">

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  73. '''

  74. # output files

  75. benchmark_LCL5 = open('LCL5_voted.vcf','w')

  76. benchmark_LCL6 = open('LCL6_voted.vcf','w')

  77. benchmark_LCL7 = open('LCL7_voted.vcf','w')

  78. benchmark_LCL8 = open('LCL8_voted.vcf','w')

  79. 

  80. all_sample_outfile = open('all_sample_information.txt','w')

  81. 

  82. # write VCF

  83. LCL5_col = '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tLCL5_benchmark_calls\n'

  84. LCL6_col = '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tLCL6_benchmark_calls\n'

  85. LCL7_col = '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tLCL7_benchmark_calls\n'

  86. LCL8_col = '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tLCL8_benchmark_calls\n'

  87. 

  88. benchmark_LCL5.write(vcf_header)

  89. benchmark_LCL5.write(LCL5_col)

  90. benchmark_LCL6.write(vcf_header)

  91. benchmark_LCL6.write(LCL6_col)

  92. benchmark_LCL7.write(vcf_header)

  93. benchmark_LCL7.write(LCL7_col)

  94. benchmark_LCL8.write(vcf_header)

  95. benchmark_LCL8.write(LCL8_col)

  96. 

  97. # all info

  98. all_info_col = 'CHROM\tPOS\tLCL5_pcr_consensus\tLCL5_pcr_free_consensus\tLCL5_mendelian_num\tLCL5_consensus_call\tLCL5_consensus_alt_seq\tLCL5_alt\tLCL5_dp\tLCL5_detected_num\tLCL6_pcr_consensus\tLCL6_pcr_free_consensus\tLCL6_mendelian_num\tLCL6_consensus_call\tLCL6_consensus_alt_seq\tLCL6_alt\tLCL6_dp\tLCL6_detected_num\tLCL7_pcr_consensus\tLCL7_pcr_free_consensus\tLCL7_mendelian_num\t LCL7_consensus_call\tLCL7_consensus_alt_seq\tLCL7_alt\tLCL7_dp\tLCL7_detected_num\tLCL8_pcr_consensus\tLCL8_pcr_free_consensus\tLCL8_mendelian_num\tLCL8_consensus_call\tLCL8_consensus_alt_seq\tLCL8_alt\tLCL8_dp\tLCL8_detected_num\n'

  99. all_sample_outfile.write(all_info_col)

  100. 

  101. # function

  102. def decide_by_rep(vcf_list,mendelian_list):

  103. 	consensus_rep = ''

  104. 	gt = [x.split(':')[0] for x in vcf_list]

  105. 	# mendelian consistent?

  106. 	mendelian_dict = Counter(mendelian_list)

  107. 	highest_mendelian = mendelian_dict.most_common(1)

  108. 	candidate_mendelian = highest_mendelian[0][0]

  109. 	freq_mendelian = highest_mendelian[0][1]

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

  111. 		con_loc = [i for i in range(len(mendelian_list)) if mendelian_list[i] == '1:1:1']

  112. 		gt_con = itemgetter(*con_loc)(gt)

  113. 		gt_num_dict = Counter(gt_con)

  114. 		highest_gt = gt_num_dict.most_common(1)

  115. 		candidate_gt = highest_gt[0][0]

  116. 		freq_gt = highest_gt[0][1]

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

  118. 			consensus_rep = candidate_gt

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

  120. 			consensus_rep = 'noGTInfo'

  121. 		else:

  122. 			consensus_rep = 'inconGT'

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

  124. 		consensus_rep = 'noMenInfo'

  125. 	else:

  126. 		consensus_rep = 'inconMen'

  127. 	return consensus_rep

  128. 

  129. def consensus_call(vcf_info_list,mendelian_list,alt_seq):

  130. 	pcr_consensus = '.'

  131. 	pcr_free_consensus = '.'

  132. 	mendelian_num = '.'

  133. 	consensus_call = '.'

  134. 	consensus_alt_seq = '.'

  135. 	# pcr

  136. 	SEQ2000 = decide_by_rep(vcf_info_list[0:3],mendelian_list[0:3])

  137. 	XTen_ARD = decide_by_rep(vcf_info_list[18:21],mendelian_list[18:21])

  138. 	XTen_NVG = decide_by_rep(vcf_info_list[21:24],mendelian_list[21:24])

  139. 	XTen_WUX = decide_by_rep(vcf_info_list[24:27],mendelian_list[24:27])

  140. 	pcr_sequence_site = [SEQ2000,XTen_ARD,XTen_NVG,XTen_WUX]

  141. 	pcr_sequence_dict = Counter(pcr_sequence_site)

  142. 	pcr_highest_sequence = pcr_sequence_dict.most_common(1)

  143. 	pcr_candidate_sequence = pcr_highest_sequence[0][0]

  144. 	pcr_freq_sequence = pcr_highest_sequence[0][1]

  145. 	if pcr_freq_sequence > 2:

  146. 		pcr_consensus = pcr_candidate_sequence

  147. 	else:

  148. 		pcr_consensus = 'inconSequenceSite'

  149. 	# pcr-free

  150. 	T7_WGE = decide_by_rep(vcf_info_list[3:6],mendelian_list[3:6])

  151. 	Nova_ARD_1 = decide_by_rep(vcf_info_list[6:9],mendelian_list[6:9])

  152. 	Nova_ARD_2 = decide_by_rep(vcf_info_list[9:12],mendelian_list[9:12])

  153. 	Nova_BRG = decide_by_rep(vcf_info_list[12:15],mendelian_list[12:15])

  154. 	Nova_WUX = decide_by_rep(vcf_info_list[15:18],mendelian_list[15:18])

  155. 	sequence_site = [T7_WGE,Nova_ARD_1,Nova_ARD_2,Nova_BRG,Nova_WUX]

  156. 	sequence_dict = Counter(sequence_site)

  157. 	highest_sequence = sequence_dict.most_common(1)

  158. 	candidate_sequence = highest_sequence[0][0]

  159. 	freq_sequence = highest_sequence[0][1]

  160. 	if freq_sequence > 3:

  161. 		pcr_free_consensus = candidate_sequence

  162. 	else:

  163. 		pcr_free_consensus = 'inconSequenceSite'

  164. 	# net alt, dp

  165. 	# alt

  166. 	AD = [x.split(':')[1] for x in vcf_info_list]

  167. 	ALT = [x.split(',')[1] for x in AD]

  168. 	ALT = [int(x) for x in ALT]

  169. 	ALL_ALT = sum(ALT)

  170. 	# dp

  171. 	DP = [x.split(':')[2] for x in vcf_info_list]

  172. 	DP = [int(x) for x in DP]

  173. 	ALL_DP = sum(DP)

  174. 	# detected number

  175. 	gt = [x.split(':')[0] for x in vcf_info_list]

  176. 	gt = [x.replace('0/0','.') for x in gt]

  177. 	gt = [x.replace('./.','.') for x in gt]

  178. 	detected_num = 27 - gt.count('.')

  179. 	# decide consensus calls

  180. 	tag = ['inconGT','noMenInfo','inconMen','inconSequenceSite','noGTInfo']

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

  182. 		consensus_call = pcr_consensus

  183. 		gt = [x.split(':')[0] for x in vcf_info_list]

  184. 		indices = [i for i, x in enumerate(gt) if x == consensus_call]

  185. 		matched_mendelian = itemgetter(*indices)(mendelian_list)

  186. 		mendelian_num = matched_mendelian.count('1:1:1')

  187. 		# Delete multiple alternative genotype to necessary expression

  188. 		alt_gt = alt_seq.split(',')

  189. 		if len(alt_gt) > 1:

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

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

  192. 			if allele1 == '0':

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

  194. 				consensus_alt_seq = allele2_seq

  195. 				consensus_call = '0/1'

  196. 			else:

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

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

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

  200. 					consensus_alt_seq = allele2_seq + ',' + allele1_seq

  201. 					consensus_call = '1/2'

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

  203. 					consensus_alt_seq = allele1_seq + ',' + allele2_seq

  204. 					consensus_call = '1/2'

  205. 				else:

  206. 					consensus_alt_seq = allele1_seq 

  207. 					consensus_call = '1/1'

  208. 		else:

  209. 			consensus_alt_seq = alt_seq

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

  211. 		consensus_call = 'filtered'				

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

  213. 		consensus_call = 'pcr-free-speicifc'				

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

  215. 		consensus_call = 'pcr-speicifc'			

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

  217. 		consensus_call = '0/0'								

  218. 	else:

  219. 		consensus_call = 'filtered'

  220. 	return pcr_consensus, pcr_free_consensus, mendelian_num, consensus_call, consensus_alt_seq, ALL_ALT, ALL_DP, detected_num

  221. 

  222. 

  223. for row in merged_df.itertuples():

  224. 	mendelian_list = [row.Quartet_DNA_BGI_SEQ2000_BGI_1_20180518,row.Quartet_DNA_BGI_SEQ2000_BGI_2_20180530,row.Quartet_DNA_BGI_SEQ2000_BGI_3_20180530, \

  225. 					row.Quartet_DNA_BGI_T7_WGE_1_20191105,row.Quartet_DNA_BGI_T7_WGE_2_20191105,row.Quartet_DNA_BGI_T7_WGE_3_20191105, \

  226. 					row.Quartet_DNA_ILM_Nova_ARD_1_20181108,row.Quartet_DNA_ILM_Nova_ARD_2_20181108,row.Quartet_DNA_ILM_Nova_ARD_3_20181108, \

  227. 					row.Quartet_DNA_ILM_Nova_ARD_4_20190111,row.Quartet_DNA_ILM_Nova_ARD_5_20190111,row.Quartet_DNA_ILM_Nova_ARD_6_20190111, \

  228. 					row.Quartet_DNA_ILM_Nova_BRG_1_20180930,row.Quartet_DNA_ILM_Nova_BRG_2_20180930,row.Quartet_DNA_ILM_Nova_BRG_3_20180930, \

  229. 					row.Quartet_DNA_ILM_Nova_WUX_1_20190917,row.Quartet_DNA_ILM_Nova_WUX_2_20190917,row.Quartet_DNA_ILM_Nova_WUX_3_20190917, \

  230. 					row.Quartet_DNA_ILM_XTen_ARD_1_20170403,row.Quartet_DNA_ILM_XTen_ARD_2_20170403,row.Quartet_DNA_ILM_XTen_ARD_3_20170403, \

  231. 					row.Quartet_DNA_ILM_XTen_NVG_1_20170329,row.Quartet_DNA_ILM_XTen_NVG_2_20170329,row.Quartet_DNA_ILM_XTen_NVG_3_20170329, \

  232. 					row.Quartet_DNA_ILM_XTen_WUX_1_20170216,row.Quartet_DNA_ILM_XTen_WUX_2_20170216,row.Quartet_DNA_ILM_XTen_WUX_3_20170216]

  233. 	lcl5_list = [row.Quartet_DNA_BGI_SEQ2000_BGI_LCL5_1_20180518,row.Quartet_DNA_BGI_SEQ2000_BGI_LCL5_2_20180530,row.Quartet_DNA_BGI_SEQ2000_BGI_LCL5_3_20180530, \

  234. 					row.Quartet_DNA_BGI_T7_WGE_LCL5_1_20191105,row.Quartet_DNA_BGI_T7_WGE_LCL5_2_20191105,row.Quartet_DNA_BGI_T7_WGE_LCL5_3_20191105, \

  235. 					row.Quartet_DNA_ILM_Nova_ARD_LCL5_1_20181108,row.Quartet_DNA_ILM_Nova_ARD_LCL5_2_20181108,row.Quartet_DNA_ILM_Nova_ARD_LCL5_3_20181108, \

  236. 					row.Quartet_DNA_ILM_Nova_ARD_LCL5_4_20190111,row.Quartet_DNA_ILM_Nova_ARD_LCL5_5_20190111,row.Quartet_DNA_ILM_Nova_ARD_LCL5_6_20190111, \

  237. 					row.Quartet_DNA_ILM_Nova_BRG_LCL5_1_20180930,row.Quartet_DNA_ILM_Nova_BRG_LCL5_2_20180930,row.Quartet_DNA_ILM_Nova_BRG_LCL5_3_20180930, \

  238. 					row.Quartet_DNA_ILM_Nova_WUX_LCL5_1_20190917,row.Quartet_DNA_ILM_Nova_WUX_LCL5_2_20190917,row.Quartet_DNA_ILM_Nova_WUX_LCL5_3_20190917, \

  239. 					row.Quartet_DNA_ILM_XTen_ARD_LCL5_1_20170403,row.Quartet_DNA_ILM_XTen_ARD_LCL5_2_20170403,row.Quartet_DNA_ILM_XTen_ARD_LCL5_3_20170403, \

  240. 					row.Quartet_DNA_ILM_XTen_NVG_LCL5_1_20170329,row.Quartet_DNA_ILM_XTen_NVG_LCL5_2_20170329,row.Quartet_DNA_ILM_XTen_NVG_LCL5_3_20170329, \

  241. 					row.Quartet_DNA_ILM_XTen_WUX_LCL5_1_20170216,row.Quartet_DNA_ILM_XTen_WUX_LCL5_2_20170216,row.Quartet_DNA_ILM_XTen_WUX_LCL5_3_20170216]

  242. 	lcl6_list = [row.Quartet_DNA_BGI_SEQ2000_BGI_LCL6_1_20180518,row.Quartet_DNA_BGI_SEQ2000_BGI_LCL6_2_20180530,row.Quartet_DNA_BGI_SEQ2000_BGI_LCL6_3_20180530, \

  243. 					row.Quartet_DNA_BGI_T7_WGE_LCL6_1_20191105,row.Quartet_DNA_BGI_T7_WGE_LCL6_2_20191105,row.Quartet_DNA_BGI_T7_WGE_LCL6_3_20191105, \

  244. 					row.Quartet_DNA_ILM_Nova_ARD_LCL6_1_20181108,row.Quartet_DNA_ILM_Nova_ARD_LCL6_2_20181108,row.Quartet_DNA_ILM_Nova_ARD_LCL6_3_20181108, \

  245. 					row.Quartet_DNA_ILM_Nova_ARD_LCL6_4_20190111,row.Quartet_DNA_ILM_Nova_ARD_LCL6_5_20190111,row.Quartet_DNA_ILM_Nova_ARD_LCL6_6_20190111, \

  246. 					row.Quartet_DNA_ILM_Nova_BRG_LCL6_1_20180930,row.Quartet_DNA_ILM_Nova_BRG_LCL6_2_20180930,row.Quartet_DNA_ILM_Nova_BRG_LCL6_3_20180930, \

  247. 					row.Quartet_DNA_ILM_Nova_WUX_LCL6_1_20190917,row.Quartet_DNA_ILM_Nova_WUX_LCL6_2_20190917,row.Quartet_DNA_ILM_Nova_WUX_LCL6_3_20190917, \

  248. 					row.Quartet_DNA_ILM_XTen_ARD_LCL6_1_20170403,row.Quartet_DNA_ILM_XTen_ARD_LCL6_2_20170403,row.Quartet_DNA_ILM_XTen_ARD_LCL6_3_20170403, \

  249. 					row.Quartet_DNA_ILM_XTen_NVG_LCL6_1_20170329,row.Quartet_DNA_ILM_XTen_NVG_LCL6_2_20170329,row.Quartet_DNA_ILM_XTen_NVG_LCL6_3_20170329, \

  250. 					row.Quartet_DNA_ILM_XTen_WUX_LCL6_1_20170216,row.Quartet_DNA_ILM_XTen_WUX_LCL6_2_20170216,row.Quartet_DNA_ILM_XTen_WUX_LCL6_3_20170216]

  251. 	lcl7_list = [row.Quartet_DNA_BGI_SEQ2000_BGI_LCL7_1_20180518,row.Quartet_DNA_BGI_SEQ2000_BGI_LCL7_2_20180530,row.Quartet_DNA_BGI_SEQ2000_BGI_LCL7_3_20180530, \

  252. 					row.Quartet_DNA_BGI_T7_WGE_LCL7_1_20191105,row.Quartet_DNA_BGI_T7_WGE_LCL7_2_20191105,row.Quartet_DNA_BGI_T7_WGE_LCL7_3_20191105, \

  253. 					row.Quartet_DNA_ILM_Nova_ARD_LCL7_1_20181108,row.Quartet_DNA_ILM_Nova_ARD_LCL7_2_20181108,row.Quartet_DNA_ILM_Nova_ARD_LCL7_3_20181108, \

  254. 					row.Quartet_DNA_ILM_Nova_ARD_LCL7_4_20190111,row.Quartet_DNA_ILM_Nova_ARD_LCL7_5_20190111,row.Quartet_DNA_ILM_Nova_ARD_LCL7_6_20190111, \

  255. 					row.Quartet_DNA_ILM_Nova_BRG_LCL7_1_20180930,row.Quartet_DNA_ILM_Nova_BRG_LCL7_2_20180930,row.Quartet_DNA_ILM_Nova_BRG_LCL7_3_20180930, \

  256. 					row.Quartet_DNA_ILM_Nova_WUX_LCL7_1_20190917,row.Quartet_DNA_ILM_Nova_WUX_LCL7_2_20190917,row.Quartet_DNA_ILM_Nova_WUX_LCL7_3_20190917, \

  257. 					row.Quartet_DNA_ILM_XTen_ARD_LCL7_1_20170403,row.Quartet_DNA_ILM_XTen_ARD_LCL7_2_20170403,row.Quartet_DNA_ILM_XTen_ARD_LCL7_3_20170403, \

  258. 					row.Quartet_DNA_ILM_XTen_NVG_LCL7_1_20170329,row.Quartet_DNA_ILM_XTen_NVG_LCL7_2_20170329,row.Quartet_DNA_ILM_XTen_NVG_LCL7_3_20170329, \

  259. 					row.Quartet_DNA_ILM_XTen_WUX_LCL7_1_20170216,row.Quartet_DNA_ILM_XTen_WUX_LCL7_2_20170216,row.Quartet_DNA_ILM_XTen_WUX_LCL7_3_20170216]

  260. 	lcl8_list = [row.Quartet_DNA_BGI_SEQ2000_BGI_LCL8_1_20180518,row.Quartet_DNA_BGI_SEQ2000_BGI_LCL8_2_20180530,row.Quartet_DNA_BGI_SEQ2000_BGI_LCL8_3_20180530, \

  261. 					row.Quartet_DNA_BGI_T7_WGE_LCL8_1_20191105,row.Quartet_DNA_BGI_T7_WGE_LCL8_2_20191105,row.Quartet_DNA_BGI_T7_WGE_LCL8_3_20191105, \

  262. 					row.Quartet_DNA_ILM_Nova_ARD_LCL8_1_20181108,row.Quartet_DNA_ILM_Nova_ARD_LCL8_2_20181108,row.Quartet_DNA_ILM_Nova_ARD_LCL8_3_20181108, \

  263. 					row.Quartet_DNA_ILM_Nova_ARD_LCL8_4_20190111,row.Quartet_DNA_ILM_Nova_ARD_LCL8_5_20190111,row.Quartet_DNA_ILM_Nova_ARD_LCL8_6_20190111, \

  264. 					row.Quartet_DNA_ILM_Nova_BRG_LCL8_1_20180930,row.Quartet_DNA_ILM_Nova_BRG_LCL8_2_20180930,row.Quartet_DNA_ILM_Nova_BRG_LCL8_3_20180930, \

  265. 					row.Quartet_DNA_ILM_Nova_WUX_LCL8_1_20190917,row.Quartet_DNA_ILM_Nova_WUX_LCL8_2_20190917,row.Quartet_DNA_ILM_Nova_WUX_LCL8_3_20190917, \

  266. 					row.Quartet_DNA_ILM_XTen_ARD_LCL8_1_20170403,row.Quartet_DNA_ILM_XTen_ARD_LCL8_2_20170403,row.Quartet_DNA_ILM_XTen_ARD_LCL8_3_20170403, \

  267. 					row.Quartet_DNA_ILM_XTen_NVG_LCL8_1_20170329,row.Quartet_DNA_ILM_XTen_NVG_LCL8_2_20170329,row.Quartet_DNA_ILM_XTen_NVG_LCL8_3_20170329, \

  268. 					row.Quartet_DNA_ILM_XTen_WUX_LCL8_1_20170216,row.Quartet_DNA_ILM_XTen_WUX_LCL8_2_20170216,row.Quartet_DNA_ILM_XTen_WUX_LCL8_3_20170216]

  269. 	# LCL5

  270. 	LCL5_pcr_consensus, LCL5_pcr_free_consensus, LCL5_mendelian_num, LCL5_consensus_call, LCL5_consensus_alt_seq, LCL5_alt, LCL5_dp, LCL5_detected_num = consensus_call(lcl5_list,mendelian_list,row.ALT)

  271. 	if LCL5_mendelian_num != '.':

  272. 		LCL5_output = row.CHROM + '\t' + str(row.POS) + '\t' + '.' + '\t' + row.REF + '\t' + LCL5_consensus_alt_seq + '\t' + '.' + '\t' + '.' + '\t' +'VOTED=' + str(LCL5_mendelian_num) + '\t' + 'GT:ALT:DP' + '\t' + LCL5_consensus_call + ':' + str(LCL5_alt) + ':' + str(LCL5_dp) +  '\n'

  273. 		benchmark_LCL5.write(LCL5_output)

  274. 	# LCL6

  275. 	LCL6_pcr_consensus, LCL6_pcr_free_consensus, LCL6_mendelian_num, LCL6_consensus_call, LCL6_consensus_alt_seq, LCL6_alt, LCL6_dp, LCL6_detected_num = consensus_call(lcl6_list,mendelian_list,row.ALT)

  276. 	if LCL6_mendelian_num != '.':

  277. 		LCL6_output = row.CHROM + '\t' + str(row.POS) + '\t' + '.' + '\t' + row.REF + '\t' + LCL6_consensus_alt_seq + '\t' + '.' + '\t' + '.' + '\t' +'VOTED=' + str(LCL6_mendelian_num) + '\t' + 'GT:ALT:DP' + '\t' + LCL6_consensus_call + ':' + str(LCL6_alt) + ':' + str(LCL6_dp) +  '\n'

  278. 		benchmark_LCL6.write(LCL6_output)

  279. 	# LCL7

  280. 	LCL7_pcr_consensus, LCL7_pcr_free_consensus, LCL7_mendelian_num, LCL7_consensus_call, LCL7_consensus_alt_seq, LCL7_alt, LCL7_dp, LCL7_detected_num = consensus_call(lcl7_list,mendelian_list,row.ALT)

  281. 	if LCL7_mendelian_num != '.':

  282. 		LCL7_output = row.CHROM + '\t' + str(row.POS) + '\t' + '.' + '\t' + row.REF + '\t' + LCL7_consensus_alt_seq + '\t' + '.' + '\t' + '.' + '\t' +'VOTED=' + str(LCL7_mendelian_num) + '\t' + 'GT:ALT:DP' + '\t' + LCL7_consensus_call + ':' + str(LCL7_alt) + ':' + str(LCL7_dp) +  '\n'

  283. 		benchmark_LCL7.write(LCL7_output)

  284. 	# LCL8

  285. 	LCL8_pcr_consensus, LCL8_pcr_free_consensus, LCL8_mendelian_num, LCL8_consensus_call, LCL8_consensus_alt_seq, LCL8_alt, LCL8_dp, LCL8_detected_num = consensus_call(lcl8_list,mendelian_list,row.ALT)

  286. 	if LCL8_mendelian_num != '.':

  287. 		LCL8_output = row.CHROM + '\t' + str(row.POS) + '\t' + '.' + '\t' + row.REF + '\t' + LCL8_consensus_alt_seq + '\t' + '.' + '\t' + '.' + '\t' +'VOTED=' + str(LCL8_mendelian_num) + '\t' + 'GT:ALT:DP' + '\t' + LCL8_consensus_call + ':' + str(LCL8_alt) + ':' + str(LCL8_dp) +  '\n'

  288. 		benchmark_LCL8.write(LCL8_output)

  289. 	# all data

  290. 	all_output = row.CHROM + '\t' + str(row.POS) + '\t' + LCL5_pcr_consensus + '\t' + LCL5_pcr_free_consensus + '\t' + str(LCL5_mendelian_num) + '\t' + LCL5_consensus_call + '\t' + LCL5_consensus_alt_seq + '\t' + str(LCL5_alt) + '\t' + str(LCL5_dp)  + '\t' + str(LCL5_detected_num) + '\t' +\

  291. 					LCL6_pcr_consensus + '\t' + LCL6_pcr_free_consensus + '\t' + str(LCL6_mendelian_num) + '\t' + LCL6_consensus_call + '\t' + LCL6_consensus_alt_seq + '\t' + str(LCL6_alt) + '\t' + str(LCL6_dp)  + '\t' + str(LCL6_detected_num) + '\t' +\

  292. 					LCL7_pcr_consensus + '\t' + LCL7_pcr_free_consensus + '\t' + str(LCL7_mendelian_num) + '\t' + LCL7_consensus_call + '\t' + LCL7_consensus_alt_seq + '\t' + str(LCL7_alt) + '\t' + str(LCL7_dp) + '\t' + str(LCL7_detected_num) + '\t' +\

  293. 					LCL8_pcr_consensus + '\t' + LCL8_pcr_free_consensus + '\t' + str(LCL8_mendelian_num) + '\t' + LCL8_consensus_call + '\t' + LCL8_consensus_alt_seq + '\t' + str(LCL8_alt) + '\t' + str(LCL8_dp) + '\t' + str(LCL8_detected_num) + '\n'

  294. 	all_sample_outfile.write(all_output)