Adding result parsers and updates throughout

Author: marchoeppner

README.md

@@ -1,4 +1,12 @@
-# Pipeline
+# GMO-Check pipeline
+
+This pipeline is being developed to detect GMO content from short-read (amplicon) data. Currently, only the detection of GABA mutations in tomato are supported. As this particular modification is characterized by a specific nucleotide insertion into the SIGAB3 gene, this is what the pipeline will currently try to check for. It does this by two independent approaches - classical variant calling on the one hand and the detection of the insertion in merged and dereplicated amplicon "ZOTUS" against a blast database containing the wild type gene sequence on the other. 
+
+The variant calling workflow will first align quality- and adapter trimmed reads against the tomato reference genome (v3.0). It then masks out bases that start and overlap with the curated primer site locations using samtools ampliconclip. No deduplication will be performed to enable the accurate determination of GMO content in mixed samples. Finally, the read alignment is analyzed with Freebayes to determine the presence of any diagnostically relevant variants. 
+
+For the amplicon-assembly approach, primer sequences are stripped from the reads using Ptrimmer. The stripped reads are then merged, filtered and reduced to unique "ZOTUS" with Vsearch. The resulting sequences are blasted against a  built-in database to check for evidence of the diagnostically relevant insertion. GMO content is determined by extracting the overall raw read count as annotated into the assembled and dereplicated reads for which a positive signal was determined with Blast versus the total number of raw reads represented in the reduced amplicon sequence set. 
+
+Using a set of 126 samples, both approaches yielded very similar estimates for % GMO content in given sample (+/- 1%). 
 
 ## Documentation 
 

assets/genomes/test/amplicon.txt

@@ -0,0 +1,3 @@
+CGAACCCTAGCAGATCGTCT	TCAAAACAACCATTAATCCTTCCCT	162	IL.S.tSIGAD3
+AAGACAATAGCCTCCACAACG	AGTCAGTACAAGACATAATAATACAAAGAG	438	N028_SiGAD3_N-term-seq2
+AGGGATATCGAAATGTAATGGAAAATTG	CAATTCAATAGAACAAAGGATGATACATTC	510	N029_SiGAD3_N-term-seq1

assets/genomes/test/genome.fa

@@ -0,0 +1,23 @@
+>SiGAD3|NM_001246898.2
+ATGGTTCTCTCAAAAACTCCTTCTGATGATTCTGTACACTCCACATTTGCTTCTCGCTATGTTCGAACTT
+CACTACCAAGGTTTGAGATGCTAGAGAAGTCTATACCAAAAGAGGCAGCATACCAAATGATTAATGATGA
+GTTAATGCTTGATGGGAATCCAAGGTTAAATTTGGCATCATTTGTAACCACATGGATGGAACCAGAATGT
+GATAAGCTTATGATGGCTTCAATTAACAAGAATTATGTTGACATGGATGAATACCCTGTCACCACTGAGC
+TTCAGAATCGATGTGTAAACATGATAGCGCGTTTATTCAATGCGCCTTTGAAAGAGGAAGAAATAGGAAT
+TGGTGTGGGGACAGTGGGGTCATCAGAGGCCATAATGTTAGCAGGGCTGGCCTTCAAGAGGAACTGGCAA
+AACAAACGCAAAGCTGAGGGAAAGCCTTATGATAAGCCCAACATTGTCACTGGTGCTAATGTTCAGGTGT
+GTTGGGAGAAATTTGCAAACTACTTTGAAGTGGAATTGAAACAAGTCAAGTTAAGTGAAGGGTACTATGT
+GATGGACCCAATCAAAGCTGTGGAAATGGTAGATGACAACACTATTTGTGTTGCTGCTATTTTGGGTTCA
+ACACTTAATGGAGAATTTGAAGATGTCAAACTCTTGAATGATCTTTTGATTGAAAAGAATAAACAAACTG
+GATGGGACACACCTATTCATGTGGATGCAGCAAGTGGTGGATTCATTGCACCATTTATCTATCCAGAGTT
+GGAATGGGATTTTAGGCTTCCTTTAGTGAAAAGTATTAATGTGAGTGGACACAAATATGGGCTTGTTTAT
+GCTGGTATTGGTTGGGTTATTTGGAGAACTAAACAAGACTTGCCTCAACAACTCATTTTTCATATCAATT
+ATCTTGGTGCTGATCAGCCTACTTTTACTCTCAATTTCTCTAAAGGTTCAAGTCAAGTCATTGCTCAATA
+TTATCAGCTTATCCGCTTGGGCTATGAGGGATATCGAAATGTAATGGAAAATTGTCGTGAAAATGCAATT
+GTGCTAAGAAAAGGACTTGAAAAAACAGGACGTTTCAATATAATCTCCAAAGATGAAGGTATACCCTTGG
+TGGCATTTTCCCTCAAAGACAATAGCCTCCACAACGAATTCGAGGTCTCTGAGACCCTCCGTAGGTTTGG
+GTGGATTGTCCCAGCCTACACTATGCCAGCTGACCTGCAACATGTTACAGTGTTGCGCGTTGTGATTAGA
+GAGGACTTCTCCCGAACCCTAGCAGATCGTCTTGTCTCTGACATCGTCAAGGTCCTCCACGAGCTCCCGA
+ATGCCAAAAAAGTGGAGGATAATTTGATGATCAATAATGAGAAGAAAACAGAAATTGAAGTTCAAAGGGC
+AATTGCTGAGTTTTGGAAGAAATATGTTTTAGCTAGGAAAGCATCTATTTGTTAGGGAAGGATTAATGGT
+TGTTTTGAAGGAATAAGTATTAATTACTAGTAACGTTTTGGTATTAATTATAAAAAATGTG

