Add GC content to basic QC stats table (#39)

Author: dfornika

bin/summarize_seqtk_fqchk.py

@@ -9,14 +9,17 @@ def parse_seqtk_fqchk_output(seqtk_fqchk_output_path, quality_threshold):
     with open(seqtk_fqchk_output_path, 'r') as f:
         reader = csv.DictReader(f)
         for row in reader:
-            num_bases_and_avg_q = {}
+            parsed_row = {}
             if row['position'] == 'ALL':
                 percent_above_header = 'percent_bases_above_q' + str(quality_threshold)
-                num_bases_and_avg_q['num_bases'] = int(row['num_bases'])
-                num_bases_and_avg_q['average_q'] = float(row['average_q'])
-                num_bases_and_avg_q[percent_above_header] = float(row[percent_above_header])
-                output.append(num_bases_and_avg_q)
-                
+                parsed_row['num_bases'] = int(row['num_bases'])
+                parsed_row['average_q'] = float(row['average_q'])
+                parsed_row[percent_above_header] = float(row[percent_above_header])
+                percent_g = float(row['percent_g'])
+                percent_c = float(row['percent_c'])
+                parsed_row['percent_gc'] = percent_g + percent_c
+                output.append(parsed_row)
+
     return output      
 
 
@@ -34,6 +37,8 @@ def main(args):
 
     total_bases = sum([x['num_bases'] for x in seqtk_fqchk_output])
 
+    overall_percent_gc = sum([x['percent_gc'] * x['num_bases'] for x in seqtk_fqchk_output]) / total_bases
+
     overall_average_q = sum([x['average_q'] * x['num_bases'] for x in seqtk_fqchk_output]) / total_bases
 
     percent_above_header = 'percent_bases_above_q' + str(quality_threshold)
@@ -42,13 +47,15 @@ def main(args):
 
     print(','.join([
         'sample_id',
+        'percent_gc',
         'total_bases',
         'average_base_quality',
         percent_above_header,
     ]))
 
     print(','.join([
         args.sample_id,
+        str(round(overall_percent_gc, 3)),
         str(total_bases),
         str(round(overall_average_q, 3)),
         str(round(overall_percent_above_threshold, 3)),

modules/combine_qc_stats.nf

@@ -15,7 +15,7 @@ process combine_qc_stats {
   printf "sample_id\\n${sample_id}\\n" > sample_id.csv
   awk -F ',' 'BEGIN {OFS=FS}; NR==1 { for (i=1; i<=NF; i++) {idx[\$i] = i} }; { print \$(idx["abundance_1_name"]), \$(idx["abundance_1_fraction_total_reads"]) }' ${species_abundance} | sed -s 's/abundance_1/most_abundant_species/g' > species_abundance_stats.csv
   awk -F ',' 'BEGIN {OFS=FS}; NR==1 { for (i=1; i<=NF; i++) {idx[\$i] = i} }; { print \$(idx["estimated_genome_size_bp"]), \$(idx["estimated_depth_coverage"]) }' ${estimated_coverage} > estimated_coverage_stats.csv
-  awk -F ',' 'BEGIN {OFS=FS}; NR==1 { for (i=1; i<=NF; i++) {idx[\$i] = i} }; { print \$(idx["total_bases"]), \$(idx["average_base_quality"]), \$(idx["percent_bases_above_q${params.seqtk_fqchk_threshold}"]) }' ${sequence_quality} > sequence_quality_stats.csv
+  awk -F ',' 'BEGIN {OFS=FS}; NR==1 { for (i=1; i<=NF; i++) {idx[\$i] = i} }; { print \$(idx["total_bases"]), \$(idx["average_base_quality"]), \$(idx["percent_bases_above_q${params.seqtk_fqchk_threshold}"]), \$(idx["percent_gc"]) }' ${sequence_quality} > sequence_quality_stats.csv
   paste -d ',' sample_id.csv species_abundance_stats.csv estimated_coverage_stats.csv sequence_quality_stats.csv > ${sample_id}_combined_qc_stats.csv
   """
 }