--- /dev/null
+module VcfPipelineHelper
+ require 'csv'
+
+ def reset_vcf_pipeline_instance(pi, input_manifest)
+ params = {
+ 'PICARD_ZIP' => '7a4073e29bfa87154b7102e75668c454+83+K@van',
+ 'GATK_BUNDLE' => '0a37aaf212464efa2a77ff9ba51c0148+10524+K@van',
+ 'GATK_TAR_BZ2' => '482ebab0408e173370c499f0b7c00878+93+K@van',
+ 'BWA' => '73be5598809c66f260fedd253c8608bd+67+K@van',
+ 'SAM' => '55d2115faa608eb95dab4f875b7511b1+72+K@van',
+ 'REGION_PADDING' => '10',
+ 'REGIONS' => 'e52c086f41c2f089d88ec2bbd45355d3+87+K@van/SeqCap_EZ_Exome_v2.hg19.bed',
+ 'STAND_CALL_CONF' => '4.0',
+ 'STAND_EMIT_CONF' => '4.0',
+ "bwa/INPUT" => input_manifest
+ }
+ pi.components = PipelineTemplate.find(pi.pipeline_uuid).components
+ pi.update_job_parameters(params)
+ pi.active = true
+ pi.success = nil
+ end
+
+ def vcf_pipeline_summary(pi)
+ stats = {}
+ collection_link = Link.
+ where(head_uuid: pi.uuid,
+ link_class: 'client-defined',
+ name: 'vcffarm-pipeline-invocation').
+ last
+ if collection_link
+ stats[:collection_uuid] = collection_link.tail_uuid
+ else
+ pi.components[:steps].each do |step|
+ if step[:name] == 'bwa'
+ step[:params].each do |param|
+ if param[:name] == 'INPUT'
+ stats[:collection_uuid] = param[:data_locator] || param[:value]
+ break
+ end
+ end
+ end
+ end
+ end
+ if stats[:collection_uuid]
+ Link.where(tail_uuid: stats[:collection_uuid],
+ head_kind: Group)[0..0].each do |c2p|
+ stats[:project_uuid] = c2p.head_uuid
+ group = Group.find stats[:project_uuid]
+ stats[:project_name] = group.name rescue nil
+ end
+ Link.where(tail_uuid: stats[:collection_uuid],
+ head_kind: Specimen)[0..0].each do |c2s|
+ stats[:specimen_uuid] = c2s.head_uuid
+ specimen = Specimen.find stats[:specimen_uuid]
+ stats[:specimen_id] = specimen.properties[:specimen_id] rescue nil
+ end
+ end
+ stats[:runtime] = {}
+ stats[:alignment_for_step] = {}
+ stats[:alignment] = {}
+ stats[:coverage] = []
+ pi.components[:steps].each do |step|
+ if step[:warehousejob]
+ if step[:name] == 'bwa' and step[:warehousejob][:starttime]
+ stats[:runtime][:started_at] = step[:warehousejob][:starttime]
+ end
+ if step[:warehousejob][:finishtime]
+ stats[:runtime][:finished_at] =
+ [ step[:warehousejob][:finishtime],
+ stats[:runtime][:finished_at] ].compact.max
+ end
+ end
+ if step[:name] == 'picard-casm' and
+ step[:complete] and
+ step[:output_data_locator]
+ tsv = IO.
+ popen("whget -r #{step[:output_data_locator]}/ -").
+ readlines.
+ collect { |x| x.strip.split "\t" }
+ casm = {}
+ head = []
+ tsv.each do |data|
+ if data.size < 4 or data[0].match /^\#/
+ next
+ elsif data[0] == 'CATEGORY' or data[1].match /[^\d\.]/
+ head = data
+ elsif data[0] == 'PAIR'
+ head.each_with_index do |name, index|
+ x = data[index]
+ if x and x.match /^\d+$/
+ x = x.to_i
+ elsif x and x.match /^\d+\.\d+$/
+ x = x.to_f
+ end
+ name = name.downcase.to_sym
+ casm[name] ||= []
+ casm[name] << x
+ end
+ end
+ end
+ stats[:picard_alignment_summary] = casm
+ end
+ if step[:name] == 'gatk-stats' and
+ step[:complete] and
+ step[:output_data_locator]
+ csv = IO.
+ popen("whget #{step[:output_data_locator]}/mincoverage_nlocus.csv").
+ readlines.
+ collect { |x| x.strip.split ',' }
+ csv.each do |depth, nlocus, percent|
+ stats[:coverage][depth.to_i] = nlocus.to_i
+ end
+ end
+ if step[:name] == 'gatk-realign' and
+ step[:complete] and
+ step[:output_data_locator]
+ logs = IO.
+ popen("whget #{step[:warehousejob][:metakey]}").
+ readlines.
+ collect(&:strip)
+ logs.each do |logline|
+ if (re = logline.match /\s(\d+) stderr INFO .* (\d+) reads were filtered out.*of (\d+) total/)
+ stats[:alignment_for_step][re[1]] ||= {}
+ stats[:alignment_for_step][re[1]][:filtered_reads] = re[2].to_i
+ stats[:alignment_for_step][re[1]][:total_reads] = re[3].to_i
+ elsif (re = logline.match /(\d+) reads.* failing BadMate/)
+ stats[:alignment][:bad_mate_reads] = re[1].to_i
+ elsif (re = logline.match /(\d+) reads.* failing MappingQualityZero/)
+ stats[:alignment][:mapq0_reads] = re[1].to_i
+ end
+ end
+ end
+ if step[:name] == 'gatk-merge-call' and
+ step[:complete] and
+ step[:output_data_locator]
+ stats[:vcf_file_name] = "#{stats[:project_name]}-#{stats[:specimen_id]}-#{step[:output_data_locator][0..31]}.vcf"
+ logs = IO.
+ popen("whget #{step[:warehousejob][:metakey]}").
+ readlines.
+ collect(&:strip)
+ logs.each do |logline|
+ if (re = logline.match /(\d+) reads were filtered out.*of (\d+) total/)
+ stats[:alignment][:filtered_reads] = re[1].to_i
+ stats[:alignment][:total_realigned_reads] = re[2].to_i
+ elsif (re = logline.match /(\d+) reads.* failing BadMate/)
+ stats[:alignment][:bad_mate_reads] = re[1].to_i
+ elsif (re = logline.match /(\d+) reads.* failing UnmappedRead/)
+ stats[:alignment][:unmapped_reads] = re[1].to_i
+ end
+ end
+
+ stats[:chromosome_calls] = {}
+ tsv = IO.
+ popen("whget #{step[:output_data_locator]}/merged.vcf | egrep -v '^#' | cut -f1 | uniq -c").
+ readlines.
+ collect { |x| x.strip.split }
+ tsv.each do |n_variants, sequence_name|
+ stats[:chromosome_calls][sequence_name] = n_variants.to_i
+ end
+
+ stats[:inferred_sex] = false
+ calls = stats[:chromosome_calls]
+ if calls['X'] and calls['X'] > 200
+ if !calls['Y']
+ stats[:inferred_sex] = 'female'
+ elsif calls['Y'] * 60 < calls['X']
+ # if Y < X/60 they are presumed to be misalignments
+ stats[:inferred_sex] = 'female'
+ elsif calls['Y'] * 25 > calls['X']
+ # if Y > X/25 we presume a Y chromosome was present
+ stats[:inferred_sex] = 'male'
+ end
+ end
+ end
+ end
+ stats[:alignment][:total_reads] = 0
+ stats[:alignment][:filtered_reads] ||= 0
+ stats[:alignment][:bad_mate_reads] ||= 0
+ stats[:alignment][:mapq0_reads] ||= 0
+ stats[:alignment_for_step].values.each do |a4s|
+ stats[:alignment][:total_reads] += (a4s[:total_reads] || 0)
+ stats[:alignment][:filtered_reads] += (a4s[:filtered_reads] || 0)
+ stats[:alignment][:bad_mate_reads] += (a4s[:bad_mate_reads] || 0)
+ stats[:alignment][:mapq0_reads] += (a4s[:mapq0_reads] || 0)
+ end
+
+ if stats[:collection_uuid]
+ csv = CSV.parse IO.
+ popen("whget #{stats[:collection_uuid]}/SampleSheet.csv -").
+ read
+ if !csv.empty?
+ pivoted = []
+ csv[0].each_with_index do |head, col|
+ pivoted << csv.collect { |row| row[col] }
+ end
+ stats[:source_data_csv_columns] = pivoted
+ end
+ end
+
+ picardas = stats[:picard_alignment_summary]
+ stats[:summary_csv_columns] =
+ [['PROJECT', stats[:project_name]],
+ ['SPECIMEN', stats[:specimen_id]],
+ ['VCF_FILE_NAME', stats[:vcf_file_name]],
+ ['INFERRED_SEX', stats[:inferred_sex]],
+ ['SOURCE_DATA', stats[:collection_uuid]],
+ ['PIPELINE_UUID', pi.pipeline_uuid],
+ ['PIPELINE_RUN_UUID', pi.uuid],
+ ['PIPELINE_RUN_START', (stats[:runtime][:started_at] rescue nil)],
+ ['PIPELINE_RUN_FINISH', (stats[:runtime][:finished_at] rescue nil)],
+ ['N_READS_RAW',
+ (n_raw = picardas[:total_reads].inject(0,:+) rescue nil)],
+ ['N_READS_MAPPED',
+ (n_mapped = picardas[:reads_aligned_in_pairs].inject(0,:+) rescue nil)],
+ ['PERCENT_READS_MAPPED',
+ (100.0 * n_mapped / n_raw rescue nil)],
+ ['N_READS_ON_TARGET',
+ (n_on_target = stats[:alignment][:total_reads] - stats[:alignment][:filtered_reads] rescue nil)],
+ ['PERCENT_READS_ON_TARGET',
+ (100.0 * n_on_target / n_raw rescue nil)],
+ ['PERCENT_TARGET_COVERAGE_1X',
+ (100.0 * stats[:coverage][1] / stats[:coverage][0] rescue nil)],
+ ['PERCENT_TARGET_COVERAGE_10X',
+ (100.0 * stats[:coverage][10] / stats[:coverage][0] rescue nil)],
+ ['PERCENT_TARGET_COVERAGE_20X',
+ (100.0 * stats[:coverage][20] / stats[:coverage][0] rescue nil)],
+ ['PERCENT_TARGET_COVERAGE_50X',
+ (100.0 * stats[:coverage][50] / stats[:coverage][0] rescue nil)],
+ ['PERCENT_TARGET_COVERAGE_100X',
+ (100.0 * stats[:coverage][100] / stats[:coverage][0] rescue nil)]]
+
+ stats
+ end
+end
projects / arvados.git / blobdiff