assets/genomes/test/primers.bed

@@ -0,0 +1,7 @@
+1	14628	14648	N028_SiGAD3_N-term-seq2_LEFT	60	+
+1	15086	15115	N030_SiGAD3_N-term-seq2_RIGHT	60	-
+1	14499	14526	N027_SiGAD3_N-term-seq1_LEFT	60	+
+1	15036	15065	N029_SiGAD3_N-term-seq1_RIGHT	60	-
+1	14765	14784	IL.S.t.SIGAD3.5s	60	+
+1	14946	14970	IL.S.t.SIGAD3.4as	60	-
+

assets/genomes/test/rules.json

@@ -0,0 +1,29 @@
+{
+    "rules": {
+        "vsearch-blast": {
+            "payload": [
+                {
+                    "format": "JSON",
+                    "name": "GABA Mutation in SIGAD3",
+                    "target": "SiGAD3|NM_001246898.2",
+                    "matcher": "AAAG-TGGA",
+                    "positive_report": "Diese Probe enthält eine GABA Mutation in SIGAD3. Nachweis erbraucht über: Amplicon Analyse.",
+                    "negative_report": "Für diese Probe konnte keine GABA Mutation in SIGAD3 nachgewiesen werden."
+                }
+            ]
+            
+        },
+        "bwa-freebayes": {
+            "payload": [
+                {
+                    "format": "VCF",
+                    "target": "1:14834",
+                    "name": "GABA Mutation in SIGAD3",
+                    "matcher": "1\t14834\t.\tGTG\tGTTG",
+                    "positive_report": "Diese Probe enthält eine GABA Mutation in SIGAD3. Nachweis erbracht über: Varianten Analyse.",
+                    "negative_report": "Für diese Probe konnte keine GABA Mutation in SIGAD3 nachgewiesen werden."
+                } 
+            ]
+        }
+    }
+}

assets/genomes/test/targets.bed

@@ -0,0 +1 @@
+1	14784	14946	SIGAD3

assets/genomes/tomato/rules.json

