jpayne@68: #!/usr/bin/env python jpayne@68: # -*- coding: utf-8 -*- jpayne@68: jpayne@68: """ jpayne@68: stand alone read qc script based on jgi-rqc-pipeline/readqc/readqc.py v8.3.6 jpayne@68: jpayne@68: Command jpayne@68: $ readqc.py -f FASTQ.GZ -o OUT_DIR --skip-blast -html jpayne@68: jpayne@68: Outputs jpayne@68: - normal QC outputs + index.html jpayne@68: jpayne@68: Created: March 15, 2018 jpayne@68: jpayne@68: Shijie Yao (syao@lbl.gov) jpayne@68: jpayne@68: Revision: jpayne@68: jpayne@68: """ jpayne@68: jpayne@68: ## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## libraries to use jpayne@68: import os jpayne@68: import sys jpayne@68: import argparse jpayne@68: import datetime jpayne@68: import shutil jpayne@68: jpayne@68: SRC_ROOT = os.path.abspath(os.path.dirname(os.path.abspath(__file__))) jpayne@68: sys.path.append(SRC_ROOT + "/lib") # common jpayne@68: jpayne@68: from readqc_constants import RQCReadQcConfig, RQCReadQc, ReadqcStats jpayne@68: from common import get_logger, get_status, append_rqc_stats, append_rqc_file, set_colors, get_subsample_rate,run_command jpayne@68: from readqc_utils import * jpayne@68: #from readqc_utils import checkpoint_step_wrapper, fast_subsample_fastq_sequences, write_unique_20_mers, illumina_read_gc jpayne@68: from rqc_fastq import get_working_read_length, read_count jpayne@68: from readqc_report import * jpayne@68: from os_utility import make_dir jpayne@68: from html_utility import html_tag, html_th, html_tr, html_link jpayne@68: from rqc_utility import get_dict_obj, pipeline_val jpayne@68: jpayne@68: VERSION = "1.0.0" jpayne@68: LOG_LEVEL = "DEBUG" jpayne@68: SCRIPT_NAME = __file__ jpayne@68: jpayne@68: PYDIR = os.path.abspath(os.path.dirname(__file__)) jpayne@68: BBDIR = os.path.join(PYDIR, os.path.pardir) jpayne@68: jpayne@68: color = {} jpayne@68: color = set_colors(color, True) jpayne@68: jpayne@68: """ jpayne@68: STEP1 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: """ jpayne@68: def do_fast_subsample_fastq_sequences(fastq, skipSubsampling, log): jpayne@68: log.info("\n\n%sSTEP1 - Subsampling reads <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: jpayne@68: statsFile = RQCReadQcConfig.CFG["stats_file"] jpayne@68: jpayne@68: status = "1_illumina_readqc_subsampling in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: log.info("1_illumina_readqc_subsampling in progress.") jpayne@68: jpayne@68: inputReadNum = 0 jpayne@68: sampleRate = 0.0 jpayne@68: jpayne@68: if not skipSubsampling: jpayne@68: # sampleRate = RQCReadQc.ILLUMINA_SAMPLE_PCTENTAGE ## 0.01 jpayne@68: inputReadNum = read_count(fastq) ## read count from the original fastq jpayne@68: assert inputReadNum > 0, "ERROR: invalid input fastq" jpayne@68: sampleRate = get_subsample_rate(inputReadNum) jpayne@68: log.info("Subsampling rate = %s", sampleRate) jpayne@68: else: jpayne@68: log.info("1_illumina_readqc_subsampling: skip subsampling. Use all the reads.") jpayne@68: sampleRate = 1.0 jpayne@68: jpayne@68: retCode = None jpayne@68: totalReadNum = 0 jpayne@68: firstSubsampledFastqFileName = "" jpayne@68: jpayne@68: sequnitFileName = os.path.basename(fastq) jpayne@68: sequnitFileName = sequnitFileName.replace(".fastq", "").replace(".gz", "") jpayne@68: jpayne@68: firstSubsampledFastqFileName = sequnitFileName + ".s" + str(sampleRate) + ".fastq" jpayne@68: jpayne@68: retCode, firstSubsampledFastqFileName, totalBaseCount, totalReadNum, subsampledReadNum, bIsPaired, readLength = fast_subsample_fastq_sequences(fastq, firstSubsampledFastqFileName, sampleRate, True, log) jpayne@68: jpayne@68: append_rqc_stats(statsFile, "SUBSAMPLE_RATE", sampleRate, log) jpayne@68: jpayne@68: if retCode in (RQCExitCodes.JGI_FAILURE, -2): jpayne@68: log.info("1_illumina_readqc_subsampling failed.") jpayne@68: status = "1_illumina_readqc_subsampling failed" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: else: jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_BASE_COUNT, totalBaseCount, log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_COUNT, totalReadNum, log) jpayne@68: jpayne@68: log.info("1_illumina_readqc_subsampling complete.") jpayne@68: status = "1_illumina_readqc_subsampling complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: jpayne@68: return status, firstSubsampledFastqFileName, totalReadNum, subsampledReadNum, bIsPaired, readLength jpayne@68: jpayne@68: jpayne@68: """ jpayne@68: STEP2 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: """ jpayne@68: def do_write_unique_20_mers(fastq, totalReadCount, log): jpayne@68: log.info("\n\n%sSTEP2 - Sampling unique 25 mers <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: jpayne@68: if totalReadCount >= RQCReadQc.ILLUMINA_MER_SAMPLE_REPORT_FRQ * 2: ## 25000 jpayne@68: log.debug("read count total in step2 = %s", totalReadCount) jpayne@68: filesFile = RQCReadQcConfig.CFG["files_file"] jpayne@68: statsFile = RQCReadQcConfig.CFG["stats_file"] jpayne@68: jpayne@68: status = "2_unique_mers_sampling in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: log.info(status) jpayne@68: jpayne@68: retCode, newDataFile, newPngPlotFile, newHtmlPlotFile = write_unique_20_mers(fastq, log) jpayne@68: jpayne@68: if retCode != RQCExitCodes.JGI_SUCCESS: jpayne@68: status = "2_unique_mers_sampling failed" jpayne@68: log.error(status) jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: else: jpayne@68: ## if no output files, skip this step. jpayne@68: if newDataFile is not None: jpayne@68: statsDict = {} jpayne@68: jpayne@68: ## in readqc_report.py jpayne@68: ## 2014.07.23 read_level_mer_sampling is updated to process new output file format from bbcountunique jpayne@68: log.info("2_unique_mers_sampling: post-processing the bbcountunique output file.") jpayne@68: read_level_mer_sampling(statsDict, newDataFile, log) jpayne@68: jpayne@68: for k, v in statsDict.items(): jpayne@68: append_rqc_stats(statsFile, k, str(v), log) jpayne@68: jpayne@68: ## outputs from bbcountunique jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_20MER_UNIQUENESS_TEXT, newDataFile, log) jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_20MER_UNIQUENESS_PLOT, newPngPlotFile, log) jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_20MER_UNIQUENESS_D3_HTML_PLOT, newHtmlPlotFile, log) jpayne@68: jpayne@68: log.info("2_unique_mers_sampling complete.") jpayne@68: status = "2_unique_mers_sampling complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: else: jpayne@68: ## if num reads < RQCReadQc.ILLUMINA_MER_SAMPLE_REPORT_FRQ = 25000 jpayne@68: ## just proceed to the next step jpayne@68: log.warning("2_unique_mers_sampling can't run it because the number of reads < %s.", RQCReadQc.ILLUMINA_MER_SAMPLE_REPORT_FRQ * 2) jpayne@68: status = "2_unique_mers_sampling complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: jpayne@68: return status jpayne@68: jpayne@68: jpayne@68: """ jpayne@68: STEP3 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: """ jpayne@68: def do_illumina_read_gc(fastq, log): jpayne@68: log.info("\n\n%sSTEP3 - Making read GC histograms <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: jpayne@68: filesFile = RQCReadQcConfig.CFG["files_file"] jpayne@68: statsFile = RQCReadQcConfig.CFG["stats_file"] jpayne@68: jpayne@68: status = "3_illumina_read_gc in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: log.info("3_illumina_read_gc in progress.") jpayne@68: jpayne@68: reformat_gchist_file, png_file, htmlFile, mean_val, stdev_val, med_val, mode_val = illumina_read_gc(fastq, log) jpayne@68: jpayne@68: if not reformat_gchist_file: jpayne@68: log.error("3_illumina_read_gc failed.") jpayne@68: status = "3_illumina_read_gc failed" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: else: jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_GC_MEAN, mean_val, log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_GC_STD, stdev_val, log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_GC_MED, med_val, log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_GC_MODE, mode_val, log) jpayne@68: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_GC_TEXT, reformat_gchist_file, log) jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_GC_PLOT, png_file, log) jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_GC_D3_HTML_PLOT, htmlFile, log) jpayne@68: jpayne@68: log.info("3_illumina_read_gc complete.") jpayne@68: status = "3_illumina_read_gc complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: jpayne@68: return status jpayne@68: jpayne@68: jpayne@68: """ jpayne@68: STEP4 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: """ jpayne@68: def do_read_quality_stats(fastq, log): jpayne@68: log.info("\n\n%sSTEP4 - Analyzing read quality <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: jpayne@68: filesFile = RQCReadQcConfig.CFG["files_file"] jpayne@68: statsFile = RQCReadQcConfig.CFG["stats_file"] jpayne@68: jpayne@68: status = "4_illumina_read_quality_stats in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: log.info("4_illumina_read_quality_stats in progress.") jpayne@68: jpayne@68: readLength = 0 jpayne@68: readLenR1 = 0 jpayne@68: readLenR2 = 0 jpayne@68: isPairedEnd = None jpayne@68: jpayne@68: if not os.path.isfile(fastq): jpayne@68: log.error("4_illumina_read_quality_stats failed. Cannot find the input fastq file") jpayne@68: status = "4_illumina_read_quality_stats failed" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: return status jpayne@68: jpayne@68: ## First figure out if it's pair-ended or not jpayne@68: ## NOTE: ssize=10000 is recommended! jpayne@68: readLength, readLenR1, readLenR2, isPairedEnd = get_working_read_length(fastq, log) jpayne@68: log.info("Read length = %s, and is_pair_ended = %s", readLength, isPairedEnd) jpayne@68: jpayne@68: if readLength == 0: jpayne@68: log.error("Failed to run get_working_read_length.") jpayne@68: status = "4_illumina_read_quality_stats failed" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: else: jpayne@68: ## Pair-ended jpayne@68: r1_r2_baseposqual_png = None ## Average Base Position Quality Plot (*.qrpt.png) jpayne@68: r1_r2_baseposqual_html = None ## Average Base Position Quality D3 Plot (*.qrpt.html) jpayne@68: r1_r2_baseposqual_txt = None ## Read 1/2 Average Base Position Quality Text (*.qhist.txt) jpayne@68: jpayne@68: r1_baseposqual_box_png = None ## Read 1 Average Base Position Quality Boxplot (*.r1.png) jpayne@68: r2_baseposqual_box_png = None ## Read 2 Average Base Position Quality Boxplot (*.r2.png) jpayne@68: r1_baseposqual_box_html = None ## Read 1 Average Base Position Quality D3 Boxplot (*.r1.html) jpayne@68: r2_baseposqual_box_html = None ## Read 2 Average Base Position Quality D3 Boxplot (*.r2.html) jpayne@68: r1_r2_baseposqual_box_txt = None ## Average Base Position Quality text jpayne@68: jpayne@68: r1_cyclenbase_png = None ## Read 1 Percent N by Read Position (*.r1.fastq.base.stats.Npercent.png) --> Read 1 Cycle N Base Percent plot jpayne@68: r2_cyclenbase_png = None ## Read 2 Percent N by Read Position (*.r2.fastq.base.stats.Npercent.png) --> Read 2 Cycle N Base Percent plot jpayne@68: jpayne@68: r1_cyclenbase_txt = None ## Read 1 Percent N by Read Position Text (*.r1.fastq.base.stats) --> Read 1 Cycle N Base Percent text jpayne@68: #r2_cyclenbase_txt = None ## Read 2 Percent N by Read Position Text (*.r2.fastq.base.stats) --> Read 2 Cycle N Base Percent text jpayne@68: r1_r2_cyclenbase_txt = None ## Merged Percent N by Read Position Text (*.r2.fastq.base.stats) --> Merged Cycle N Base Percent text jpayne@68: jpayne@68: r1_cyclenbase_html = None jpayne@68: r2_cyclenbase_html = None jpayne@68: jpayne@68: r1_nuclcompfreq_png = None ## Read 1 Nucleotide Composition Frequency Plot (*.r1.stats.png) jpayne@68: r2_nuclcompfreq_png = None ## Read 2 Nucleotide Composition Frequency Plot (*.r2.stats.png) jpayne@68: r1_nuclcompfreq_html = None jpayne@68: r2_nuclcompfreq_html = None jpayne@68: jpayne@68: ## Single-ended jpayne@68: #se_baseposqual_txt = None ## Average Base Position Quality Text (*.qrpt) jpayne@68: #se_nuclcompfreq_png = None ## Cycle Nucleotide Composition (*.stats.png) jpayne@68: jpayne@68: if isPairedEnd: jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_LENGTH_1, readLenR1, log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_LENGTH_2, readLenR2, log) jpayne@68: jpayne@68: else: jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_LENGTH_1, readLength, log) jpayne@68: jpayne@68: ## Average Base Position Quality Plot/Text using qhist.txt jpayne@68: r1_r2_baseposqual_txt, r1_r2_baseposqual_png, r1_r2_baseposqual_html = gen_average_base_position_quality_plot(fastq, isPairedEnd, log) ## .reformat.qhist.txt jpayne@68: jpayne@68: log.debug("Outputs: %s %s %s", r1_r2_baseposqual_png, r1_r2_baseposqual_html, r1_r2_baseposqual_txt) jpayne@68: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_QUAL_POS_PLOT_MERGED, r1_r2_baseposqual_png, log) jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_QUAL_POS_PLOT_MERGED_D3_HTML_PLOT, r1_r2_baseposqual_html, log) jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_QUAL_POS_QRPT_1, r1_r2_baseposqual_txt, log) ## for backward compatibility jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_QUAL_POS_QRPT_2, r1_r2_baseposqual_txt, log) ## for backward compatibility jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_QUAL_POS_QRPT_MERGED, r1_r2_baseposqual_txt, log) jpayne@68: jpayne@68: ## Average Base Position Quality Plot/Text using bqhist.txt jpayne@68: r1_r2_baseposqual_box_txt, r1_baseposqual_box_png, r2_baseposqual_box_png, r1_baseposqual_box_html, r2_baseposqual_box_html = gen_average_base_position_quality_boxplot(fastq, log) jpayne@68: jpayne@68: log.debug("Read qual outputs: %s %s %s %s %s", r1_r2_baseposqual_box_txt, r1_baseposqual_box_png, r1_baseposqual_box_html, r2_baseposqual_box_png, r2_baseposqual_box_html) jpayne@68: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_QUAL_POS_QRPT_BOXPLOT_1, r1_baseposqual_box_png, log) jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_QUAL_POS_QRPT_D3_HTML_BOXPLOT_1, r1_baseposqual_box_html, log) jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_QUAL_POS_QRPT_BOXPLOT_2, r2_baseposqual_box_png, log) jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_QUAL_POS_QRPT_D3_HTML_BOXPLOT_2, r2_baseposqual_box_html, log) jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_QUAL_POS_QRPT_BOXPLOT_TEXT, r1_r2_baseposqual_box_txt, log) jpayne@68: jpayne@68: ## ---------------------------------------------------------------------------------------------------- jpayne@68: ## compute Q20 of the two reads jpayne@68: q20Read1 = None jpayne@68: q20Read2 = None jpayne@68: jpayne@68: ## using bqhist.txt jpayne@68: if r1_r2_baseposqual_box_txt: jpayne@68: q20Read1 = q20_score_new(r1_r2_baseposqual_box_txt, 1, log) jpayne@68: if isPairedEnd: jpayne@68: q20Read2 = q20_score_new(r1_r2_baseposqual_box_txt, 2, log) jpayne@68: jpayne@68: log.debug("q20 for read 1 = %s", q20Read1) jpayne@68: log.debug("q20 for read 2 = %s", q20Read2) jpayne@68: jpayne@68: if q20Read1 is not None: jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_Q20_READ1, q20Read1, log) jpayne@68: else: jpayne@68: log.error("Failed to get q20 read 1 from %s", r1_r2_baseposqual_box_txt) jpayne@68: status = "4_illumina_read_quality_stats failed" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: return status jpayne@68: jpayne@68: if isPairedEnd: jpayne@68: if q20Read2 is not None: jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_Q20_READ2, q20Read2, log) jpayne@68: else: jpayne@68: log.error("Failed to get q20 read 2 from %s", r1_r2_baseposqual_box_txt) jpayne@68: status = "4_illumina_read_quality_stats failed" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: return status jpayne@68: jpayne@68: r1_r2_cyclenbase_txt, r1_nuclcompfreq_png, r1_nuclcompfreq_html, r2_nuclcompfreq_png, r2_nuclcompfreq_html = gen_cycle_nucleotide_composition_plot(fastq, readLength, isPairedEnd, log) jpayne@68: jpayne@68: log.debug("gen_cycle_nucleotide_composition_plot() ==> %s %s %s %s %s", r1_cyclenbase_txt, r1_nuclcompfreq_png, r1_nuclcompfreq_html, r2_nuclcompfreq_png, r2_nuclcompfreq_html) jpayne@68: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_BASE_COUNT_TEXT_1, r1_r2_cyclenbase_txt, log) jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_BASE_COUNT_TEXT_2, r1_r2_cyclenbase_txt, log) # reformat.sh generates a single merged output file. jpayne@68: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_BASE_COUNT_PLOT_1, r1_nuclcompfreq_png, log) jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_BASE_COUNT_D3_HTML_PLOT_1, r1_nuclcompfreq_html, log) jpayne@68: jpayne@68: if r2_nuclcompfreq_png: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_BASE_COUNT_PLOT_2, r2_nuclcompfreq_png, log) jpayne@68: if r2_nuclcompfreq_html: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_BASE_COUNT_D3_HTML_PLOT_2, r2_nuclcompfreq_html, log) jpayne@68: jpayne@68: ### --------------------------------------------------------------------------------------------------- jpayne@68: ## using bhist.txt jpayne@68: r1_cyclenbase_txt, r1_cyclenbase_png, r1_cyclenbase_html, r2_cyclenbase_png, r2_cyclenbase_html = gen_cycle_n_base_percent_plot(fastq, readLength, isPairedEnd, log) jpayne@68: jpayne@68: log.debug("Outputs: %s %s %s", r1_cyclenbase_txt, r1_cyclenbase_png, r1_cyclenbase_html) jpayne@68: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_BASE_PERCENTAGE_TEXT_1, r1_cyclenbase_txt, log) jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_BASE_PERCENTAGE_TEXT_2, r1_cyclenbase_txt, log) jpayne@68: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_BASE_PERCENTAGE_PLOT_1, r1_cyclenbase_png, log) jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_BASE_PERCENTAGE_D3_HTML_PLOT_1, r1_cyclenbase_html, log) jpayne@68: jpayne@68: if r2_cyclenbase_png: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_BASE_PERCENTAGE_PLOT_2, r2_cyclenbase_png, log) jpayne@68: if r2_cyclenbase_html: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_BASE_PERCENTAGE_D3_HTML_PLOT_2, r2_cyclenbase_html, log) jpayne@68: jpayne@68: log.info("4_illumina_read_quality_stats complete.") jpayne@68: status = "4_illumina_read_quality_stats complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: jpayne@68: return status jpayne@68: jpayne@68: jpayne@68: jpayne@68: """ jpayne@68: STEP5 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: """ jpayne@68: def do_write_base_quality_stats(fastq, log): jpayne@68: log.info("\n\n%sSTEP5 - Calculating base quality statistics for reads <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: jpayne@68: filesFile = RQCReadQcConfig.CFG["files_file"] jpayne@68: statsFile = RQCReadQcConfig.CFG["stats_file"] jpayne@68: jpayne@68: status = "5_illumina_read_quality_stats in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: log.info("5_illumina_read_quality_stats in progress.") jpayne@68: jpayne@68: reformatObqhistFile = write_avg_base_quality_stats(fastq, log) ## *.reformat.obqhist.txt jpayne@68: jpayne@68: if not reformatObqhistFile: jpayne@68: log.error("5_illumina_read_quality_stats failed.") jpayne@68: status = "5_illumina_read_quality_stats failed" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: else: jpayne@68: ## Generate qual scores and plots of read level QC jpayne@68: statsDict = {} jpayne@68: jpayne@68: retCode = base_level_qual_stats(statsDict, reformatObqhistFile, log) jpayne@68: jpayne@68: if retCode != RQCExitCodes.JGI_SUCCESS: jpayne@68: log.error("5_illumina_read_quality_stats failed.") jpayne@68: status = "5_illumina_read_quality_stats failed" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: else: jpayne@68: for k, v in statsDict.items(): jpayne@68: append_rqc_stats(statsFile, k, str(v), log) jpayne@68: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_BASE_QUALITY_STATS, reformatObqhistFile, log) jpayne@68: jpayne@68: log.info("5_illumina_read_quality_stats complete.") jpayne@68: status = "5_illumina_read_quality_stats complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: jpayne@68: return status jpayne@68: jpayne@68: jpayne@68: jpayne@68: """ jpayne@68: STEP6 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: """ jpayne@68: def do_illumina_count_q_score(fastq, log): jpayne@68: log.info("\n\n%sSTEP6 - Generating quality score histogram <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: jpayne@68: filesFile = RQCReadQcConfig.CFG["files_file"] jpayne@68: statsFile = RQCReadQcConfig.CFG["stats_file"] jpayne@68: jpayne@68: status = "6_illumina_count_q_score in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: log.info("6_illumina_count_q_score in progress.") jpayne@68: jpayne@68: qhistTxtFile, qhistPngFile, qhistHtmlPlotFile = illumina_count_q_score(fastq, log) ## *.obqhist.txt jpayne@68: jpayne@68: if not qhistTxtFile: jpayne@68: log.error("6_illumina_count_q_score failed.") jpayne@68: status = "6_illumina_count_q_score failed" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: else: jpayne@68: ## save qscores in statsFile jpayne@68: qscore = {} jpayne@68: read_level_qual_stats(qscore, qhistTxtFile, log) jpayne@68: jpayne@68: for k, v in qscore.items(): jpayne@68: append_rqc_stats(statsFile, k, str(v), log) jpayne@68: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_QHIST_TEXT, qhistTxtFile, log) jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_QHIST_PLOT, qhistPngFile, log) jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_QHIST_D3_HTML_PLOT, qhistHtmlPlotFile, log) jpayne@68: jpayne@68: log.info("6_illumina_count_q_score complete.") jpayne@68: status = "6_illumina_count_q_score complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: jpayne@68: return status jpayne@68: jpayne@68: jpayne@68: jpayne@68: """ jpayne@68: STEP7 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: """ jpayne@68: ## 20140903 removed jpayne@68: ##def do_illumina_calculate_average_quality(fastq, log): jpayne@68: jpayne@68: jpayne@68: jpayne@68: """ jpayne@68: STEP8 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: """ jpayne@68: def do_illumina_find_common_motifs(fastq, log): jpayne@68: log.info("\n\n%sSTEP8 - Locating N stutter motifs in sequence reads <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: jpayne@68: filesFile = RQCReadQcConfig.CFG["files_file"] jpayne@68: statsFile = RQCReadQcConfig.CFG["stats_file"] jpayne@68: jpayne@68: status = "8_illumina_find_common_motifs in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: log.info(status) jpayne@68: jpayne@68: retCode, statDataFile = illumina_find_common_motifs(fastq, log) jpayne@68: jpayne@68: log.info("nstutter statDataFile name = %s", statDataFile) jpayne@68: jpayne@68: if retCode != RQCExitCodes.JGI_SUCCESS: jpayne@68: log.error("8_illumina_find_common_motifs failed.") jpayne@68: status = "8_illumina_find_common_motifs failed" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: else: jpayne@68: ## read_level_stutter jpayne@68: ## ex) jpayne@68: ##688 N---------------------------------------------------------------------------------------------------------------------------- jpayne@68: ##------------------------- jpayne@68: ##346 NNNNNNNNNNNNN---------------------------------------------------------------------------------------------------------------- jpayne@68: ##------------------------- jpayne@68: ##53924 ------------N---------------------------------------------------------------------------------------------------------------- jpayne@68: ##------------------------- jpayne@68: ##sum pct patterNs past 0.1 == 15.9245930330268 ( 54958 / 345114 * 100 ) jpayne@68: jpayne@68: with open(statDataFile, "r") as stutFH: jpayne@68: lines = stutFH.readlines() jpayne@68: jpayne@68: ## if no motifs are detected the file is empty jpayne@68: if not lines: jpayne@68: log.warning("The *.nstutter.stat file is not available in function read_level_stutter(). The function still returns JGI_SUCCESS.") jpayne@68: jpayne@68: else: jpayne@68: assert lines[-1].find("patterNs") != -1 jpayne@68: t = lines[-1].split() jpayne@68: ## ["sum", "pct", "patterNs", "past", "0.1", "==", "15.9245930330268", "(", "54958", "/", "345114", "*", "100", ")"] jpayne@68: percent = "%.2f" % float(t[6]) jpayne@68: jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_N_FREQUENCE, percent, log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_N_PATTERN, str("".join(lines[:-1])), log) jpayne@68: jpayne@68: ## NOTE: ??? jpayne@68: append_rqc_file(filesFile, "find_common_motifs.dataFile", statDataFile, log) jpayne@68: jpayne@68: log.info("8_illumina_find_common_motifs complete.") jpayne@68: status = "8_illumina_find_common_motifs complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: jpayne@68: return status jpayne@68: jpayne@68: jpayne@68: """ jpayne@68: STEP11 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: """ jpayne@68: def do_illumina_detect_read_contam(fastq, bpToCut, log): jpayne@68: log.info("\n\n%sSTEP11 - Detect read contam <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: jpayne@68: filesFile = RQCReadQcConfig.CFG["files_file"] jpayne@68: statsFile = RQCReadQcConfig.CFG["stats_file"] jpayne@68: jpayne@68: status = "11_illumina_detect_read_contam in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: log.info("11_illumina_detect_read_contam in progress.") jpayne@68: jpayne@68: jpayne@68: ######### jpayne@68: ## seal jpayne@68: ######### jpayne@68: retCode2, outFileDict2, ratioResultDict2, contamStatDict = illumina_detect_read_contam3(fastq, bpToCut, log) ## seal version jpayne@68: jpayne@68: if retCode2 != RQCExitCodes.JGI_SUCCESS: jpayne@68: log.error("11_illumina_detect_read_contam seal version failed.") jpayne@68: status = "11_illumina_detect_read_contam failed" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: else: jpayne@68: for k, v in outFileDict2.items(): jpayne@68: append_rqc_file(filesFile, k + " seal", str(v), log) jpayne@68: append_rqc_file(filesFile, k, str(v), log) jpayne@68: jpayne@68: for k, v in ratioResultDict2.items(): jpayne@68: append_rqc_stats(statsFile, k + " seal", str(v), log) jpayne@68: append_rqc_stats(statsFile, k, str(v), log) jpayne@68: jpayne@68: ## contamination stat jpayne@68: for k, v in contamStatDict.items(): jpayne@68: append_rqc_stats(statsFile, k, str(v), log) jpayne@68: jpayne@68: log.info("11_illumina_detect_read_contam seal version complete.") jpayne@68: status = "11_illumina_detect_read_contam complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: jpayne@68: return status jpayne@68: jpayne@68: jpayne@68: jpayne@68: """ jpayne@68: STEP13 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: """ jpayne@68: ## Removed!! jpayne@68: ##def do_illumina_read_megablast(firstSubsampledFastqFileName, skipSubsampling, subsampledReadNum, log, blastDbPath=None): jpayne@68: jpayne@68: jpayne@68: jpayne@68: """ jpayne@68: New STEP13 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: """ jpayne@68: def do_illumina_subsampling_read_blastn(firstSubsampledFastqFileName, skipSubsampling, subsampledReadNum, log): jpayne@68: log.info("\n\n%sSTEP13 - Run subsampling for Blast search <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: jpayne@68: subsampeldFastqFile = None jpayne@68: totalReadNum = 0 jpayne@68: readNumToReturn = 0 jpayne@68: jpayne@68: status = "13_illumina_subsampling_read_megablast in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: log.info("13_illumina_subsampling_read_megablast in progress.") jpayne@68: jpayne@68: if subsampledReadNum == 0: jpayne@68: cmd = " ".join(["grep", "-c", "'^+'", firstSubsampledFastqFileName]) jpayne@68: jpayne@68: stdOut, _, exitCode = run_command(cmd, True, log) jpayne@68: jpayne@68: if exitCode != 0: jpayne@68: log.error("Failed to run grep cmd") jpayne@68: return RQCExitCodes.JGI_FAILURE, None, None, None jpayne@68: jpayne@68: else: jpayne@68: readNum = int(stdOut) jpayne@68: jpayne@68: else: jpayne@68: readNum = subsampledReadNum jpayne@68: jpayne@68: log.info("Subsampled read number = %d", readNum) jpayne@68: jpayne@68: sampl_per = RQCReadQc.ILLUMINA_SAMPLE_PCTENTAGE ## 0.01 jpayne@68: max_count = RQCReadQc.ILLUMINA_SAMPLE_COUNT ## 50000 jpayne@68: jpayne@68: ## TODO jpayne@68: ## Use "samplereadstarget" option in reformat.sh jpayne@68: jpayne@68: if skipSubsampling: jpayne@68: log.debug("No subsampling for megablast. Use fastq file, %s (readnum = %s) as query for megablast.", firstSubsampledFastqFileName, readNum) jpayne@68: subsampeldFastqFile = firstSubsampledFastqFileName jpayne@68: readNumToReturn = readNum jpayne@68: jpayne@68: else: jpayne@68: if readNum > max_count: jpayne@68: log.debug("Run the 2nd subsampling for running megablast.") jpayne@68: secondSubsamplingRate = float(max_count) / readNum jpayne@68: log.debug("SecondSubsamplingRate=%s, max_count=%s, readNum=%s", secondSubsamplingRate, max_count, readNum) jpayne@68: log.info("Second subsampling of Reads after Percent Subsampling reads = %s with new percent subsampling %f.", readNum, secondSubsamplingRate) jpayne@68: jpayne@68: secondSubsampledFastqFile = "" jpayne@68: # dataFile = "" jpayne@68: jpayne@68: sequnitFileName = os.path.basename(firstSubsampledFastqFileName) jpayne@68: sequnitFileName = sequnitFileName.replace(".fastq", "").replace(".gz", "") jpayne@68: jpayne@68: secondSubsampledFastqFile = sequnitFileName + ".s" + str(sampl_per) + ".s" + str(secondSubsamplingRate) + ".n" + str(max_count) + ".fastq" jpayne@68: jpayne@68: ## ex) fq_sub_sample.pl -f .../7601.1.77813.CTTGTA.s0.01.fastq -o .../7601.1.77813.CTTGTA.s0.01.stats -r 0.0588142575171 > .../7601.1.77813.CTTGTA.s0.01.s0.01.s0.0588142575171.n50000.fastq jpayne@68: ## ex) reformat.sh jpayne@68: retCode, secondSubsampledFastqFile, totalBaseCount, totalReadNum, subsampledReadNum, _, _ = fast_subsample_fastq_sequences(firstSubsampledFastqFileName, secondSubsampledFastqFile, secondSubsamplingRate, False, log) jpayne@68: jpayne@68: if retCode != RQCExitCodes.JGI_SUCCESS: jpayne@68: log.error("Second subsampling failed.") jpayne@68: return RQCExitCodes.JGI_FAILURE, None, None, None jpayne@68: jpayne@68: else: jpayne@68: log.info("Second subsampling complete.") jpayne@68: log.info("Second Subsampling Total Base Count = %s.", totalBaseCount) jpayne@68: log.info("Second Subsampling Total Reads = %s.", totalReadNum) jpayne@68: log.info("Second Subsampling Sampled Reads = %s.", subsampledReadNum) jpayne@68: jpayne@68: if subsampledReadNum == 0: jpayne@68: log.warning("Too small first subsampled fastq file. Skip the 2nd sampling.") jpayne@68: secondSubsampledFastqFile = firstSubsampledFastqFileName jpayne@68: jpayne@68: subsampeldFastqFile = secondSubsampledFastqFile jpayne@68: readNumToReturn = subsampledReadNum jpayne@68: jpayne@68: else: jpayne@68: log.debug("The readNum is smaller than max_count=%s. The 2nd sampling skipped.", str(max_count)) jpayne@68: subsampeldFastqFile = firstSubsampledFastqFileName jpayne@68: totalReadNum = readNum jpayne@68: readNumToReturn = readNum jpayne@68: jpayne@68: log.info("13_illumina_subsampling_read_megablast complete.") jpayne@68: status = "13_illumina_subsampling_read_megablast complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: jpayne@68: return status, subsampeldFastqFile, totalReadNum, readNumToReturn jpayne@68: jpayne@68: jpayne@68: jpayne@68: ##""" jpayne@68: ##STEP14 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## jpayne@68: ##""" jpayne@68: ##def do_illumina_read_blastn_refseq_microbial(subsampeldFastqFile, subsampledReadNum, log, blastDbPath=None): jpayne@68: ## 12212015 sulsj REMOVED! jpayne@68: jpayne@68: jpayne@68: """ jpayne@68: New STEP14 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: """ jpayne@68: def do_illumina_read_blastn_refseq_archaea(subsampeldFastqFile, subsampledReadNum, log): jpayne@68: log.info("\n\n%sSTEP14 - Run read blastn against refseq archaea <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: jpayne@68: statsFile = RQCReadQcConfig.CFG["stats_file"] jpayne@68: retCode = None jpayne@68: jpayne@68: log.info("14_illumina_read_blastn_refseq_archaea in progress.") jpayne@68: status = "14_illumina_read_blastn_refseq_archaea in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: retCode = illumina_read_blastn_refseq_archaea(subsampeldFastqFile, log) jpayne@68: jpayne@68: if retCode == RQCExitCodes.JGI_FAILURE: jpayne@68: log.error("14_illumina_read_blastn_refseq_archaea failed.") jpayne@68: status = "14_illumina_read_blastn_refseq_archaea failed" jpayne@68: jpayne@68: elif retCode == -143: ## timeout jpayne@68: log.warning("14_illumina_read_blastn_refseq_archaea timeout.") jpayne@68: status = "14_illumina_read_blastn_refseq_archaea complete" jpayne@68: jpayne@68: else: jpayne@68: ## read number used in blast search? jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READS_NUMBER, subsampledReadNum, log) jpayne@68: jpayne@68: ret2 = read_megablast_hits("refseq.archaea", log) jpayne@68: jpayne@68: if ret2 != RQCExitCodes.JGI_SUCCESS: jpayne@68: log.error("Errors in read_megablast_hits() of refseq.microbial") jpayne@68: log.error("14_illumina_read_blastn_refseq_archaea reporting failed.") jpayne@68: status = "14_illumina_read_blastn_refseq_archaea failed" jpayne@68: jpayne@68: else: jpayne@68: log.info("14_illumina_read_blastn_refseq_archaea complete.") jpayne@68: status = "14_illumina_read_blastn_refseq_archaea complete" jpayne@68: jpayne@68: jpayne@68: return status jpayne@68: jpayne@68: jpayne@68: """ jpayne@68: STEP15 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: """ jpayne@68: def do_illumina_read_blastn_refseq_bacteria(subsampeldFastqFile, log): jpayne@68: log.info("\n\n%sSTEP15 - Run read blastn against refseq bacteria <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: jpayne@68: #statsFile = RQCReadQcConfig.CFG["stats_file"] jpayne@68: #retCode = None jpayne@68: jpayne@68: log.info("15_illumina_read_blastn_refseq_bacteria in progress.") jpayne@68: status = "15_illumina_read_blastn_refseq_bacteria in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: retCode = illumina_read_blastn_refseq_bacteria(subsampeldFastqFile, log) jpayne@68: jpayne@68: if retCode == RQCExitCodes.JGI_FAILURE: jpayne@68: log.error("15_illumina_read_blastn_refseq_bacteria failed.") jpayne@68: status = "15_illumina_read_blastn_refseq_bacteria failed" jpayne@68: jpayne@68: elif retCode == -143: ## timeout jpayne@68: log.warning("15_illumina_read_blastn_refseq_bacteria timeout.") jpayne@68: status = "15_illumina_read_blastn_refseq_bacteria complete" jpayne@68: jpayne@68: else: jpayne@68: ret2 = read_megablast_hits("refseq.bacteria", log) jpayne@68: jpayne@68: if ret2 != RQCExitCodes.JGI_SUCCESS: jpayne@68: log.error("Errors in read_megablast_hits() of refseq.microbial") jpayne@68: log.error("15_illumina_read_blastn_refseq_bacteria reporting failed.") jpayne@68: status = "15_illumina_read_blastn_refseq_bacteria failed" jpayne@68: jpayne@68: else: jpayne@68: log.info("15_illumina_read_blastn_refseq_bacteria complete.") jpayne@68: status = "15_illumina_read_blastn_refseq_bacteria complete" jpayne@68: jpayne@68: jpayne@68: return status jpayne@68: jpayne@68: jpayne@68: jpayne@68: """ jpayne@68: STEP16 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: """ jpayne@68: def do_illumina_read_blastn_nt(subsampeldFastqFile, log): jpayne@68: log.info("\n\n%sSTEP16 - Run read blastn against nt <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: jpayne@68: retCode = None jpayne@68: jpayne@68: log.info("16_illumina_read_blastn_nt in progress.") jpayne@68: status = "16_illumina_read_blastn_nt in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: retCode = illumina_read_blastn_nt(subsampeldFastqFile, log) jpayne@68: jpayne@68: if retCode == RQCExitCodes.JGI_FAILURE: jpayne@68: log.error("16_illumina_read_blastn_nt failed.") jpayne@68: status = "16_illumina_read_blastn_nt failed" jpayne@68: jpayne@68: elif retCode == -143: ## timeout jpayne@68: log.warning("16_illumina_read_blastn_nt timeout.") jpayne@68: status = "16_illumina_read_blastn_nt complete" jpayne@68: jpayne@68: else: jpayne@68: ret2 = read_megablast_hits("nt", log) jpayne@68: jpayne@68: if ret2 != RQCExitCodes.JGI_SUCCESS: jpayne@68: log.error("Errors in read_megablast_hits() of nt") jpayne@68: log.error("16_illumina_read_blastn_nt reporting failed.") jpayne@68: status = "16_illumina_read_blastn_nt failed" jpayne@68: jpayne@68: else: jpayne@68: log.info("16_illumina_read_blastn_nt complete.") jpayne@68: status = "16_illumina_read_blastn_nt complete" jpayne@68: jpayne@68: jpayne@68: return status jpayne@68: jpayne@68: jpayne@68: jpayne@68: """ jpayne@68: STEP17 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: do_illumina_multiplex_statistics jpayne@68: jpayne@68: Demultiplexing analysis for pooled lib jpayne@68: jpayne@68: """ jpayne@68: def do_illumina_multiplex_statistics(fastq, log, isMultiplexed=None): jpayne@68: log.info("\n\n%sSTEP17 - Run Multiplex statistics analysis <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: jpayne@68: filesFile = RQCReadQcConfig.CFG["files_file"] jpayne@68: jpayne@68: log.info("17_multiplex_statistics in progress.") jpayne@68: status = "17_multiplex_statistics in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: log.debug("fastq file: %s", fastq) jpayne@68: jpayne@68: retCode, demultiplexStatsFile, detectionPngPlotFile, detectionHtmlPlotFile = illumina_generate_index_sequence_detection_plot(fastq, log, isMultiplexed=isMultiplexed) jpayne@68: jpayne@68: if retCode != RQCExitCodes.JGI_SUCCESS: jpayne@68: log.error("17_multiplex_statistics failed.") jpayne@68: status = "17_multiplex_statistics failed" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: else: jpayne@68: log.info("17_multiplex_statistics complete.") jpayne@68: status = "17_multiplex_statistics complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: if detectionPngPlotFile is not None: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_DEMULTIPLEX_STATS_PLOT, detectionPngPlotFile, log) jpayne@68: jpayne@68: if detectionHtmlPlotFile is not None: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_DEMULTIPLEX_STATS_D3_HTML_PLOT, detectionHtmlPlotFile, log) jpayne@68: jpayne@68: if demultiplexStatsFile is not None: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_DEMULTIPLEX_STATS, demultiplexStatsFile, log) jpayne@68: jpayne@68: jpayne@68: return status jpayne@68: jpayne@68: jpayne@68: jpayne@68: """ jpayne@68: STEP18 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: do_end_of_read_illumina_adapter_check jpayne@68: jpayne@68: """ jpayne@68: def do_end_of_read_illumina_adapter_check(firstSubsampledFastqFileName, log): jpayne@68: log.info("\n\n%sSTEP18 - Run end_of_read_illumina_adapter_check analysis <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: jpayne@68: filesFile = RQCReadQcConfig.CFG["files_file"] jpayne@68: plotFile = None jpayne@68: dataFile = None jpayne@68: jpayne@68: log.info("18_end_of_read_illumina_adapter_check in progress.") jpayne@68: status = "18_end_of_read_illumina_adapter_check in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: log.debug("sampled fastq file: %s", firstSubsampledFastqFileName) jpayne@68: jpayne@68: retCode, dataFile, plotFile, htmlFile = end_of_read_illumina_adapter_check(firstSubsampledFastqFileName, log) jpayne@68: jpayne@68: if retCode != RQCExitCodes.JGI_SUCCESS: jpayne@68: log.error("18_end_of_read_illumina_adapter_check failed.") jpayne@68: status = "18_end_of_read_illumina_adapter_check failed" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: else: jpayne@68: log.info("18_end_of_read_illumina_adapter_check complete.") jpayne@68: status = "18_end_of_read_illumina_adapter_check complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: if plotFile is not None: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_END_OF_READ_ADAPTER_CHECK_PLOT, plotFile, log) jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_END_OF_READ_ADAPTER_CHECK_D3_HTML_PLOT, htmlFile, log) jpayne@68: jpayne@68: if dataFile is not None: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_END_OF_READ_ADAPTER_CHECK_DATA, dataFile, log) jpayne@68: jpayne@68: jpayne@68: return status jpayne@68: jpayne@68: jpayne@68: jpayne@68: """ jpayne@68: STEP19 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: do_insert_size_analysis jpayne@68: jpayne@68: """ jpayne@68: def do_insert_size_analysis(fastq, log): jpayne@68: log.info("\n\n%sSTEP19 - Run insert size analysis <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: jpayne@68: filesFile = RQCReadQcConfig.CFG["files_file"] jpayne@68: statsFile = RQCReadQcConfig.CFG["stats_file"] jpayne@68: jpayne@68: plotFile = None jpayne@68: dataFile = None jpayne@68: jpayne@68: log.info("19_insert_size_analysis in progress.") jpayne@68: status = "19_insert_size_analysis in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: log.debug("fastq file used: %s", fastq) jpayne@68: jpayne@68: retCode, dataFile, plotFile, htmlFile, statsDict = insert_size_analysis(fastq, log) ## by bbmerge.sh jpayne@68: jpayne@68: if retCode != RQCExitCodes.JGI_SUCCESS: jpayne@68: log.error("19_insert_size_analysis failed.") jpayne@68: status = "19_insert_size_analysis failed" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: else: jpayne@68: if plotFile is not None: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_HISTO_PLOT, plotFile, log) jpayne@68: jpayne@68: if dataFile is not None: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_HISTO_DATA, dataFile, log) jpayne@68: jpayne@68: if htmlFile is not None: jpayne@68: append_rqc_file(filesFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_HISTO_D3_HTML_PLOT, htmlFile, log) jpayne@68: jpayne@68: if statsDict: jpayne@68: try: jpayne@68: ## -------------------------------------------------------------------------------------------------------- jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_AVG_INSERT, statsDict["avg_insert"], log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_STD_INSERT, statsDict["std_insert"], log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_MODE_INSERT, statsDict["mode_insert"], log) jpayne@68: ## -------------------------------------------------------------------------------------------------------- jpayne@68: jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_TOTAL_TIME, statsDict["total_time"], log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_NUM_READS, statsDict["num_reads"], log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_JOINED_NUM, statsDict["joined_num"], log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_JOINED_PERC, statsDict["joined_perc"], log) jpayne@68: jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_AMBIGUOUS_NUM, statsDict["ambiguous_num"], log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_AMBIGUOUS_PERC, statsDict["ambiguous_perc"], log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_NO_SOLUTION_NUM, statsDict["no_solution_num"], log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_NO_SOLUTION_PERC, statsDict["no_solution_perc"], log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_TOO_SHORT_NUM, statsDict["too_short_num"], log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_TOO_SHORT_PERC, statsDict["too_short_perc"], log) jpayne@68: jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_INSERT_RANGE_START, statsDict["insert_range_start"], log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_INSERT_RANGE_END, statsDict["insert_range_end"], log) jpayne@68: jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_90TH_PERC, statsDict["perc_90th"], log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_50TH_PERC, statsDict["perc_50th"], log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_10TH_PERC, statsDict["perc_10th"], log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_75TH_PERC, statsDict["perc_75th"], log) jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_READ_INSERT_SIZE_25TH_PERC, statsDict["perc_25th"], log) jpayne@68: jpayne@68: except KeyError: jpayne@68: log.error("19_insert_size_analysis failed (KeyError).") jpayne@68: status = "19_insert_size_analysis failed" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: return status jpayne@68: jpayne@68: log.info("19_insert_size_analysis complete.") jpayne@68: status = "19_insert_size_analysis complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: jpayne@68: return status jpayne@68: jpayne@68: jpayne@68: jpayne@68: jpayne@68: """ jpayne@68: STEP21 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: do_sketch_vs_nt_refseq_silva jpayne@68: jpayne@68: """ jpayne@68: def do_sketch_vs_nt_refseq_silva(fasta, log): jpayne@68: log.info("\n\n%sSTEP21 - Run sketch vs nt, refseq, silva <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: jpayne@68: status = "21_sketch_vs_nt_refseq_silva in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: sketchOutDir = "sketch" jpayne@68: sketchOutPath = os.path.join(outputPath, sketchOutDir) jpayne@68: make_dir(sketchOutPath) jpayne@68: change_mod(sketchOutPath, "0755") jpayne@68: jpayne@68: seqUnitName = os.path.basename(fasta) jpayne@68: seqUnitName = file_name_trim(seqUnitName) jpayne@68: jpayne@68: filesFile = RQCReadQcConfig.CFG["files_file"] jpayne@68: jpayne@68: comOptions = "ow=t colors=f printtaxa=t depth depth2 unique2 merge" jpayne@68: jpayne@68: ## NT ########################## jpayne@68: sketchOutFile = os.path.join(sketchOutPath, seqUnitName + ".sketch_vs_nt.txt") jpayne@68: sendSketchShCmd = os.path.join(BBDIR, 'sendsketch.sh') jpayne@68: cmd = "%s in=%s out=%s %s nt" % (sendSketchShCmd, fasta, sketchOutFile, comOptions) jpayne@68: stdOut, stdErr, exitCode = run_sh_command(cmd, True, log, True) jpayne@68: if exitCode != 0: jpayne@68: log.error("Failed to run : %s, stdout : %s, stderr: %s", cmd, stdOut, stdErr) jpayne@68: status = "21_sketch_vs_nt_refseq_silva failed" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: return status jpayne@68: jpayne@68: append_rqc_file(filesFile, "sketch_vs_nt_output", sketchOutFile, log) jpayne@68: jpayne@68: ## Refseq ########################## jpayne@68: sketchOutFile = os.path.join(sketchOutPath, seqUnitName + ".sketch_vs_refseq.txt") jpayne@68: cmd = "%s in=%s out=%s %s refseq" % (sendSketchShCmd, fasta, sketchOutFile, comOptions) jpayne@68: stdOut, stdErr, exitCode = run_sh_command(cmd, True, log, True) jpayne@68: if exitCode != 0: jpayne@68: log.error("Failed to run : %s, stdout : %s, stderr: %s", cmd, stdOut, stdErr) jpayne@68: status = "21_sketch_vs_refseq failed" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: return status jpayne@68: jpayne@68: append_rqc_file(filesFile, "sketch_vs_refseq_output", sketchOutFile, log) jpayne@68: jpayne@68: ## Silva ########################## jpayne@68: sketchOutFile = os.path.join(sketchOutPath, seqUnitName + ".sketch_vs_silva.txt") jpayne@68: cmd = "%s in=%s out=%s %s silva" % (sendSketchShCmd, fasta, sketchOutFile, comOptions) jpayne@68: stdOut, stdErr, exitCode = run_sh_command(cmd, True, log, True) jpayne@68: if exitCode != 0: jpayne@68: log.error("Failed to run : %s, stdout : %s, stderr: %s", cmd, stdOut, stdErr) jpayne@68: status = "21_sketch_vs_silva failed" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: return status jpayne@68: jpayne@68: append_rqc_file(filesFile, "sketch_vs_silva_output", sketchOutFile, log) jpayne@68: jpayne@68: log.info("21_sketch_vs_nt_refseq_silva complete.") jpayne@68: status = "21_sketch_vs_nt_refseq_silva complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: return status jpayne@68: jpayne@68: jpayne@68: jpayne@68: """ jpayne@68: STEP22 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: do_illumina_read_level_report_postprocessing jpayne@68: jpayne@68: """ jpayne@68: ## Removed!! jpayne@68: #def do_illumina_read_level_report(fastq, firstSubsampledFastqFileName, log): jpayne@68: jpayne@68: jpayne@68: """ jpayne@68: STEP23 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: jpayne@68: do_cleanup_readqc jpayne@68: jpayne@68: """ jpayne@68: def do_cleanup_readqc(log): jpayne@68: log.info("\n\n%sSTEP23 - Cleanup <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: jpayne@68: status = "23_cleanup_readqc in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: ## Get option jpayne@68: skipCleanup = RQCReadQcConfig.CFG["skip_cleanup"] jpayne@68: jpayne@68: if not skipCleanup: jpayne@68: retCode = cleanup_readqc(log) jpayne@68: else: jpayne@68: log.warning("File cleaning is skipped.") jpayne@68: retCode = RQCExitCodes.JGI_SUCCESS jpayne@68: jpayne@68: if retCode != RQCExitCodes.JGI_SUCCESS: jpayne@68: log.error("23_cleanup_readqc failed.") jpayne@68: status = "23_cleanup_readqc failed" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: else: jpayne@68: log.info("23_cleanup_readqc complete.") jpayne@68: status = "23_cleanup_readqc complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: jpayne@68: return status jpayne@68: jpayne@68: jpayne@68: def stetch_section_note(table_header): jpayne@68: overview = [ jpayne@68: 'BBTools sketch.sh uses a technique called MinHash to rapidly compare large sequences. ', jpayne@68: 'The result is similar to BLAST, a list of hits from a query sequence to various reference sequences, ', jpayne@68: 'sorted by similarity but the mechanisms are very different. ', jpayne@68: 'For more information, see http://bbtools.jgi.doe.gov' jpayne@68: ] jpayne@68: jpayne@68: legend = { jpayne@68: 'WKID' : 'Weighted Kmer IDentity, which is the kmer identity compensating for differences in size. So, comparing human chr1 to the full human genome would yield 100% WKID but approximately 10% KID.', jpayne@68: 'KID' : 'Kmer IDentity, equal to matches/length; this is the fraction of shared kmers.', jpayne@68: 'ANI' : 'Average Nucleotide Identity, derived from WKID and kmer length.', jpayne@68: 'Complt' : 'Genome completeness (percent of the reference represented in the query). Derived from WKID and KID.', jpayne@68: 'Contam' : 'Contamination (percent of the query that does not match this reference, but matches some other reference).', jpayne@68: 'Depth' : 'Per-kmer depth in sketches, indicating how many times that kmer was seen in the sequence.', jpayne@68: 'Depth2' : 'Repeat-compensated depth', jpayne@68: 'Matches' : 'The number of shared kmers between query and ref.', jpayne@68: 'Unique' : 'The number of shared kmers between query and ref, and no other ref.', jpayne@68: 'Unique2' : '??', jpayne@68: 'Depth2' : 'Repeat-compensated depth', jpayne@68: 'noHit' : 'Number of kmers that did not hit any reference sequence. Though constant for a query, it will be reported differently for different references based on the relative size of the reference and query (if the reference is bigger than the query, it will report all of them).', jpayne@68: 'TaxID' : 'NCBI taxonomic id, when available.', jpayne@68: 'gSize' : 'Estimate of genomic size (number of unique kmers in the genome). This is based on the smallest hash value in the list. This is affected by blacklists or whitelists, and by using an assembly versus raw reads.', jpayne@68: 'gSeqs' : 'Number of sequences used in the sketch.', jpayne@68: 'taxName' : 'NCBI\'s name for that taxID. If there is no taxID, the sequence name will be used.', jpayne@68: } jpayne@68: jpayne@68: html = '' jpayne@68: for name in table_header: jpayne@68: if name != 'taxonomy': jpayne@68: html += html_tag('li', '%s: %s' % (html_tag('span', name, {'class': 'notice'}), html_tag('span', legend.get(name), {'class': 'small-italic'}))) jpayne@68: html = html_tag('span', ''.join(overview), {'class': 'notice'}) + '

