Merge pull request #2 from BCCDC-PHL/add-bracken

Add bracken

Author: dfornika

README.md

@@ -4,6 +4,14 @@ A generic pipeline that can be run routinely on all Illumina sequence runs, rega
 * Sequence quality information
 * Possible contamination
 
+## Analyses
+
+* Parse run-level QC statistics from the 'InterOp' directory and write to `.csv` and `.json` format.
+* [FastQC](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/): sample-level sequence quality metrics
+* [Kraken2](https://github.com/DerrickWood/kraken2) + [Bracken](https://github.com/jenniferlu717/Bracken): Taxonomic classification
+of reads. Estimation of relative abundances of taxonomic groups (genus, species) in each sample.
+* [MultiQC](https://github.com/ewels/MultiQC): Collect several QC metrics into a single interactive HTML report.
+
 ## Usage
 
 ```

assets/multiqc_config_base.yaml

@@ -0,0 +1,12 @@
+show_analysis_paths: False
+fn_clean_exts:
+  - "_"
+table_columns_visible:
+  Kraken: False
+top_modules:
+  - "interop"
+  - "fastqc"
+remove_sections:
+  - fastqc_per_base_sequence_content
+  - fastqc_per_base_n_content
+  - fastqc_sequence_length_distribution

bin/bracken_top_n_linelist.py

@@ -0,0 +1,52 @@
+#!/usr/bin/env python
+
+import argparse
+import csv
+import sys
+import re
+import json
+
+
+def parse_bracken_report(bracken_report_path):
+    bracken_report_lines = []
+    with open(bracken_report_path, 'r') as f:
+        reader = csv.DictReader(f, dialect='excel-tab')
+        for row in reader:
+            bracken_report_lines.append(row)
+
+    return bracken_report_lines
+        
+
+def main(args):
+    bracken_report = parse_bracken_report(args.bracken_report)
+
+    bracken_report_sorted = sorted(bracken_report, key=lambda k: k['fraction_total_reads'], reverse=True) 
+    
+    output_fields = ['sample_id', 'taxonomy_level']
+    output_line = {
+        'sample_id': args.sample_id,
+        'taxonomy_level': bracken_report_sorted[0]['taxonomy_lvl']
+    }
+    
+    for n in range(args.top_n):
+        num = str(n + 1)
+        name_field = 'abundance_' + num + '_name'
+        output_line[name_field] = bracken_report_sorted[n]['name']
+        output_fields.append(name_field)
+        fraction_total_reads_field = 'abundance_' + num + '_fraction_total_reads'
+        output_line[fraction_total_reads_field] = bracken_report_sorted[n]['fraction_total_reads']
+        output_fields.append(fraction_total_reads_field)
+        
+
+    csv.register_dialect('unix-csv-quote-minimal', delimiter=',', doublequote=False, lineterminator='\n', quoting=csv.QUOTE_MINIMAL)
+    writer = csv.DictWriter(sys.stdout, fieldnames=output_fields, dialect='unix-csv-quote-minimal')
+    writer.writeheader()
+    writer.writerow(output_line)
+        
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser()
+    parser.add_argument('bracken_report')
+    parser.add_argument('-s', '--sample-id')
+    parser.add_argument('-n', '--top-n', type=int)
+    args = parser.parse_args()
+    main(args)

bin/get_read_length.py

@@ -0,0 +1,33 @@
+#!/usr/bin/env python
+
+import argparse
+import json
+
+
+def main(args):
+  with open(args.sample_sheet_json, 'r') as f:
+    sample_sheet = json.load(f)
+  actual_read_length = sample_sheet['reads'][0]
+
+  if actual_read_length < 60:
+      read_length = 50
+  elif actual_read_length < 90:
+      read_length = 75
+  elif actual_read_length < 125:
+    read_length = 100
+  elif actual_read_length < 175:
+    read_length = 150
+  elif actual_read_length < 225:
+    read_length = 200
+  elif actual_read_length < 275:
+    read_length = 250
+  else:
+    read_length = 300
+  print(read_length)
+
+
+if __name__ == '__main__':
+  parser = argparse.ArgumentParser()
+  parser.add_argument('sample_sheet_json')
+  args = parser.parse_args()
+  main(args)

bin/kraken_parser.py

@@ -0,0 +1,64 @@
+#!/usr/bin/env python3
+
+import argparse
+import json
+import sys
+
+
+def main(args):
+
+    taxonomic_levels = {'U', 'D', 'K', 'P', 'C', 'O', 'F', 'G', 'S'}
+    # assert args.taxonomic_level in rank_codes, "Rank must be one of [(U)nclassified, (D)omain, (K)ingdom, (P)hylum, (C)lass, (O)rder, (F)amily, (G)enus, or (S)pecies.]"
+
+    total_reads = 0
+    unclassified_reads = 0
+    total_reads_reported = 0
+    other_reads = 0
+
+    headers = [
+        "percent_reads_in_clade",
+        "num_reads_in_clade",
+        "ncbi_taxonomy_id",
+        "taxon_name",
+    ]
+
+    print('\t'.join(headers))
+
+    with open(args.kraken_report, 'r') as f:
+        for line in f:
+            record = {}
+            split_line = [record.strip() for record in line.strip().split("\t")]
+
+            record['percent_reads_clade'] = float(split_line[0])
+            record['num_reads_clade'] = int(split_line[1])
+            record['num_reads_taxon'] = int(split_line[2])
+            record['taxonomic_level'] = split_line[3]
+            record['ncbi_taxonomy_id'] = split_line[4]
+            record['clade_name'] = split_line[5]
+
+            if record['clade_name'] == "unclassified":
+                unclassified_reads = record['num_reads_clade']
+                total_reads += record['num_reads_clade']
+            elif record['clade_name'] == "root":
+                total_reads += record['num_reads_clade']
+
+            # record['percent_reads_in_clade'] = record['num_reads_clade'] / float(total_reads) * 100
+
+            if record['taxonomic_level'] == args.taxonomic_level and record['percent_reads_clade'] >= args.threshold_percent:
+                print('\t'.join([str(record['percent_reads_clade']), str(record['num_reads_clade']), str(record['ncbi_taxonomy_id']), record['clade_name']]))
+                total_reads_reported += record['num_reads_clade']
+
+
+    print(str('%.3f' % (unclassified_reads / float(total_reads) * 100)) + '\t' + str(unclassified_reads) + '\t\t' + 'unclassified')
+    total_reads_reported += unclassified_reads
+    other_reads = total_reads - total_reads_reported
+    print(str('%.3f' % (other_reads / float(total_reads) * 100)) + '\t' + str(other_reads) + '\t\t' + 'other')
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser()
+    parser.add_argument('kraken_report')
+    parser.add_argument('-l', '--taxonomic_level')
+    parser.add_argument('-t', '--threshold_percent', type=float)
+    args = parser.parse_args()
+    main(args)

bin/parse_run_summary.py

@@ -8,24 +8,28 @@
 def parse_read_summary(summary_path):
     read_summary_headers = []
     read_summary_lines = []
-    # Basic approach to parsing text between two specific lines
-    # described here: https://stackoverflow.com/a/7559542/780188
+
+    replaced_fields = {'%>=q30': 'percent_greater_than_q30',
+                       '%_occupied': 'percent_occupied'}
+
     with open(summary_path) as summary:
         for line in summary:
             if re.match("^Level", line):
                 read_summary_headers = re.split("\s*,", line.rstrip())
                 read_summary_headers = [
                     x.lower().replace(" ", "_") for x in read_summary_headers
                 ]
-                read_summary_headers = [
-                    x.replace("%>=q30", "percent_greater_than_q30") for x in read_summary_headers
-                ]
+                for idx, header in enumerate(read_summary_headers):
+                    if header in replaced_fields:
+                        read_summary_headers[idx] = replaced_fields[header]
+                
                 break
         for line in summary:
             if re.match("^Total", line):
-                read_summary_lines.append(re.split("\s*,", line.rstrip()))
+                read_summary_lines.append(re.split(",", line.rstrip()))
                 break
-            read_summary_lines.append(re.split("\s*,", line.rstrip()))
+            else:
+                read_summary_lines.append(re.split(",", line.rstrip()))
 
     read_summary = []
     for line in read_summary_lines:
@@ -41,14 +45,15 @@ def parse_read_summary(summary_path):
 
     return read_summary
 
+
 def parse_read_summary_detail(summary_path):
     headers = [
         'lane',
         'surface',
         'tiles',
         'density',
         'clusters_passing_filter',
-        'legacy_pasing_prephasing_rate',
+        'legacy_phasing_prephasing_rate',
         'phasing_slope_offset',
         'prephasing_slope_offset',
         'reads',
@@ -64,6 +69,31 @@ def parse_read_summary_detail(summary_path):
         'percent_occupied',
         'intensity_at_cycle_1',
     ]
+    average_stdev_fields = [
+        'aligned',
+        'clusters_passing_filter',
+        'density',
+        'error',
+        'error_100',
+        'error_75',
+        'error_35',
+        'intensity_at_cycle_1',
+        'percent_occupied',
+    ]
+    slash_fields = { 'legacy_phasing_prephasing_rate': {'numerator_field': 'legacy_phasing_rate',
+                                                        'denominator_field': 'legacy_prephasing_rate'},
+                     'phasing_slope_offset': {'numerator_field': 'phasing_slope',
+                                              'denominator_field': 'phasing_offset'},
+                     'prephasing_slope_offset': {'numerator_field': 'prephasing_slope',
+                                                 'denominator_field': 'prephasing_offset'},
+    }
+    float_fields = [
+        'percent_greater_than_q30',
+        'reads',
+        'reads_passing_filter',
+        'yield',
+    ]
+    
     lines_by_read = {
         'read_1': [],
         'read_i1': [],
@@ -73,15 +103,13 @@ def parse_read_summary_detail(summary_path):
     with open(summary_path) as summary:
         current_read = None
         for line in summary:
-            if re.match("^Read 1$", line):
+            if re.match("^Read 1\n$", line):
                 current_read = 'read_1'
-            elif re.match("^Read 2 \(I\)$", line):
+            elif re.match("^Read 2 \(I\)\n$", line):
                 current_read = 'read_i1'
-            elif re.match("^Read 3 \(I\)$", line):
+            elif re.match("^Read 3 \(I\)\n$", line):
                 current_read = 'read_i2'
-            elif re.match("^Read 4$", line):
-                current_read = 'read_2'
-            elif re.match("^Read 4$", line):
+            elif re.match("^Read 4$\n", line):
                 current_read = 'read_2'
             elif re.match("^Extracted", line) or re.match("^Called", line) or re.match("^Scored", line):
                 current_read = None
@@ -91,15 +119,33 @@ def parse_read_summary_detail(summary_path):
                 lines_by_read[current_read].append(read_line_dict)
             else:
                 pass
+
+        for field in average_stdev_fields:
+            string_value = read_line_dict[field]
+            [average, stdev] = [float(value) for value in string_value.split(' +/- ')]
+            read_line_dict[field] = { 'average': average,
+                                      'stdev': stdev }
+
+        for field, num_denom in slash_fields.items():
+            string_value = read_line_dict[field]
+            numerator_field = num_denom['numerator_field']
+            denominator_field = num_denom['denominator_field']
+            [numerator, denominator] = [float(value) for value in string_value.split(' / ')]
+            read_line_dict[numerator_field] = numerator
+            read_line_dict[denominator_field] = denominator
+            read_line_dict.pop(field, None)
 
     return lines_by_read
 
+
 def main(args):
     read_summary = parse_read_summary(args.summary)
     read_summary_detail = parse_read_summary_detail(args.summary)
 
+    output = {'read_summary': read_summary,
+              'read_details': read_summary_detail}
     # print(json.dumps(read_summary))
-    print(json.dumps(read_summary_detail))
+    print(json.dumps(output))
 
 if __name__ == '__main__':
     parser = argparse.ArgumentParser()

main.nf

@@ -7,14 +7,19 @@ include { multiqc } from './modules/multiqc.nf'
 include { interop_summary } from './modules/interop.nf'
 include { parse_sample_sheet } from './modules/sample-sheet.nf'
 include { kraken2 } from './modules/kraken2.nf'
+include { bracken } from './modules/bracken.nf'
+include { abundance_top_n } from './modules/bracken.nf'
 
 workflow {
   ch_fastq = Channel.fromFilePairs( "${params.run_dir}/Data/Intensities/BaseCalls/*_{R1,R2}_*.fastq.gz" )
   ch_sample_sheet = Channel.fromPath( "${params.run_dir}/SampleSheet.csv" )
+  ch_multiqc_config = Channel.fromPath( "${projectDir}/assets/multiqc_config_base.yaml" )
   ch_run_dir = Channel.fromPath(params.run_dir)
   ch_run_id = Channel.fromPath(params.run_dir).map{ it -> it.baseName }
   ch_kraken2_db = Channel.fromPath(params.kraken2_db)
-  
+  ch_bracken_db = Channel.fromPath(params.bracken_db)
+  ch_taxonomic_levels = Channel.from('Genus', 'Species')
+
   main:
     interop_summary(ch_run_id.combine(ch_run_dir))
 
@@ -23,15 +28,17 @@ workflow {
     fastqc(ch_fastq)
 
     kraken2(ch_fastq.combine(ch_kraken2_db))
+    bracken(kraken2.out.combine(ch_bracken_db).combine(parse_sample_sheet.out).combine(ch_taxonomic_levels))
+
+    abundance_top_n(bracken.out)
 
-    // Line below is just composing the run_id and the list of fastqc_outdirs into a new list. There must be a better way(?)
-    // "run_id" + ["fastqc_outdir1", "fastqc_outdir2", ...] => ["run_id", ["fastqc_outdir1", "fastqc_outdir2"]]
+    abundance_top_n.out.filter{ it[2] == 'Genus' }.map{ it -> it[1] }.collectFile(keepHeader: true, sort: { it.text }, name: "top_3_abundances_genus.csv", storeDir: "${params.outdir}/abundance_top_n")
+    abundance_top_n.out.filter{ it[2] == 'Species' }.map{ it -> it[1] }.collectFile(keepHeader: true, sort: { it.text }, name: "top_5_abundances_species.csv", storeDir: "${params.outdir}/abundance_top_n")
 
     ch_fastqc_collected = fastqc.out.map{ it -> [it[1], it[2]] }.collect()
-    ch_kraken2_collected = kraken2.out.collect()
-    ch_all_qc_outputs = ch_fastqc_collected.combine(ch_kraken2_collected)
+    ch_bracken_species_multiqc_collected = bracken.out.filter{ it[4] == 'Species' }.map{ it -> it[2] }.collect()
+    ch_all_qc_outputs = interop_summary.out.map{ it -> it.drop(1) }.combine(ch_fastqc_collected).combine(ch_bracken_species_multiqc_collected)
 
     ch_all_qc_outputs_with_run_id = ch_run_id.combine(ch_all_qc_outputs).map{ it -> [it[0], it.drop(1)] }
-    multiqc(ch_all_qc_outputs_with_run_id)
-
+    multiqc(ch_multiqc_config.combine(ch_all_qc_outputs_with_run_id))
 }