@@ -3,10 +3,12 @@
         "vsearch-blast": {
             "payload": [
                 {
-                    "format": "XML",
+                    "format": "JSON",
+                    "name": "GABA Mutation in SIGAD3",
                     "target": "SiGAD3|NM_001246898.2",
                     "matcher": "AAAG-TGGA",
-                    "yields": "Diese Probe enthält eine GABA Mutation in SIGAD3. Nachweis erbraucht über: Amplicon Analyse."
+                    "positive_report": "Diese Probe enthält eine GABA Mutation in SIGAD3. Nachweis erbraucht über: Amplicon Analyse.",
+                    "negative_report": "Für diese Probe konnte keine GABA Mutation in SIGAD3 nachgewiesen werden."
                 }
             ]
 
@@ -16,8 +18,10 @@
                 {
                     "format": "VCF",
                     "target": "1:14834",
+                    "name": "GABA Mutation in SIGAD3",
                     "matcher": "1\t14834\t.\tGTG\tGTTG",
-                    "yields": "Diese Probe enthält eine GABA Mutation in SIGAD3. Nachweis erbracht über: Varianten Analyse."    
+                    "positive_report": "Diese Probe enthält eine GABA Mutation in SIGAD3. Nachweis erbracht über: Varianten Analyse.",
+                    "negative_report": "Für diese Probe konnte keine GABA Mutation in SIGAD3 nachgewiesen werden."
                 } 
             ]
         }

bin/analyze_blast.rb

@@ -2,7 +2,7 @@
 
 require 'optparse'
 require 'ostruct'
-require 'nokogiri'
+require 'json'
 
 ### Define modules and classes here
 
@@ -11,18 +11,83 @@
 opts = OptionParser.new()
 opts.banner = "Reads Fastq files from a folder and writes a sample sheet to STDOUT"
 opts.separator ""
-opts.on("-b","--blast", "=BLAST","Blast report to read") {|argument| options.vcf = argument }
+opts.on("-b","--blast", "=BLAST","Blast report to read") {|argument| options.blast = argument }
 opts.on("-j","--json", "=JSON","JSON to read") {|argument| options.json = argument }
+opts.on("-s","--sample", "=SAMPLE","Sample name") {|argument| options.sample = argument }
 opts.on("-h","--help","Display the usage information") {
     puts opts
     exit
 }
 
 opts.parse! 
 
+min_coverage = 100
+
+output = { "sample" => options.sample, "matches" => [] }
+
 date = Time.now.strftime("%Y-%m-%d")
 
-file = File.open(options.blast)
+json = JSON.parse(IO.readlines(options.json).join)
+
+rules = json["rules"]["vsearch-blast"]["payload"]
+
+blast = JSON.parse(IO.readlines(options.blast).join)
+
+reports = blast["BlastOutput2"]
+
+findings = []
+total_cov = 0
+carrier_cov = 0
+
+rules.each do |rule|
+
+    total_cov = 0
+    query_cov = 0
+
+    rule_name = rule["name"]
+    rule_string = rule["matcher"]
+
+    has_matched = false
+    reports.each do |r|
+
+        report = r["report"]
+        results = report["results"]["search"]
+
+        query_string = results["query_title"]
+
+        query,coverage = query_string.split(";")
+        coverage = coverage.gsub("size=", "").to_i
+    
+        total_cov += coverage
+    
+        hits = results["hits"]
+    
+        hits.each do |hit|
+    
+            target = hit["description"][0]["title"]
+    
+            hsps = hit["hsps"]
+    
+            hsps.each do |hsp|
+                target_seq = hsp["hseq"]
+
+                if target_seq.include?(rule_string)
+                    has_matched = true
+                    carrier_cov += coverage
+                end
+    
+            end
+        end
+
+    end
+
+    if has_matched
+        perc = (carrier_cov.to_f / total_cov.to_f) * 100
+        output["matches"] << { "rule" => rule_name , "Befund" => rule["positive_report"], "Anteil Variante %" => perc.round(2), "Abdeckung Referenzallel" => total_cov-carrier_cov, "Abdeckung Variantenallel" => carrier_cov }
+    else
+        output["matches"] << { "rule" => rule_name , "Befund" => rule["negative_report"] }
+    end
 
-xml = Nokogiri::XML(file)
+end
 
+puts output.to_json

bin/analyze_vcf.rb

@@ -72,6 +72,7 @@ def parse_vcf(file)
 opts.separator ""
 opts.on("-v","--vcf", "=VCF","VCF to read") {|argument| options.vcf = argument }
 opts.on("-j","--json", "=JSON","JSON to read") {|argument| options.json = argument }