Column Legend
' + html_tag('ul', html) jpayne@68: jpayne@68: return html jpayne@68: jpayne@68: jpayne@68: def sketch_table(fname): jpayne@68: # print('DEBUG : %s' % fname) jpayne@68: jpayne@68: html = '' jpayne@68: jpayne@68: if os.path.isfile(fname): jpayne@68: title = None jpayne@68: header = None jpayne@68: data = [] jpayne@68: jpayne@68: with open(fname, 'r') as fh: jpayne@68: for line in fh: jpayne@68: line = line.strip() jpayne@68: if line == '': jpayne@68: continue jpayne@68: if line.startswith('Query:'): jpayne@68: title = line jpayne@68: elif line.startswith('WKID'): jpayne@68: header = line.split('\t') jpayne@68: else: jpayne@68: data.append(line) jpayne@68: if title and header and data: jpayne@68: html += html_th(header) jpayne@68: for line in data: jpayne@68: row = line.split('\t') jpayne@68: html += html_tr(row) jpayne@68: jpayne@68: # html = html_tag('p', title) + html_tag('table', html, attrs={'class': 'data'}) jpayne@68: html = html_tag('table', html, attrs={'class': 'data'}) jpayne@68: return html, header jpayne@68: jpayne@68: def do_html_contam_art_first_n_pb_tr(stats, files, odir, filepath_prefix): jpayne@68: temp = os.path.join(PYDIR, 'template/readqc_artifacts.html') jpayne@68: html = '' jpayne@68: with open(temp, 'r') as fh: jpayne@68: html = fh.read() jpayne@68: jpayne@68: ## the optional artifact contam jpayne@68: artifact_tr = '' jpayne@68: artifact_type = '50' jpayne@68: artifact_val = pipeline_val('illumina read percent contamination artifact 50bp seal', {'type': 'raw'}, stats, files, filepath_prefix) jpayne@68: artifact_file = pipeline_val('artifact_50bp.seal.stats seal', {'type': 'file'}, stats, files, filepath_prefix) jpayne@68: if not artifact_file: jpayne@68: artifact_file = pipeline_val('artifact_20bp.seal.stats seal', {'type': 'file'}, stats, files, filepath_prefix) jpayne@68: artifact_type = '20' jpayne@68: artifact_val = pipeline_val('illumina read percent contamination artifact 20bp seal', {'type': 'raw'}, stats, files, filepath_prefix) jpayne@68: if artifact_file: jpayne@68: html = html.replace('[_CONTAM-ART-FIRST-BP_]', artifact_type) jpayne@68: html = html.replace('[_CONTAM-ART-FIRST-BP-SEAL_]', artifact_file) jpayne@68: html = html.replace('[_CONTAM-ART-FIRST-BP-SEAL-PCT_]', artifact_val) jpayne@68: jpayne@68: return html jpayne@68: jpayne@68: jpayne@68: def do_html_body(odir, filepath_prefix): jpayne@68: print('do_html_body - %s' % filepath_prefix) jpayne@68: temp = os.path.join(PYDIR, 'template/readqc_body_template.html') jpayne@68: statsf = os.path.join(odir, 'readqc_stats.txt') jpayne@68: if not os.path.isfile(statsf): jpayne@68: print('ERROR : file not found: %s' % statsf) jpayne@68: stats = get_dict_obj(statsf) jpayne@68: jpayne@68: filesf = os.path.join(odir, 'readqc_files.txt') jpayne@68: if not os.path.isfile(statsf): jpayne@68: print('ERROR : file not found: %s' % statsf) jpayne@68: files = get_dict_obj(os.path.join(odir, 'readqc_files.txt')) jpayne@68: jpayne@68: ## key (key name in readqc_stats or readqc_files) : {token (space holder in html template), type (value format)} jpayne@68: tok_map = { jpayne@68: ## Average Base Quality section jpayne@68: 'overall bases Q score mean' : {'token' : '[_BASE-QUALITY-SCORE_]', 'type': 'float', 'filter': 1}, jpayne@68: 'overall bases Q score std' : {'token' : '[_BASE-QUALITY-SCORE-STD_]', 'type': 'float', 'filter': 1}, jpayne@68: 'Q30 bases Q score mean' : {'token' : '[_Q30-BASE-QUALITY-SCORE_]', 'type': 'float', 'filter': 1}, jpayne@68: 'Q30 bases Q score std' : {'token' : '[_Q30-BASE-QUALITY-SCORE-STD_]', 'type': 'float', 'filter': 1}, jpayne@68: 'base C30' : {'token' : '[_COUNT-OF-BAESE-Q30_]', 'type': 'bigint'}, jpayne@68: 'base C25' : {'token' : '[_COUNT-OF-BAESE-Q25_]', 'type': 'bigint'}, jpayne@68: 'base C20' : {'token' : '[_COUNT-OF-BAESE-Q20_]', 'type': 'bigint'}, jpayne@68: 'base C15' : {'token' : '[_COUNT-OF-BAESE-Q15_]', 'type': 'bigint'}, jpayne@68: 'base C10' : {'token' : '[_COUNT-OF-BAESE-Q10_]', 'type': 'bigint'}, jpayne@68: 'base C5' : {'token' : '[_COUNT-OF-BAESE-Q5_]', 'type': 'bigint'}, jpayne@68: 'base Q30' : {'token' : '[_PCT-OF-BAESE-Q30_]', 'type': 'raw'}, jpayne@68: 'base Q25' : {'token' : '[_PCT-OF-BAESE-Q25_]', 'type': 'raw'}, jpayne@68: 'base Q20' : {'token' : '[_PCT-OF-BAESE-Q20_]', 'type': 'raw'}, jpayne@68: 'base Q15' : {'token' : '[_PCT-OF-BAESE-Q15_]', 'type': 'raw'}, jpayne@68: 'base Q10' : {'token' : '[_PCT-OF-BAESE-Q10_]', 'type': 'raw'}, jpayne@68: 'base Q5' : {'token' : '[_PCT-OF-BAESE-Q5_]', 'type': 'raw'}, jpayne@68: jpayne@68: ## Average Read Quality section jpayne@68: 'read Q30' : {'token' : '[_PCT-OF-READS-Q30_]', 'type': 'raw'}, jpayne@68: 'read Q25' : {'token' : '[_PCT-OF-READS-Q25_]', 'type': 'raw'}, jpayne@68: 'read Q20' : {'token' : '[_PCT-OF-READS-Q20_]', 'type': 'raw'}, jpayne@68: 'read Q15' : {'token' : '[_PCT-OF-READS-Q15_]', 'type': 'raw'}, jpayne@68: 'read Q10' : {'token' : '[_PCT-OF-READS-Q10_]', 'type': 'raw'}, jpayne@68: 'read Q5' : {'token' : '[_PCT-OF-READS-Q5_]', 'type': 'raw'}, jpayne@68: jpayne@68: 'SUBSAMPLE_RATE' : {'token' : '[_SUBSAMPLE-RATE_]', 'type': 'pct', 'filter': 1}, jpayne@68: 'read base quality stats' : {'token' : '[_AVG-READ-QUAL-HISTO-DATA_]', 'type': 'file', 'filter': 'link', 'label':'data file'}, jpayne@68: 'ILLUMINA_READ_QHIST_D3_HTML_PLOT' : {'token' : '[_AVG-READ-QUAL-HOSTO-D3_]', 'type': 'file', 'filter': 'link', 'label':'interactive plot'}, jpayne@68: 'read qhist plot' : {'token' : '[_AVG-READ-QUALITY-HISTOGRAM_]', 'type': 'file'}, jpayne@68: jpayne@68: ## Average Base Position Quality section jpayne@68: 'read q20 read1' : {'token' : '[_READ_Q20_READ1_]', 'type': 'raw'}, jpayne@68: 'read q20 read2' : {'token' : '[_READ_Q20_READ2_]', 'type': 'raw'}, jpayne@68: 'read qual pos qrpt 1' : {'token' : '[_AVG-BASE-POS-QUAL-HISTO-DATA_]', 'type': 'file', 'filter': 'link', 'label':'data file'}, jpayne@68: 'ILLUMINA_READ_QUAL_POS_PLOT_MERGED_D3_HTML_PLOT' : {'token' : '[_AVG-BASE-POS-QUAL-HISTO-D3_]', 'type': 'file', 'filter': 'link', 'label':'interactive plot'}, jpayne@68: 'read qual pos plot merged' : {'token' : '[_AVG-BASE-POSITION-QUALITY_]', 'type': 'file'}, jpayne@68: jpayne@68: ## Insert Size jpayne@68: 'ILLUMINA_READ_INSERT_SIZE_JOINED_PERC' : {'token' : '[_PCT-READS-JOINED_]', 'type': 'float', 'filter': 1}, jpayne@68: 'ILLUMINA_READ_INSERT_SIZE_AVG_INSERT' : {'token' : '[_PCT-READS-JOINED-AVG_]', 'type': 'raw'}, jpayne@68: 'ILLUMINA_READ_INSERT_SIZE_STD_INSERT' : {'token' : '[_PCT-READS-JOINED-STDDEV_]', 'type': 'raw'}, jpayne@68: 'ILLUMINA_READ_INSERT_SIZE_MODE_INSERT' : {'token' : '[_PCT-READS-JOINED-MODE_]', 'type': 'raw'}, jpayne@68: jpayne@68: 'ILLUMINA_READ_INSERT_SIZE_HISTO_DATA' : {'token' : '[_INSERT-SIZE-HISTO-DATA_]', 'type': 'file', 'filter': 'link', 'label':'data file'}, jpayne@68: 'ILLUMINA_READ_INSERT_SIZE_HISTO_D3_HTML_PLOT' : {'token' : '[_INSERT-SIZE-HISTO-D3_]', 'type': 'file', 'filter': 'link', 'label':'interactive plot'}, jpayne@68: 'ILLUMINA_READ_INSERT_SIZE_HISTO_PLOT' : {'token' : '[_INSERT-SIZE-HISTOGRAM_]', 'type': 'file'}, jpayne@68: jpayne@68: ## Read GC jpayne@68: 'read GC mean' : {'token' : '[_READ-GC-AVG_]', 'type': 'float', 'filter': 1}, jpayne@68: 'read GC std' : {'token' : '[_READ-GC-STDDEV_]', 'type': 'float', 'filter': 1}, jpayne@68: 'read GC text hist' : {'token' : '[_READ-QC-HISTO-DATA_]', 'type': 'file', 'filter': 'link', 'label':'data file'}, jpayne@68: 'ILLUMINA_READ_GC_D3_HTML_PLOT' : {'token' : '[_READ-QC-HISTO-D3_]', 'type': 'file', 'filter': 'link', 'label':'interactive plot'}, jpayne@68: 'read GC plot' : {'token' : '[_READ-GC-HIST_]', 'type': 'file'}, jpayne@68: jpayne@68: ## Cycle Nucleotide Composition jpayne@68: 'read base count text 1' : {'token' : '[_NUC-COMP-FREQ-R1-DATA_]', 'type': 'file', 'filter': 'link', 'label':'data file'}, jpayne@68: 'read base count text 2' : {'token' : '[_NUC-COMP-FREQ-R2-DATA_]', 'type': 'file', 'filter': 'link', 'label':'data file'}, jpayne@68: 'ILLUMINA_READ_BASE_COUNT_D3_HTML_PLOT_1' : {'token' : '[_NUC-COMP-FREQ-R1-D3_]', 'type': 'file', 'filter': 'link', 'label':'interactive plot'}, jpayne@68: 'ILLUMINA_READ_BASE_COUNT_D3_HTML_PLOT_2' : {'token' : '[_NUC-COMP-FREQ-R2-D3_]', 'type': 'file', 'filter': 'link', 'label':'interactive plot'}, jpayne@68: 'read base count plot 1' : {'token' : '[_CYCLE-NUCL-COMPOSITION-READ1_]', 'type': 'file'}, jpayne@68: 'read base count plot 2' : {'token' : '[_CYCLE-NUCL-COMPOSITION-READ2_]', 'type': 'file'}, jpayne@68: jpayne@68: ## Percentage of Common Contaminants jpayne@68: 'illumina read percent contamination artifact seal' : {'token' : '[_CONTAM-ART-SEA-PCT_]', 'type': 'floatstr'}, jpayne@68: 'artifact.seal.stats seal' : {'token' : '[_CONTAM-ART-SEAL_]', 'type': 'file'}, jpayne@68: jpayne@68: 'illumina read percent contamination DNA spikein seal' : {'token' : '[_DNA-SPIKEIN-SEAL_PCT_]', 'type': 'floatstr'}, jpayne@68: 'DNA_spikein.seal.stats seal' : {'token' : '[_DNA-SPIKEIN-SEAL_]', 'type': 'file'}, jpayne@68: jpayne@68: 'illumina read percent contamination RNA spikein seal' : {'token' : '[_RNA-SPIKEIN-SEAL_PCT_]', 'type': 'floatstr'}, jpayne@68: 'RNA_spikein.seal.stats seal' : {'token' : '[_RNA-SPIKEIN-SEAL_]', 'type': 'file'}, jpayne@68: jpayne@68: 'illumina read percent contamination fosmid seal' : {'token' : '[_CONTAM-FOSMID-SEAL-PCT_]', 'type': 'floatstr'}, jpayne@68: 'fosmid.seal.stats seal' : {'token' : '[_CONTAM-FOSMID-SEAL_]', 'type': 'file'}, jpayne@68: jpayne@68: 'illumina read percent contamination fosmid seal' : {'token' : '[_CONTAM-FOSMID-SEAL-PCT_]', 'type': 'floatstr'}, jpayne@68: 'fosmid.seal.stats seal' : {'token' : '[_CONTAM-FOSMID-SEAL_]', 'type': 'file'}, jpayne@68: jpayne@68: 'illumina read percent contamination mitochondrion seal' : {'token' : '[_CONTAM-MITO-SEAL-PCT_]', 'type': 'floatstr'}, jpayne@68: 'mitochondrion.seal.stats seal' : {'token' : '[_CONTAM-MITO-SEAL_]', 'type': 'file'}, jpayne@68: jpayne@68: 'illumina read percent contamination plastid seal' : {'token' : '[_CONTAM-CHLO-SEAL-PCT_]', 'type': 'floatstr'}, jpayne@68: 'plastid.seal.stats seal' : {'token' : '[_CONTAM-CHLO-SEAL_]', 'type': 'file'}, jpayne@68: jpayne@68: 'illumina read percent contamination phix seal' : {'token' : '[_CONTAM-PHIX-SEAL-PCT_]', 'type': 'floatstr'}, jpayne@68: 'phix.seal.stats seal' : {'token' : '[_CONTAM-PHIX-SEAL_]', 'type': 'file'}, jpayne@68: jpayne@68: 'illumina read percent contamination rrna seal' : {'token' : '[_CONTAM-RRNA-SEAL-PCT_]', 'type': 'floatstr'}, jpayne@68: 'rrna.seal.stats seal' : {'token' : '[_CONTAM-RRNA-SEAL_]', 'type': 'file'}, jpayne@68: jpayne@68: 'illumina read percent contamination microbes seal' : {'token' : '[_CONTAM-NON-SYN-SEAL-PCT_]', 'type': 'floatstr'}, jpayne@68: 'microbes.seal.stats seal' : {'token' : '[_CONTAM-NON-SYN-SEAL_]', 'type': 'file'}, jpayne@68: jpayne@68: # 'illumina read percent contamination synthetic seal' : {'token' : '[_CONTAM-SYN-SEAL-PCT_]', 'type': 'floatstr'}, jpayne@68: 'synthetic.seal.stats seal' : {'token' : '[_CONTAM-SYN-SEAL_]', 'type': 'file'}, jpayne@68: jpayne@68: 'illumina read percent contamination adapters seal' : {'token' : '[_CONTAM-ADAPTER-PCT_]', 'type': 'floatstr'}, jpayne@68: 'adapters.seal.stats seal' : {'token' : '[_CONTAM-ADAPTER-SEAL_]', 'type': 'file'}, jpayne@68: jpayne@68: ## Sketch vs NT jpayne@68: 'sketch_vs_nt_output' : {'token' : '[_SKETCH-VS-NT_]', 'type': 'file'}, jpayne@68: # 'sketch_vs_nt_output' : {'token' : '[_SKETCH-VS-NT-BASE_]', 'type': 'file', 'filter': 'base'}, jpayne@68: } jpayne@68: jpayne@68: html = '' jpayne@68: with open(temp, 'r') as fh: jpayne@68: html = fh.read() jpayne@68: jpayne@68: for key in tok_map: jpayne@68: dat = tok_map[key] jpayne@68: val = pipeline_val(key, dat, stats, files, filepath_prefix) jpayne@68: # print('key=%s; %s; ====type=%s' % (key, val, type(val))) jpayne@68: html = html.replace(dat['token'], val) jpayne@68: jpayne@68: artifact_tr = do_html_contam_art_first_n_pb_tr(stats, files, odir, filepath_prefix) jpayne@68: html = html.replace('[_CONTAN-ART-SEAL-FIRST-BP_]', artifact_tr) jpayne@68: jpayne@68: ## synthetic contam jpayne@68: contam_syn_file = pipeline_val('synthetic.seal.stats seal', {'type': 'file', 'filter': 'full'}, stats, files, filepath_prefix) jpayne@68: if contam_syn_file and os.path.isfile(contam_syn_file): jpayne@68: contam_syn_pct_key = 'illumina read percent contamination synthetic seal' jpayne@68: cmd = 'grep "#Matched" %s' % contam_syn_file jpayne@68: stdOut, stdErr, exitCode = run_sh_command(cmd, True) jpayne@68: if exitCode == 0: jpayne@68: toks = stdOut.strip().split() jpayne@68: if len(toks) == 3: jpayne@68: html = html.replace('[_CONTAM-SYN-SEAL-PCT_]', toks[2][:-1]) jpayne@68: jpayne@68: ###--- Sketch jpayne@68: sketch_html = '' jpayne@68: ## add Sketch vs NT if file exists: jpayne@68: sketch_nt_file = pipeline_val('sketch_vs_nt_output', {'type': 'file', 'filter': 'full'}, stats, files, filepath_prefix) jpayne@68: jpayne@68: if os.path.isfile(sketch_nt_file): jpayne@68: basename = os.path.basename(sketch_nt_file) jpayne@68: # href = pipeline_val('sketch_vs_nt_output', {'type': 'file'}, stats, files, filepath_prefix) jpayne@68: sketch_html = html_tag('h4', 'Sketch vs. NT') #+ '
' + 'Raw file:%s' % html_link(href, basename) jpayne@68: sketch_tab, sketch_table_header = sketch_table(sketch_nt_file) jpayne@68: sketch_html += sketch_tab jpayne@68: jpayne@68: jpayne@68: ## add Sketch vs refseq if file exists? jpayne@68: if sketch_html: jpayne@68: sketch_html = html_tag('h2', 'Read Sketch', attrs={'class': 'section-title'}) \ jpayne@68: + stetch_section_note(sketch_table_header) \ jpayne@68: + html_tag('div', sketch_html, {'class': 'section'}) jpayne@68: jpayne@68: html = html.replace('[_SKETCH-TABLE_]', sketch_html) jpayne@68: jpayne@68: # copy the image file jpayne@68: imgdir = os.path.join(PYDIR, 'images') jpayne@68: todir = os.path.join(odir, 'images') jpayne@68: if os.path.isdir(todir): jpayne@68: shutil.rmtree(todir) jpayne@68: shutil.copytree(imgdir, todir, False, None) jpayne@68: jpayne@68: return html jpayne@68: jpayne@68: jpayne@68: def do_html(odir, infile): jpayne@68: # odir = os.path.abspath(odir) ## DO NOT convert!! The relative file path in html need this original odir string matching jpayne@68: jpayne@68: screen('Output dir - %s' % odir) jpayne@68: fname = os.path.basename(infile) jpayne@68: screen('Create HTML page in %s for %s ..' % (odir, fname)) jpayne@68: temp = os.path.join(PYDIR, 'template/template.html') jpayne@68: jpayne@68: with open(temp, 'r') as fh: jpayne@68: html = fh.read() jpayne@68: html = html.replace('[_PAGE-TITLE_]', 'Read QC Report') jpayne@68: html = html.replace('[_REPORT-TITLE_]', 'BBTools Read QC Report') jpayne@68: html = html.replace('[_INPUT-FILE-NAME_]', fname) jpayne@68: html = html.replace('[_REPORT-DATE_]', '{:%Y-%m-%d %H:%M:%S}'.format(datetime.datetime.now())) jpayne@68: jpayne@68: hbody =do_html_body(odir, odir) jpayne@68: html = html.replace('[_REPORT-BODY_]', hbody) jpayne@68: jpayne@68: ## write the html to file jpayne@68: idxfile = os.path.join(odir, 'index.html') jpayne@68: with open(idxfile, 'w') as fh2: jpayne@68: fh2.write(html) jpayne@68: screen('HTML index file written to %s' % idxfile) jpayne@68: jpayne@68: # copy the css file jpayne@68: cssdir = os.path.join(PYDIR, 'css') jpayne@68: todir = os.path.join(odir, 'css') jpayne@68: if os.path.isdir(todir): jpayne@68: shutil.rmtree(todir) jpayne@68: shutil.copytree(cssdir, todir, False, None) jpayne@68: jpayne@68: def screen(txt): jpayne@68: print('.. %s' % txt) jpayne@68: ## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## Main Program jpayne@68: if __name__ == "__main__": jpayne@68: desc = "RQC ReadQC Pipeline" jpayne@68: parser = argparse.ArgumentParser(description=desc) jpayne@68: parser.add_argument("-f", "--fastq", dest="fastq", help="Set input fastq file (full path to fastq)", required=True) jpayne@68: parser.add_argument("-o", "--output-path", dest="outputPath", help="Set output path to write to") jpayne@68: parser.add_argument("-x", "--cut", dest="cutLenOfFirstBases", help="Set Read cut length (bp) for read contamination detection") jpayne@68: parser.add_argument("-v", "--version", action="version", version=VERSION) jpayne@68: jpayne@68: ## For running on Crays jpayne@68: parser.add_argument("-l", "--lib-name", dest="libName", help="Set library name", required=False) jpayne@68: parser.add_argument("-m", "--is-multiplexed", dest="isMultiplexed", help="Set multiplexed data", required=False) jpayne@68: parser.add_argument("-r", "--is-rna", dest="isRna", help="Set RNA data", required=False) jpayne@68: jpayne@68: ## produce html when processing is done jpayne@68: parser.add_argument("-html", "--html", action="store_true", help="Create html file", dest="html", default=False, required=False) jpayne@68: jpayne@68: ## Toggle options jpayne@68: parser.add_argument("-b", "--skip-blast", action="store_true", help="Skip the blast run", dest="skipBlast", default=False, required=False) jpayne@68: jpayne@68: parser.add_argument("-bn", "--skip-blast-nt", action="store_true", help="Skip Blast run against nt", dest="skipBlastNt", default=False, required=False) jpayne@68: parser.add_argument("-br", "--skip-blast-refseq", action="store_true", help="Skip Blast run against refseq.archaea and refseq.bateria", dest="skipBlastRefseq", default=False, required=False) jpayne@68: parser.add_argument("-c", "--skip-cleanup", action="store_true", help="Skip temporary file cleanup", dest="skipCleanup", default=False, required=False) jpayne@68: parser.add_argument("-p", "--pooled-analysis", action="store_true", help="Enable pooled analysis (demultiplexing)", dest="pooledAnalysis", default=False, required=False) jpayne@68: parser.add_argument("-s", "--skip-subsample", action="store_true", help="Skip subsampling.", dest="skipSubsample", default=False, required=False) jpayne@68: parser.add_argument("-z", "--skip-localization", action="store_true", help="Skip database localization", dest="skipDbLocalization", default=False, required=False) jpayne@68: jpayne@68: options = parser.parse_args() jpayne@68: jpayne@68: outputPath = None # output path, defaults to current working directory jpayne@68: fastq = None # full path to input fastq.gz jpayne@68: jpayne@68: status = "start" jpayne@68: nextStepToDo = 0 jpayne@68: jpayne@68: if options.outputPath: jpayne@68: outputPath = options.outputPath jpayne@68: else: jpayne@68: outputPath = os.getcwd() jpayne@68: jpayne@68: if options.fastq: jpayne@68: fastq = options.fastq jpayne@68: jpayne@68: ## create output_directory if it doesn't exist jpayne@68: if not os.path.isdir(outputPath): jpayne@68: os.makedirs(outputPath) jpayne@68: jpayne@68: libName = None ## for illumina_sciclone_analysis() jpayne@68: isRna = None ## for illumina_sciclone_analysis() jpayne@68: isMultiplexed = None ## for illumina_generate_index_sequence_detection_plot() jpayne@68: jpayne@68: bSkipBlast = False jpayne@68: bSkipBlastRefseq = False ## refseq.archaea and refseq.bacteria jpayne@68: bSkipBlastNt = False jpayne@68: bPooledAnalysis = False jpayne@68: jpayne@68: ## RQC-743 Need to specify the first cutting bp for read contam detection (eps. for smRNA ) jpayne@68: firstBptoCut = 50 jpayne@68: if options.cutLenOfFirstBases: jpayne@68: firstBptoCut = options.cutLenOfFirstBases jpayne@68: jpayne@68: if options.libName: jpayne@68: libName = options.libName jpayne@68: jpayne@68: jpayne@68: ## switches jpayne@68: if options.isRna: jpayne@68: isRna = options.isRna jpayne@68: if options.isMultiplexed: jpayne@68: isMultiplexed = options.isMultiplexed jpayne@68: if options.skipBlast: jpayne@68: bSkipBlast = options.skipBlast jpayne@68: if options.skipBlastRefseq: jpayne@68: bSkipBlastRefseq = options.skipBlastRefseq jpayne@68: if options.skipBlastNt: jpayne@68: bSkipBlastNt = options.skipBlastNt jpayne@68: if options.pooledAnalysis: jpayne@68: bPooledAnalysis = options.pooledAnalysis jpayne@68: jpayne@68: jpayne@68: skipSubsampling = options.skipSubsample jpayne@68: jpayne@68: ## Set readqc config jpayne@68: RQCReadQcConfig.CFG["status_file"] = os.path.join(outputPath, "readqc_status.log") jpayne@68: RQCReadQcConfig.CFG["files_file"] = os.path.join(outputPath, "readqc_files.tmp") jpayne@68: RQCReadQcConfig.CFG["stats_file"] = os.path.join(outputPath, "readqc_stats.tmp") jpayne@68: RQCReadQcConfig.CFG["output_path"] = outputPath jpayne@68: RQCReadQcConfig.CFG["skip_cleanup"] = options.skipCleanup jpayne@68: RQCReadQcConfig.CFG["skip_localization"] = options.skipDbLocalization jpayne@68: RQCReadQcConfig.CFG["log_file"] = os.path.join(outputPath, "readqc.log") jpayne@68: jpayne@68: screen("Started readqc pipeline, writing log to: %s" % (RQCReadQcConfig.CFG["log_file"])) jpayne@68: jpayne@68: log = get_logger("readqc", RQCReadQcConfig.CFG["log_file"], LOG_LEVEL, False, True) jpayne@68: log.info("=================================================================") jpayne@68: log.info(" Read Qc Analysis (version %s)", VERSION) jpayne@68: log.info("=================================================================") jpayne@68: jpayne@68: log.info("Starting %s with %s", SCRIPT_NAME, fastq) jpayne@68: jpayne@68: if os.path.isfile(RQCReadQcConfig.CFG["status_file"]): jpayne@68: status = get_status(RQCReadQcConfig.CFG["status_file"], log) jpayne@68: else: jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: jpayne@68: if not os.path.isfile(fastq): jpayne@68: log.error("%s not found, aborting!", fastq) jpayne@68: exit(2) jpayne@68: elif status != "complete": jpayne@68: ## check for fastq file jpayne@68: # if os.path.isfile(fastq): jpayne@68: log.info("Found %s, starting processing.", fastq) jpayne@68: log.info("Latest status = %s", status) jpayne@68: if status != 'start': jpayne@68: nextStepToDo = int(status.split("_")[0]) jpayne@68: if status.find("complete") != -1: jpayne@68: nextStepToDo += 1 jpayne@68: log.info("Next step to do = %s", nextStepToDo) jpayne@68: jpayne@68: if status != 'complete': jpayne@68: bDone = False jpayne@68: cycle = 0 jpayne@68: totalReadNum = 0 jpayne@68: firstSubsampledFastqFileName = "" jpayne@68: secondSubsampledFastqFile = "" jpayne@68: totalReadCount = 0 jpayne@68: subsampledReadNum = 0 jpayne@68: bIsPaired = False jpayne@68: readLength = 0 jpayne@68: firstSubsampledLogFile = os.path.join(outputPath, "subsample", "first_subsampled.txt") jpayne@68: secondSubsampledLogFile = os.path.join(outputPath, "subsample", "second_subsampled.txt") jpayne@68: jpayne@68: while not bDone: jpayne@68: jpayne@68: cycle += 1 jpayne@68: jpayne@68: if bPooledAnalysis: jpayne@68: nextStepToDo = 17 jpayne@68: status = "16_illumina_read_blastn_nt complete" jpayne@68: jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## 1. fast_subsample_fastq_sequences jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: if nextStepToDo == 1 or status == "start": jpayne@68: status, firstSubsampledFastqFileName, totalReadCount, subsampledReadNum, bIsPaired, readLength = do_fast_subsample_fastq_sequences(fastq, skipSubsampling, log) jpayne@68: jpayne@68: if status.endswith("failed"): jpayne@68: log.error("Subsampling failed.") jpayne@68: sys.exit(-1) jpayne@68: jpayne@68: ## if not skipSubsampling and subsampledReadNum == 0: ## too small input file jpayne@68: if not skipSubsampling and subsampledReadNum < RQCReadQc.ILLUMINA_MIN_NUM_READS: ## min=1000 jpayne@68: log.info("Too small input fastq file. Skip subsampling: total number of reads = %s, sampled number of reads = %s", totalReadCount, subsampledReadNum) jpayne@68: skipSubsampling = True jpayne@68: status, firstSubsampledFastqFileName, totalReadCount, subsampledReadNum, bIsPaired, readLength = do_fast_subsample_fastq_sequences(fastq, skipSubsampling, log) jpayne@68: jpayne@68: if status.endswith("failed"): jpayne@68: log.error("Subsampling failed.") jpayne@68: sys.exit(-1) jpayne@68: jpayne@68: subsampledReadNum = totalReadCount jpayne@68: if subsampledReadNum >= RQCReadQc.ILLUMINA_MIN_NUM_READS: jpayne@68: status = "1_illumina_readqc_subsampling complete" jpayne@68: jpayne@68: ## Still too low number of reads -> record ILLUMINA_TOO_SMALL_NUM_READS and quit jpayne@68: ## jpayne@68: if subsampledReadNum == 0 or subsampledReadNum < RQCReadQc.ILLUMINA_MIN_NUM_READS: ## min=1000 jpayne@68: log.info("Too small number of reads (< %s). Stop processing.", RQCReadQc.ILLUMINA_MIN_NUM_READS) jpayne@68: print("WARNING : Too small number of reads (< %s). Stop processing." % RQCReadQc.ILLUMINA_MIN_NUM_READS) jpayne@68: log.info("Completed %s: %s", SCRIPT_NAME, fastq) jpayne@68: jpayne@68: ## Add ILLUMINA_TOO_SMALL_NUM_READS=1 to stats to be used in reporting jpayne@68: statsFile = RQCReadQcConfig.CFG["stats_file"] jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_TOO_SMALL_NUM_READS, "1", log) jpayne@68: jpayne@68: ## move rqc-files.tmp to rqc-files.txt jpayne@68: newFilesFile = os.path.join(outputPath, "readqc_files.txt") jpayne@68: newStatsFile = os.path.join(outputPath, "readqc_stats.txt") jpayne@68: jpayne@68: cmd = "mv %s %s " % (RQCReadQcConfig.CFG["files_file"], newFilesFile) jpayne@68: log.info("mv cmd: %s", cmd) jpayne@68: run_sh_command(cmd, True, log) jpayne@68: jpayne@68: cmd = "mv %s %s " % (RQCReadQcConfig.CFG["stats_file"], newStatsFile) jpayne@68: log.info("mv cmd: %s", cmd) jpayne@68: run_sh_command(cmd, True, log) jpayne@68: jpayne@68: exit(0) jpayne@68: jpayne@68: if status == "1_illumina_readqc_subsampling failed": jpayne@68: bDone = True jpayne@68: jpayne@68: jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## 2. write_unique_20_mers jpayne@68: ## NOTE: fastq = orig input fastq, jpayne@68: ## totalReadCount = total reads in the orig input fastq jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: if nextStepToDo == 2 or status == "1_illumina_readqc_subsampling complete": jpayne@68: ## Cope with restarting jpayne@68: ## save subsamples file in "first_subsampled.txt" so that the file jpayne@68: ## name can be read when started jpayne@68: if not os.path.isfile(firstSubsampledLogFile): jpayne@68: make_dir(os.path.join(outputPath, "subsample")) jpayne@68: with open(firstSubsampledLogFile, "w") as SAMP_FH: jpayne@68: SAMP_FH.write(firstSubsampledFastqFileName + " " + str(totalReadCount) + " " + str(subsampledReadNum) + "\n") jpayne@68: jpayne@68: if firstSubsampledFastqFileName == "" and options.skipSubsample: jpayne@68: if firstSubsampledFastqFileName == "" and os.path.isfile(firstSubsampledLogFile): jpayne@68: with open(firstSubsampledLogFile, "r") as SAMP_FH: jpayne@68: l = SAMP_FH.readline().strip() jpayne@68: t = l.split() jpayne@68: assert len(t) == 2 or len(t) == 3 jpayne@68: firstSubsampledFastqFileName = t[0] jpayne@68: totalReadCount = t[1] jpayne@68: subsampledReadNum = t[2] jpayne@68: log.debug("firstSubsampledFastqFileName=%s, totalReadCount=%s, subsampledReadNum=%s", firstSubsampledFastqFileName, totalReadCount, subsampledReadNum) jpayne@68: jpayne@68: else: jpayne@68: nextStepToDo = 1 jpayne@68: continue jpayne@68: jpayne@68: ## TODO: move getting totalReadNum in the func ??? jpayne@68: ## jpayne@68: log.debug("firstSubsampledFastqFileName=%s, totalReadCount=%s, subsampledReadNum=%s", firstSubsampledFastqFileName, totalReadCount, subsampledReadNum) jpayne@68: jpayne@68: if totalReadCount is None or totalReadCount == 0: jpayne@68: nextStepToDo = 1 jpayne@68: continue jpayne@68: else: jpayne@68: status = do_write_unique_20_mers(fastq, totalReadCount, log) jpayne@68: jpayne@68: if status == "2_unique_mers_sampling failed": jpayne@68: bDone = True jpayne@68: jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## 3. generate read GC histograms: Illumina_read_gc jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: if nextStepToDo == 3 or status == "2_unique_mers_sampling complete": jpayne@68: ## Read the subsampled file name jpayne@68: if firstSubsampledFastqFileName == "" and os.path.isfile(firstSubsampledLogFile): jpayne@68: with open(firstSubsampledLogFile, "r") as SAMP_FH: jpayne@68: firstSubsampledFastqFileName = SAMP_FH.readline().strip().split()[0] jpayne@68: status = do_illumina_read_gc(firstSubsampledFastqFileName, log) jpayne@68: jpayne@68: if status == "3_illumina_read_gc failed": jpayne@68: bDone = True jpayne@68: jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## 4. read_quality_stats jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: if nextStepToDo == 4 or status == "3_illumina_read_gc complete": jpayne@68: if firstSubsampledFastqFileName == "" and os.path.isfile(firstSubsampledLogFile): jpayne@68: with open(firstSubsampledLogFile, "r") as SAMP_FH: jpayne@68: firstSubsampledFastqFileName = SAMP_FH.readline().strip().split()[0] jpayne@68: status = do_read_quality_stats(firstSubsampledFastqFileName, log) jpayne@68: jpayne@68: if status == "4_illumina_read_quality_stats failed": jpayne@68: bDone = True jpayne@68: jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## 5. write_base_quality_stats jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: if nextStepToDo == 5 or status == "4_illumina_read_quality_stats complete": jpayne@68: if firstSubsampledFastqFileName == "" and os.path.isfile(firstSubsampledLogFile): jpayne@68: with open(firstSubsampledLogFile, "r") as SAMP_FH: jpayne@68: firstSubsampledFastqFileName = SAMP_FH.readline().strip().split()[0] jpayne@68: status = do_write_base_quality_stats(firstSubsampledFastqFileName, log) jpayne@68: jpayne@68: if status == "5_illumina_read_quality_stats failed": jpayne@68: bDone = True jpayne@68: jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## 6. illumina_count_q_score jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: if nextStepToDo == 6 or status == "5_illumina_read_quality_stats complete": jpayne@68: if firstSubsampledFastqFileName == "" and os.path.isfile(firstSubsampledLogFile): jpayne@68: with open(firstSubsampledLogFile, "r") as SAMP_FH: jpayne@68: firstSubsampledFastqFileName = SAMP_FH.readline().strip().split()[0] jpayne@68: status = do_illumina_count_q_score(firstSubsampledFastqFileName, log) jpayne@68: jpayne@68: if status == "6_illumina_count_q_score failed": jpayne@68: bDone = True jpayne@68: jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## 7. illumina_calculate_average_quality jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: if nextStepToDo == 7 or status == "6_illumina_count_q_score complete": jpayne@68: ## Let's skip this step. (20140902) jpayne@68: log.info("\n\n%sSTEP7 - Skipping 21mer analysis <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: jpayne@68: status = "7_illumina_calculate_average_quality in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: status = "7_illumina_calculate_average_quality complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## 8. illumina_find_common_motifs jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: if nextStepToDo == 8 or status == "7_illumina_calculate_average_quality complete": jpayne@68: if firstSubsampledFastqFileName == "" and os.path.isfile(firstSubsampledLogFile): jpayne@68: with open(firstSubsampledLogFile, "r") as SAMP_FH: jpayne@68: firstSubsampledFastqFileName = SAMP_FH.readline().strip().split()[0] jpayne@68: status = do_illumina_find_common_motifs(firstSubsampledFastqFileName, log) jpayne@68: jpayne@68: if status == "8_illumina_find_common_motifs failed": jpayne@68: bDone = True jpayne@68: jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## 10. illumina_run_tagdust jpayne@68: ## jpayne@68: ## NOTE: This step will be skipped. No need to run. jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: # if nextStepToDo == 10 or status == "9_illumina_run_dedupe complete": jpayne@68: if nextStepToDo == 9 or status == "8_illumina_find_common_motifs complete": jpayne@68: ## 20131023 skip this step jpayne@68: log.info("\n\n%sSTEP10 - Skipping tag dust <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: status = "10_illumina_run_tagdust in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: status = "10_illumina_run_tagdust complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## 11. illumina_detect_read_contam jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: if nextStepToDo == 11 or status == "10_illumina_run_tagdust complete": jpayne@68: if firstSubsampledFastqFileName == "" and os.path.isfile(firstSubsampledLogFile): jpayne@68: with open(firstSubsampledLogFile, "r") as SAMP_FH: jpayne@68: firstSubsampledFastqFileName = SAMP_FH.readline().strip().split()[0] jpayne@68: status = do_illumina_detect_read_contam(firstSubsampledFastqFileName, firstBptoCut, log) jpayne@68: jpayne@68: if status == "11_illumina_detect_read_contam failed": jpayne@68: bDone = True jpayne@68: jpayne@68: jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## 13. illumina_subsampling_read_megablast jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: if nextStepToDo == 12 or status == "11_illumina_detect_read_contam complete": jpayne@68: if not bSkipBlast: jpayne@68: if firstSubsampledFastqFileName == "" and os.path.isfile(firstSubsampledLogFile): jpayne@68: with open(firstSubsampledLogFile, "r") as SAMP_FH: jpayne@68: l = SAMP_FH.readline().strip() jpayne@68: firstSubsampledFastqFileName = l.split()[0] jpayne@68: totalReadCount = int(l.split()[1]) jpayne@68: subsampledReadNum = int(l.split()[2]) jpayne@68: jpayne@68: status, secondSubsampledFastqFile, second_read_cnt_total, second_read_cnt_sampled = do_illumina_subsampling_read_blastn(firstSubsampledFastqFileName, skipSubsampling, subsampledReadNum, log) jpayne@68: log.debug("status=%s, secondSubsampledFastqFile=%s, second_read_cnt_total=%s, second_read_cnt_sampled=%s", status, secondSubsampledFastqFile, second_read_cnt_total, second_read_cnt_sampled) jpayne@68: jpayne@68: else: jpayne@68: statsFile = RQCReadQcConfig.CFG["stats_file"] jpayne@68: append_rqc_stats(statsFile, ReadqcStats.ILLUMINA_SKIP_BLAST, "1", log) jpayne@68: log.info("\n\n%sSTEP13 - Run subsampling for Blast search <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: log.info("13_illumina_subsampling_read_megablast skipped.\n") jpayne@68: status = "13_illumina_subsampling_read_megablast complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: if status == "13_illumina_subsampling_read_megablast failed": jpayne@68: bDone = True jpayne@68: jpayne@68: #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: # 14. illumina_read_blastn_refseq_archaea jpayne@68: #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: if nextStepToDo == 14 or status == "13_illumina_subsampling_read_megablast complete": jpayne@68: if not bSkipBlast and not bSkipBlastRefseq: jpayne@68: if secondSubsampledFastqFile == "" and not os.path.isfile(secondSubsampledLogFile): jpayne@68: nextStepToDo = 13 jpayne@68: continue jpayne@68: jpayne@68: if secondSubsampledFastqFile != "" and not os.path.isfile(secondSubsampledLogFile): jpayne@68: make_dir(os.path.join(outputPath, "subsample")) jpayne@68: with open(secondSubsampledLogFile, "w") as SAMP_FH: jpayne@68: SAMP_FH.write(secondSubsampledFastqFile + " " + str(second_read_cnt_total) + " " + str(second_read_cnt_sampled) + "\n") jpayne@68: else: jpayne@68: with open(secondSubsampledLogFile, "r") as SAMP_FH: jpayne@68: l = SAMP_FH.readline().strip() jpayne@68: secondSubsampledFastqFile = l.split()[0] jpayne@68: second_read_cnt_total = int(l.split()[1]) jpayne@68: second_read_cnt_sampled = int(l.split()[2]) jpayne@68: jpayne@68: status = do_illumina_read_blastn_refseq_archaea(secondSubsampledFastqFile, second_read_cnt_sampled, log) jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: else: jpayne@68: log.info("\n\n%sSTEP14 - Run illumina_read_blastn_refseq_archaea analysis <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: log.info("14_illumina_read_blastn_refseq_archaea skipped.\n") jpayne@68: status = "14_illumina_read_blastn_refseq_archaea in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: status = "14_illumina_read_blastn_refseq_archaea complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: if status == "14_illumina_read_blastn_refseq_archaea failed": jpayne@68: bDone = True jpayne@68: jpayne@68: #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: # 15. illumina_read_blastn_refseq_bacteria jpayne@68: #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: if nextStepToDo == 15 or status == "14_illumina_read_blastn_refseq_archaea complete": jpayne@68: if not bSkipBlast and not bSkipBlastRefseq: jpayne@68: if secondSubsampledFastqFile == "" and not os.path.isfile(secondSubsampledLogFile): jpayne@68: nextStepToDo = 13 jpayne@68: continue jpayne@68: jpayne@68: if secondSubsampledFastqFile != "" and not os.path.isfile(secondSubsampledLogFile): jpayne@68: make_dir(os.path.join(outputPath, "subsample")) jpayne@68: with open(secondSubsampledLogFile, "w") as SAMP_FH: jpayne@68: SAMP_FH.write(secondSubsampledFastqFile + " " + str(second_read_cnt_total) + " " + str(second_read_cnt_sampled) + "\n") jpayne@68: else: jpayne@68: with open(secondSubsampledLogFile, "r") as SAMP_FH: jpayne@68: l = SAMP_FH.readline().strip() jpayne@68: secondSubsampledFastqFile = l.split()[0] jpayne@68: second_read_cnt_total = int(l.split()[1]) jpayne@68: second_read_cnt_sampled = int(l.split()[2]) jpayne@68: jpayne@68: status = do_illumina_read_blastn_refseq_bacteria(secondSubsampledFastqFile, log) jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: else: jpayne@68: log.info("\n\n%sSTEP15 - Run illumina_read_blastn_refseq_bacteria analysis <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: log.info("15_illumina_read_blastn_refseq_bacteria skipped.\n") jpayne@68: status = "15_illumina_read_blastn_refseq_bacteria in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: status = "15_illumina_read_blastn_refseq_bacteria complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: if status == "15_illumina_read_blastn_refseq_bacteria failed": jpayne@68: bDone = True jpayne@68: jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## 16. illumina_read_blastn_nt jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: if nextStepToDo == 16 or status == "15_illumina_read_blastn_refseq_bacteria complete": jpayne@68: if not bSkipBlast and not bSkipBlastNt: jpayne@68: if secondSubsampledFastqFile == "" and not os.path.isfile(secondSubsampledLogFile): jpayne@68: nextStepToDo = 13 jpayne@68: continue jpayne@68: jpayne@68: if secondSubsampledFastqFile == "" and os.path.isfile(secondSubsampledLogFile): jpayne@68: with open(secondSubsampledLogFile, "r") as SAMP_FH: jpayne@68: l = SAMP_FH.readline().strip() jpayne@68: secondSubsampledFastqFile = l.split()[0] jpayne@68: second_read_cnt_total = int(l.split()[1]) jpayne@68: second_read_cnt_sampled = int(l.split()[2]) jpayne@68: jpayne@68: status = do_illumina_read_blastn_nt(secondSubsampledFastqFile, log) jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: else: jpayne@68: log.info("\n\n%sSTEP16 - Run illumina_read_blastn_nt analysis <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: log.info("16_illumina_read_blastn_nt skipped.\n") jpayne@68: status = "16_illumina_read_blastn_nt in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: status = "16_illumina_read_blastn_nt complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: if status == "16_illumina_read_blastn_nt failed": jpayne@68: bDone = True jpayne@68: jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## 17. multiplex_statistics jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: if nextStepToDo == 17 or status == "16_illumina_read_blastn_nt complete": jpayne@68: status = do_illumina_multiplex_statistics(fastq, log, isMultiplexed=isMultiplexed) jpayne@68: jpayne@68: if status == "17_multiplex_statistics failed": jpayne@68: bDone = True jpayne@68: jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## 18. end_of_read_illumina_adapter_check jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: if nextStepToDo == 18 or status == "17_multiplex_statistics complete": jpayne@68: if firstSubsampledFastqFileName == "" and os.path.isfile(firstSubsampledLogFile): jpayne@68: with open(firstSubsampledLogFile, "r") as SAMP_FH: jpayne@68: firstSubsampledFastqFileName = SAMP_FH.readline().strip().split()[0] jpayne@68: status = do_end_of_read_illumina_adapter_check(firstSubsampledFastqFileName, log) jpayne@68: jpayne@68: if status == "18_end_of_read_illumina_adapter_check failed": jpayne@68: bDone = True jpayne@68: jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## 19. insert size analysis (bbmerge.sh) jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: if nextStepToDo == 19 or status == "18_end_of_read_illumina_adapter_check complete": jpayne@68: status = do_insert_size_analysis(fastq, log) jpayne@68: jpayne@68: if status == "19_insert_size_analysis failed": jpayne@68: bDone = True jpayne@68: jpayne@68: jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## 21. sketch vs nt, refseq, silva jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: # if nextStepToDo == 21 or status == "20_gc_divergence_analysis complete": jpayne@68: if nextStepToDo == 20 or status == "19_insert_size_analysis complete": jpayne@68: status = do_sketch_vs_nt_refseq_silva(fastq, log) jpayne@68: jpayne@68: if status == "21_sketch_vs_nt_refseq_silva failed": jpayne@68: bDone = True jpayne@68: jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## 22. postprocessing & reporting jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: if nextStepToDo == 22 or status == "21_sketch_vs_nt_refseq_silva complete": jpayne@68: ## 20131023 skip this step jpayne@68: log.info("\n\n%sSTEP22 - Run illumina_readqc_report_postprocess: mv rqc-*.tmp to rqc-*.txt <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<%s\n", color['pink'], color['']) jpayne@68: status = "22_illumina_readqc_report_postprocess in progress" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: ## move rqc-files.tmp to rqc-files.txt jpayne@68: newFilesFile = os.path.join(outputPath, "readqc_files.txt") jpayne@68: newStatsFile = os.path.join(outputPath, "readqc_stats.txt") jpayne@68: jpayne@68: cmd = "mv %s %s " % (RQCReadQcConfig.CFG["files_file"], newFilesFile) jpayne@68: log.info("mv cmd: %s", cmd) jpayne@68: run_sh_command(cmd, True, log) jpayne@68: jpayne@68: cmd = "mv %s %s " % (RQCReadQcConfig.CFG["stats_file"], newStatsFile) jpayne@68: log.info("mv cmd: %s", cmd) jpayne@68: run_sh_command(cmd, True, log) jpayne@68: jpayne@68: log.info("22_illumina_readqc_report_postprocess complete.") jpayne@68: status = "22_illumina_readqc_report_postprocess complete" jpayne@68: checkpoint_step_wrapper(status) jpayne@68: jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: ## 23. Cleanup jpayne@68: ##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ jpayne@68: if nextStepToDo == 23 or status == "22_illumina_readqc_report_postprocess complete": jpayne@68: status = do_cleanup_readqc(log) jpayne@68: jpayne@68: if status == "23_cleanup_readqc failed": jpayne@68: bDone = True jpayne@68: jpayne@68: if status == "23_cleanup_readqc complete": jpayne@68: status = "complete" ## FINAL COMPLETE! jpayne@68: bDone = True jpayne@68: jpayne@68: ## don't cycle more than 10 times ... jpayne@68: if cycle > 10: jpayne@68: bDone = True jpayne@68: jpayne@68: if status != "complete": jpayne@68: log.info("Status %s", status) jpayne@68: else: jpayne@68: log.info("\n\nCompleted %s: %s", SCRIPT_NAME, fastq) jpayne@68: checkpoint_step_wrapper("complete") jpayne@68: log.info("Pipeline processing Done.") jpayne@68: print("Pipeline processing Done.") jpayne@68: jpayne@68: else: jpayne@68: log.info("Pipeline processing is already complete, skip.") jpayne@68: print("Pipeline processing is already complete, skip.") jpayne@68: jpayne@68: if status == 'complete' and options.html: jpayne@68: do_html(outputPath, fastq) jpayne@68: jpayne@68: jpayne@68: exit(0) jpayne@68: jpayne@68: jpayne@68: ## EOF