quartet-wgs-germline-qc/codescripts/extract_tables.py

import json
import pandas as pd
import sys, argparse, os

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")


parser.add_argument('-quality', '--quality_yield', type=str, help='*.quality_yield.txt',  required=True)
parser.add_argument('-depth', '--wgs_metrics', type=str, help='*deduped_WgsMetricsAlgo.txt',  required=True)
parser.add_argument('-aln', '--aln_metrics', type=str, help='*_deduped_aln_metrics.txt',  required=True)
parser.add_argument('-is', '--is_metrics', type=str, help='*_deduped_is_metrics.txt',  required=True)

parser.add_argument('-fastqc', '--fastqc', type=str, help='multiqc_fastqc.txt',  required=True)
parser.add_argument('-fastqscreen', '--fastqscreen', type=str, help='multiqc_fastq_screen.txt',  required=True)
parser.add_argument('-hap', '--happy', type=str, help='multiqc_happy_data.json',  required=True)

args = parser.parse_args()

# Rename input:
quality_yield_file = args.quality_yield
wgs_metrics_file = args.wgs_metrics
aln_metrics_file = args.aln_metrics
is_metrics_file = args.is_metrics

fastqc_file = args.fastqc
fastqscreen_file = args.fastqscreen
hap_file = args.happy

#############################################
# fastqc
fastqc = pd.read_table(fastqc_file)

#fastqc = dat.loc[:, dat.columns.str.startswith('FastQC')]
#fastqc.insert(loc=0, column='Sample', value=dat['Sample'])
#fastqc_stat = fastqc.dropna()

# qulimap
#qualimap = dat.loc[:, dat.columns.str.startswith('QualiMap')]
#qualimap.insert(loc=0, column='Sample', value=dat['Sample'])
#qualimap_stat = qualimap.dropna()

# fastqc
#dat = pd.read_table(fastqc_file)

#fastqc_module = dat.loc[:, "per_base_sequence_quality":"kmer_content"]
#fastqc_module.insert(loc=0, column='Sample', value=dat['Sample'])
#fastqc_all = pd.merge(fastqc_stat,fastqc_module,  how='outer', left_on=['Sample'], right_on = ['Sample'])

# fastqscreen
dat = pd.read_table(fastqscreen_file)
fastqscreen = dat.loc[:, dat.columns.str.endswith('percentage')]
dat['Sample'] = [i.replace('_screen','') for i in dat['Sample']]
fastqscreen.insert(loc=0, column='Sample', value=dat['Sample'])

# pre-alignment
pre_alignment_dat = pd.merge(fastqc,fastqscreen,how="outer",left_on=['Sample'],right_on=['Sample'])
pre_alignment_dat['FastQC_mqc-generalstats-fastqc-total_sequences'] = pre_alignment_dat['FastQC_mqc-generalstats-fastqc-total_sequences']/1000000
del pre_alignment_dat['FastQC_mqc-generalstats-fastqc-percent_fails']
del pre_alignment_dat['FastQC_mqc-generalstats-fastqc-avg_sequence_length']
del pre_alignment_dat['ERCC percentage']
del pre_alignment_dat['Phix percentage']
del pre_alignment_dat['Mouse percentage']
pre_alignment_dat = pre_alignment_dat.round(2)
pre_alignment_dat.columns = ['Sample','%Dup','%GC','Total Sequences (million)','%Human','%EColi','%Adapter','%Vector','%rRNA','%Virus','%Yeast','%Mitoch','%No hits']
pre_alignment_dat.to_csv('pre_alignment.txt',sep="\t",index=0)



############################
dat = pd.read_table(aln_metrics_file)
dat['PCT_ALIGNED_READS'] = dat["PF_READS_ALIGNED"]/dat["TOTAL_READS"]
aln_metrics = dat[["Sample", "PCT_ALIGNED_READS","PF_MISMATCH_RATE"]]
aln_metrics = aln_metrics * 100

dat = pd.read_table(is_metrics_file)
is_metrics = dat[['Sample', 'MEDIAN_INSERT_SIZE']]

dat = pd.read_table(quality_yield_file)
dat['%Q20'] = dat['Q20_BASES']/dat['TOTAL_BASES']
dat['%Q30'] = dat['Q30_BASES']/dat['TOTAL_BASES']
quality_yield = dat[['Sample','%Q20','%Q30']]
quality_yield = quality_yield * 100

dat = pd.read_table(wgs_metrics_file)
wgs_metrics = dat[['Sample','MEDIAN_COVERAGE','PCT_1X', 'PCT_5X', 'PCT_10X','PCT_30X']]
wgs_metrics['PCT_1X'] = wgs_metrics['PCT_1X'] * 100
wgs_metrics['PCT_5X'] = wgs_metrics['PCT_5X'] * 100
wgs_metrics['PCT_10X'] = wgs_metrics['PCT_10X'] * 100
wgs_metrics['PCT_30X'] = wgs_metrics['PCT_30X'] * 100

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

# benchmark
with open(hap_file) as hap_json:
	happy = json.load(hap_json)
dat =pd.DataFrame.from_records(happy)
dat = dat.loc[:, dat.columns.str.endswith('ALL')]
dat_transposed = dat.T
benchmark = dat_transposed.loc[:,['sample_id','METRIC.Precision','METRIC.Recall']]
benchmark.columns = ['Sample','Precision','Recall']

#output
fastqc_all.to_csv('fastqc.final.result.txt',sep="\t",index=0)
fastqscreen.to_csv('fastqscreen.final.result.txt',sep="\t",index=0)
qualimap_stat.to_csv('qualimap.final.result.txt',sep="\t",index=0)
benchmark.to_csv('benchmark.final.result.txt',sep="\t",index=0)