+opts.on("-s","--sample", "=SAMPLE","Sample name") {|argument| options.sample = argument }
 opts.on("-h","--help","Display the usage information") {
     puts opts
     exit
@@ -81,43 +82,58 @@ def parse_vcf(file)
 
 date = Time.now.strftime("%Y-%m-%d")
 
+result = { "sample" => options.sample, "matches" => [] }
+
 json = JSON.parse(IO.readlines(options.json).join)
 
 rules = json["rules"]["bwa-freebayes"]["payload"]
 
 vcf = parse_vcf(options.vcf)
 
-vcf.each do |entry|
-
-    allele = entry.allele_string
+rules.each do  |rule|
+    
+    this_match = {}
 
-    sample_name = entry.sample_names[0]
-    puts ">>>" + sample_name + "<<<"
+    rule_name = rule["name"]
+    rule_string = rule["matcher"]
 
     has_matched = false
 
-    rules.each do |rule|
-        string = rule["matcher"]
-        if string == allele
+    vcf.each do |entry|
+        
+        allele_string = entry.allele_string
+        this_sample = entry.samples[0]
+
+        # A match, presumably
+        if rule_string == allele_string
+
             has_matched = true
-            
-            puts rule["yields"]
-            sample = entry.samples[0]
-            genotype = sample["GT"]
+
+            this_match["rule"] = rule_name
+            this_match["Befund"] = rule["positive_report"]
+
+            genotype = this_sample["GT"]
+
             if genotype == "0/0"
-                puts "Varianten Frequenz unter Detektierungsschwelle!"
+                this_match["Anmerkung"] = "Variantenfrequenz unter Call-Schwelle!"
             end
-            rcov,acov = sample["AD"].split(",")
-            perc = (acov.to_f / rcov.to_f)*100.0
-            puts "\tGenotyp: #{sample["GT"]}\tAnteil: #{perc.round(2)}%\tRef: #{rcov}\tAlt: #{acov}\t"
+
+            rcov,acov = this_sample["AD"].split(",")
+            cov_sum = acov.to_i + rcov.to_i
+            perc = (acov.to_f / cov_sum.to_f)*100.0
+            this_match["Anteil Variante %"] = perc.round(2)
+            this_match["Abdeckung Referenzallel"] = rcov
+            this_match["Abdeckung Variantenallel"] = acov
+
+            result["matches"] << this_match
 
         end
     end
-    
-    if !has_matched
-        puts "Keine GABA Mutation nachgewiesen!"
+
+    unless has_matched
+        result["matches"] << { "rule" => rule_name, "Befund" => rule["negative_report"]}
     end
 
-    puts "==============================================================================="
+end
 
-end
+puts result.to_json

bin/reports_to_table.rb

@@ -0,0 +1,70 @@
+#!/bin/env ruby
+
+require 'optparse'
+require 'ostruct'
+require 'json'
+require 'csv'
+
+### Define modules and classes here
+
+def parse_json(filename)
+
+    return JSON.parse(IO.readlines(filename).join)
+
+end
+
+### Get the script arguments and open relevant files
+options = OpenStruct.new()
+opts = OptionParser.new()
+opts.banner = "Reads reports and makes a table"
+opts.separator ""
+opts.on("-h","--help","Display the usage information") {
+    puts opts
+    exit
+}
+
+opts.parse! 
+
+files = Dir["*.json"]
+
+rows = []
+
+header = [ "Sample", "Blast", "Freebayes" ]
+
+rows << header
+
+files.group_by{|f| f.split(".")[0..-3].join}.each do |group,reports|
+
+    blast = reports.find {|r| r.include?("blast.json")}
+    freebaytes = reports.find { |r| r.include?("freebayes.json")}
+
+    this_data = []
+
+    sample = group
+    this_data << sample
+
+    if blast
+        json = parse_json(blast)
+        matches = json["matches"]
+        this_data << matches[0]["Befund"]
+    else
+        this_data << ""
+    end
+
+    if freebayes
+        json = parse_json(freebayes)
+        matches = json["matches"]
+        this_data << matches[0]["Befund"]
+    else
+        this_data << ""
+    end
+
+    rows << this_data
+
+end
+
+File.write("summary_mqc.csv", rows.map(&:to_csv).join)
+
+
+
+