renluyao
/
quartet-wes-germline-data-generation-qc
Obserwuj 1
Polub 0
Forkuj 0
Kod Zgłoszenia 0 Oczekujące zmiany 0 Wydania 0 Wiki Aktywność
Nie możesz wybrać więcej, niż 25 tematów Tematy muszą się zaczynać od litery lub cyfry, mogą zawierać myślniki ('-') i mogą mieć do 35 znaków.
1 Commit
1 Gałąź
Gałąź: master
Gałęzie Tagi
${ item.name }
Utwórz gałąź ${ searchTerm }
z 'master'
${ noResults }
quartet-wes-germline-data-g.../codescripts/extract_tables.py

157 lines
7.8KB
Czysty Bezpośredni odnośnik Wina Historia 

  1. import json

  2. import pandas as pd

  3. from functools import reduce

  4. import sys, argparse, os

  5. 

  6. parser = argparse.ArgumentParser(description="This script is to get information from multiqc and sentieon, output the raw fastq, bam and variants calling (precision and recall) quality metrics")

  7. 

  8. 

  9. parser.add_argument('-quality', '--quality_yield', type=str, help='*.quality_yield.txt',  required=True)

  10. parser.add_argument('-depth', '--wgs_metrics', type=str, help='*deduped_WgsMetricsAlgo.txt',  required=True)

  11. parser.add_argument('-aln', '--aln_metrics', type=str, help='*_deduped_aln_metrics.txt',  required=True)

  12. parser.add_argument('-is', '--is_metrics', type=str, help='*_deduped_is_metrics.txt',  required=True)

  13. 

  14. parser.add_argument('-fastqc', '--fastqc', type=str, help='multiqc_fastqc.txt',  required=True)

  15. parser.add_argument('-fastqscreen', '--fastqscreen', type=str, help='multiqc_fastq_screen.txt',  required=True)

  16. parser.add_argument('-hap', '--happy', type=str, help='multiqc_happy_data.json',  required=True)

  17. 

  18. parser.add_argument('-project', '--project_name', type=str, help='project_name',  required=True)

  19. 

  20. args = parser.parse_args()

  21. 

  22. # Rename input:

  23. quality_yield_file = args.quality_yield

  24. wgs_metrics_file = args.wgs_metrics

  25. aln_metrics_file = args.aln_metrics

  26. is_metrics_file = args.is_metrics

  27. 

  28. fastqc_file = args.fastqc

  29. fastqscreen_file = args.fastqscreen

  30. hap_file = args.happy

  31. 

  32. project_name = args.project_name

  33. 

  34. #############################################

  35. # fastqc

  36. fastqc = pd.read_table(fastqc_file)

  37. 

  38. #fastqc = dat.loc[:, dat.columns.str.startswith('FastQC')]

  39. #fastqc.insert(loc=0, column='Sample', value=dat['Sample'])

  40. #fastqc_stat = fastqc.dropna()

  41. 

  42. # qulimap

  43. #qualimap = dat.loc[:, dat.columns.str.startswith('QualiMap')]

  44. #qualimap.insert(loc=0, column='Sample', value=dat['Sample'])

  45. #qualimap_stat = qualimap.dropna()

  46. 

  47. # fastqc

  48. #dat = pd.read_table(fastqc_file)

  49. 

  50. #fastqc_module = dat.loc[:, "per_base_sequence_quality":"kmer_content"]

  51. #fastqc_module.insert(loc=0, column='Sample', value=dat['Sample'])

  52. #fastqc_all = pd.merge(fastqc_stat,fastqc_module,  how='outer', left_on=['Sample'], right_on = ['Sample'])

  53. 

  54. # fastqscreen

  55. dat = pd.read_table(fastqscreen_file)

  56. fastqscreen = dat.loc[:, dat.columns.str.endswith('percentage')]

  57. dat['Sample'] = [i.replace('_screen','') for i in dat['Sample']]

  58. fastqscreen.insert(loc=0, column='Sample', value=dat['Sample'])

  59. 

  60. # pre-alignment

  61. pre_alignment_dat = pd.merge(fastqc,fastqscreen,how="outer",left_on=['Sample'],right_on=['Sample'])

  62. pre_alignment_dat['FastQC_mqc-generalstats-fastqc-total_sequences'] = pre_alignment_dat['FastQC_mqc-generalstats-fastqc-total_sequences']/1000000

  63. del pre_alignment_dat['FastQC_mqc-generalstats-fastqc-percent_fails']

  64. del pre_alignment_dat['FastQC_mqc-generalstats-fastqc-avg_sequence_length']

  65. del pre_alignment_dat['ERCC percentage']

  66. del pre_alignment_dat['Phix percentage']

  67. del pre_alignment_dat['Mouse percentage']

  68. pre_alignment_dat = pre_alignment_dat.round(2)

  69. pre_alignment_dat.columns = ['Sample','%Dup','%GC','Total Sequences (million)','%Human','%EColi','%Adapter','%Vector','%rRNA','%Virus','%Yeast','%Mitoch','%No hits']

  70. pre_alignment_dat.to_csv('pre_alignment.txt',sep="\t",index=0)

  71. 

  72. 

  73. 

  74. ############################

  75. dat = pd.read_table(aln_metrics_file,index_col=False)

  76. dat['PCT_ALIGNED_READS'] = dat["PF_READS_ALIGNED"]/dat["TOTAL_READS"]

  77. aln_metrics = dat[["Sample", "PCT_ALIGNED_READS","PF_MISMATCH_RATE"]]

  78. aln_metrics = aln_metrics * 100

  79. aln_metrics['Sample'] = [x[-1] for x in aln_metrics['Sample'].str.split('/')]

  80. 

  81. dat = pd.read_table(is_metrics_file,index_col=False)

  82. is_metrics = dat[['Sample', 'MEDIAN_INSERT_SIZE']]

  83. is_metrics['Sample'] = [x[-1] for x in is_metrics['Sample'].str.split('/')]

  84. 

  85. dat = pd.read_table(quality_yield_file,index_col=False)

  86. dat['%Q20'] = dat['Q20_BASES']/dat['TOTAL_BASES']

  87. dat['%Q30'] = dat['Q30_BASES']/dat['TOTAL_BASES']

  88. quality_yield = dat[['Sample','%Q20','%Q30']]

  89. quality_yield = quality_yield * 100

  90. quality_yield['Sample'] = [x[-1] for x in quality_yield['Sample'].str.split('/')]

  91. 

  92. dat = pd.read_table(wgs_metrics_file,index_col=False)

  93. wgs_metrics = dat[['Sample','MEAN_COVERAGE','MEDIAN_COVERAGE','PCT_1X', 'PCT_5X', 'PCT_10X','PCT_30X']]

  94. wgs_metrics['PCT_1X'] = wgs_metrics['PCT_1X'] * 100

  95. wgs_metrics['PCT_5X'] = wgs_metrics['PCT_5X'] * 100

  96. wgs_metrics['PCT_10X'] = wgs_metrics['PCT_10X'] * 100

  97. wgs_metrics['PCT_30X'] = wgs_metrics['PCT_30X'] * 100

  98. wgs_metrics['Sample'] = [x[-1] for x in wgs_metrics['Sample'].str.split('/')]

  99. 

  100. data_frames = [aln_metrics, is_metrics, quality_yield, wgs_metrics]

  101. post_alignment_dat = reduce(lambda  left,right: pd.merge(left,right,on=['Sample'],how='outer'), data_frames)

  102. post_alignment_dat.columns = ['Sample', '%Mapping', '%Mismatch Rate', 'Mendelian Insert Size','%Q20', '%Q30', 'Mean Coverage' ,'Median Coverage', 'PCT_1X', 'PCT_5X', 'PCT_10X','PCT_30X']

  103. post_alignment_dat = post_alignment_dat.round(2)

  104. post_alignment_dat.to_csv('post_alignment.txt',sep="\t",index=0)

  105. 

  106. 

  107. 

  108. #########################################

  109. # variants calling

  110. with open(hap_file) as hap_json:

  111. 	happy = json.load(hap_json)

  112. dat =pd.DataFrame.from_records(happy)

  113. dat = dat.loc[:, dat.columns.str.endswith('ALL')]

  114. dat_transposed = dat.T

  115. dat_transposed = dat_transposed.loc[:,['sample_id','TRUTH.FN','QUERY.TOTAL','QUERY.FP','QUERY.UNK','METRIC.Precision','METRIC.Recall','METRIC.F1_Score']]

  116. dat_transposed['QUERY.TP'] = dat_transposed['QUERY.TOTAL'].astype(int) - dat_transposed['QUERY.UNK'].astype(int) - dat_transposed['QUERY.FP'].astype(int)

  117. dat_transposed['QUERY'] =dat_transposed['QUERY.TOTAL'].astype(int) - dat_transposed['QUERY.UNK'].astype(int)

  118. indel = dat_transposed[['INDEL' in s for s in dat_transposed.index]]

  119. snv = dat_transposed[['SNP' in s for s in dat_transposed.index]]

  120. indel.reset_index(drop=True, inplace=True)

  121. snv.reset_index(drop=True, inplace=True)

  122. benchmark = pd.concat([snv, indel], axis=1)

  123. benchmark = benchmark[["sample_id", 'QUERY.TOTAL', 'QUERY','QUERY.TP','QUERY.FP','TRUTH.FN','METRIC.Precision', 'METRIC.Recall','METRIC.F1_Score']]

  124. benchmark.columns = ['Sample','sample_id','SNV number','INDEL number','SNV query','INDEL query','SNV TP','INDEL TP','SNV FP','INDEL FP','SNV FN','INDEL FN','SNV precision','INDEL precision','SNV recall','INDEL recall','SNV F1','INDEL F1']

  125. benchmark = benchmark[['Sample','SNV number','INDEL number','SNV query','INDEL query','SNV TP','INDEL TP','SNV FP','INDEL FP','SNV FN','INDEL FN','SNV precision','INDEL precision','SNV recall','INDEL recall','SNV F1','INDEL F1']]

  126. benchmark['SNV precision'] = benchmark['SNV precision'].astype(float)

  127. benchmark['INDEL precision'] = benchmark['INDEL precision'].astype(float)

  128. benchmark['SNV recall'] = benchmark['SNV recall'].astype(float)

  129. benchmark['INDEL recall'] = benchmark['INDEL recall'].astype(float)

  130. benchmark['SNV F1'] = benchmark['SNV F1'].astype(float)

  131. benchmark['INDEL F1'] = benchmark['INDEL F1'].astype(float)

  132. benchmark = benchmark.round(2)

  133. benchmark.to_csv('variants.calling.qc.txt',sep="\t",index=0)

  134. 

  135. #all_rep = [x.split('_')[6] for x in benchmark['Sample']]

  136. #if all_rep.count(all_rep[0]) == 4:

  137. #	rep = list(set(all_rep))

  138. #	columns = ['Family','Average Precision','Average Recall','Precison SD','Recall SD']

  139. #	df_ = pd.DataFrame(columns=columns)

  140. #	for i in rep:

  141. #		string = "_" + i + "_"

  142. #		sub_dat = benchmark[benchmark['Sample'].str.contains('_1_')]

  143. #		mean = list(sub_dat.mean(axis = 0, skipna = True))

  144. #		sd = list(sub_dat.std(axis = 0, skipna = True))

  145. #		family_name = project_name + "." + i + ".SNV"

  146. #		df_ = df_.append({'Family': family_name, 'Average Precision': mean[0], 'Average Recall': mean[2], 'Precison SD': sd[0], 'Recall SD': sd[2] }, ignore_index=True)

  147. #		family_name = project_name + "." + i + ".INDEL"

  148. #		df_ = df_.append({'Family': family_name, 'Average Precision': mean[1], 'Average Recall': mean[3], 'Precison SD': sd[1], 'Recall SD': sd[3] }, ignore_index=True)

  149. #	df_ = df_.round(2)

  150. #	df_.to_csv('precision.recall.txt',sep="\t",index=0)

  151. #else:

  152. #	pass

  153. 

  154. 

  155.