<commit_before>"""Quality control and summary metrics for next-gen alignments and analysis.
"""
import collections
import contextlib
import csv
import os
import glob
import shutil
import subprocess

import pandas as pd
import lxml.html
import yaml
from datetime import datetime
from collections import defaultdict

try:
    from fadapa import Fadapa
except ImportError:
    Fadapa = None
import pybedtools
import pysam
import toolz as tz
import toolz.dicttoolz as dtz

from bcbio import bam, utils
from bcbio.distributed.transaction import file_transaction, tx_tmpdir
from bcbio.log import logger
from bcbio.pipeline import config_utils, run_info
from bcbio.install import _get_data_dir
from bcbio.provenance import do
import bcbio.pipeline.datadict as dd
from bcbio.variation import bedutils
from bcbio.variation import coverage as cov
from bcbio.ngsalign.postalign import dedup_bam
from bcbio.rnaseq import gtf

# ## High level functions to generate summary

def generate_parallel(samples, run_parallel):
    """Provide parallel preparation of summary information for alignment and variant calling.
    """
    samples = run_parallel("pipeline_summary", samples)
    if any(x[0].get("analysis", "").lower().startswith(("standard", "variant")) for x in samples):
        samples = run_parallel("coverage_report", samples)
    qsign_info = run_parallel("qsignature_summary", [samples])
    samples = run_parallel("qc_report_summary", [samples])
    samples = run_parallel("multiqc_summary", [samples])
    summary_file = write_project_summary(samples, qsign_info)
    out = []
    for data in samples:
        if "summary" not in data[0]:
            data[0]["summary"] = {}
        data[0]["summary"]["project"] = summary_file
        if qsign_info:
            data[0]["summary"]["mixup_check"] = qsign_info[0]["out_dir"]
        out.append(data)
    out = _add_researcher_summary(out, summary_file)
    return out

def pipeline_summary(data):
    """Provide summary information on processing sample.
    """
    data = utils.to_single_data(data)
    work_bam = data.get("align_bam")
    if data["analysis"].lower().startswith("smallrna-seq"):
        work_bam = data["clean_fastq"]
        data["summary"] = _run_qc_tools(work_bam, data)
    elif data["analysis"].lower().startswith("chip-seq"):
        work_bam = data["raw_bam"]
        data["summary"] = _run_qc_tools(work_bam, data)
    elif dd.get_ref_file(data) is not None and work_bam and work_bam.endswith(".bam"):
        logger.info("Generating summary files: %s" % dd.get_sample_name(data))
        data["summary"] = _run_qc_tools(work_bam, data)
    return [[data]]

def prep_pdf(qc_dir, config):
    """Create PDF from HTML summary outputs in QC directory.

    Requires wkhtmltopdf installed: http://www.msweet.org/projects.php?Z1
    Thanks to: https://www.biostars.org/p/16991/

    Works around issues with CSS conversion on CentOS by adjusting CSS.
    """
    html_file = os.path.join(qc_dir, "fastqc", "fastqc_report.html")
    html_fixed = "%s-fixed%s" % os.path.splitext(html_file)
    try:
        topdf = config_utils.get_program("wkhtmltopdf", config)
    except config_utils.CmdNotFound:
        topdf = None
    if topdf and utils.file_exists(html_file):
        out_file = "%s.pdf" % os.path.splitext(html_file)[0]
        if not utils.file_exists(out_file):
            cmd = ("sed 's/div.summary/div.summary-no/' %s | sed 's/div.main/div.main-no/' > %s"
                   % (html_file, html_fixed))
            do.run(cmd, "Fix fastqc CSS to be compatible with wkhtmltopdf")
            cmd = [topdf, html_fixed, out_file]
            do.run(cmd, "Convert QC HTML to PDF")
        return out_file

def _run_qc_tools(bam_file, data):
    """Run a set of third party quality control tools, returning QC directory and metrics.

        :param bam_file: alignments in bam format
        :param data: dict with all configuration information

        :returns: dict with output of different tools
    """
    metrics = {}
    to_run = []
    if "fastqc" not in tz.get_in(("config", "algorithm", "tools_off"), data, []):
        to_run.append(("fastqc", _run_fastqc))
    if data["analysis"].lower().startswith("rna-seq"):
        to_run += [("samtools", _run_samtools_stats)]
        if gtf.is_qualimap_compatible(dd.get_gtf_file(data)):
            to_run.append(("qualimap", _rnaseq_qualimap))
    elif data["analysis"].lower().startswith("chip-seq"):
        bam_file = bam_file.replace(".unique", "")
        to_run += [("samtools", _run_samtools_stats)]
    elif not data["analysis"].lower().startswith("smallrna-seq"):
        to_run += [("samtools", _run_samtools_stats), ("gemini", _run_gemini_stats)]

    if data["analysis"].lower().startswith(("standard", "variant", "variant2")):
        to_run.append(["qsignature", _run_qsignature_generator])
        if any([tool in tz.get_in(("config", "algorithm", "tools_on"), data, []) for tool in ["qualimap", "qualimap_full"]]):
            to_run.append(("qualimap", _run_qualimap))
    qc_dir = utils.safe_makedir(os.path.join(data["dirs"]["work"], "qc", data["description"]))
    metrics = {}
    qc_out = {}
    for program_name, qc_fn in to_run:
        cur_qc_dir = os.path.join(qc_dir, program_name)
        out = qc_fn(bam_file, data, cur_qc_dir)
        qc_files = None
        if out and isinstance(out, dict):
            metrics.update(out)
        elif out and isinstance(out, basestring) and os.path.exists(out):
            qc_files = {"base": out, "secondary": []}
        if not qc_files:
            qc_files = _organize_qc_files(program_name, cur_qc_dir)
        if qc_files:
            qc_out[program_name] = qc_files
    if data['config']["algorithm"].get("kraken", None):
        if data["analysis"].lower().startswith("smallrna-seq"):
            logger.info("Kraken is not compatible with srnaseq pipeline yet.")
        else:
            ratio = bam.get_aligned_reads(bam_file, data)
            cur_metrics = _run_kraken(data, ratio)
            metrics.update(cur_metrics)

    bam.remove("%s-downsample%s" % os.path.splitext(bam_file))

    metrics["Name"] = dd.get_sample_name(data)
    metrics["Quality format"] = dd.get_quality_format(data).lower()
    return {"qc": qc_out, "metrics": metrics}

def _organize_qc_files(program, qc_dir):
    """Organize outputs from quality control runs into a base file and secondary outputs.

    Provides compatibility with CWL output. Returns None if no files created during processing.
    """
    base_files = {"fastqc": "fastqc_report.html", "qualimap": "qualimapReport.html"}
    if os.path.exists(qc_dir):
        out_files = []
        for fname in [os.path.join(qc_dir, x) for x in os.listdir(qc_dir)]:
            if os.path.isfile(fname):
                out_files.append(fname)
            elif os.path.isdir(fname) and not fname == "tx":
                out_files.extend([os.path.join(fname, x) for x in os.listdir(fname)])
        if len(out_files) > 0:
            if len(out_files) == 1:
                base = out_files[0]
                secondary = []
            else:
                base = None
                if program in base_files:
                    base_choices = [x for x in out_files if x.endswith("/%s" % base_files[program])]
                    if len(base_choices) == 1:
                        base = base_choices[0]
                if not base:
                    base = out_files[0]
                secondary = [x for x in out_files if x != base]
            return {"base": base, "secondary": secondary}

# ## Generate project level QC summary for quickly assessing large projects

def write_project_summary(samples, qsign_info=None):
    """Write project summary information on the provided samples.
    write out dirs, genome resources,

    """
    work_dir = samples[0][0]["dirs"]["work"]
    out_file = os.path.join(work_dir, "project-summary.yaml")
    upload_dir = (os.path.join(work_dir, samples[0][0]["upload"]["dir"])<commit_msg>QC: avoid passing along transactional directories<commit_after>"""Quality control and summary metrics for next-gen alignments and analysis.
"""
import collections
import contextlib
import csv
import os
import glob
import shutil
import subprocess

import pandas as pd
import lxml.html
import yaml
from datetime import datetime
from collections import defaultdict

try:
    from fadapa import Fadapa
except ImportError:
    Fadapa = None
import pybedtools
import pysam
import toolz as tz
import toolz.dicttoolz as dtz

from bcbio import bam, utils
from bcbio.distributed.transaction import file_transaction, tx_tmpdir
from bcbio.log import logger
from bcbio.pipeline import config_utils, run_info
from bcbio.install import _get_data_dir
from bcbio.provenance import do
import bcbio.pipeline.datadict as dd
from bcbio.variation import bedutils
from bcbio.variation import coverage as cov
from bcbio.ngsalign.postalign import dedup_bam
from bcbio.rnaseq import gtf

# ## High level functions to generate summary

def generate_parallel(samples, run_parallel):
    """Provide parallel preparation of summary information for alignment and variant calling.
    """
    samples = run_parallel("pipeline_summary", samples)
    if any(x[0].get("analysis", "").lower().startswith(("standard", "variant")) for x in samples):
        samples = run_parallel("coverage_report", samples)
    qsign_info = run_parallel("qsignature_summary", [samples])
    samples = run_parallel("qc_report_summary", [samples])
    samples = run_parallel("multiqc_summary", [samples])
    summary_file = write_project_summary(samples, qsign_info)
    out = []
    for data in samples:
        if "summary" not in data[0]:
            data[0]["summary"] = {}
        data[0]["summary"]["project"] = summary_file
        if qsign_info:
            data[0]["summary"]["mixup_check"] = qsign_info[0]["out_dir"]
        out.append(data)
    out = _add_researcher_summary(out, summary_file)
    return out

def pipeline_summary(data):
    """Provide summary information on processing sample.
    """
    data = utils.to_single_data(data)
    work_bam = data.get("align_bam")
    if data["analysis"].lower().startswith("smallrna-seq"):
        work_bam = data["clean_fastq"]
        data["summary"] = _run_qc_tools(work_bam, data)
    elif data["analysis"].lower().startswith("chip-seq"):
        work_bam = data["raw_bam"]
        data["summary"] = _run_qc_tools(work_bam, data)
    elif dd.get_ref_file(data) is not None and work_bam and work_bam.endswith(".bam"):
        logger.info("Generating summary files: %s" % dd.get_sample_name(data))
        data["summary"] = _run_qc_tools(work_bam, data)
    return [[data]]

def prep_pdf(qc_dir, config):
    """Create PDF from HTML summary outputs in QC directory.

    Requires wkhtmltopdf installed: http://www.msweet.org/projects.php?Z1
    Thanks to: https://www.biostars.org/p/16991/

    Works around issues with CSS conversion on CentOS by adjusting CSS.
    """
    html_file = os.path.join(qc_dir, "fastqc", "fastqc_report.html")
    html_fixed = "%s-fixed%s" % os.path.splitext(html_file)
    try:
        topdf = config_utils.get_program("wkhtmltopdf", config)
    except config_utils.CmdNotFound:
        topdf = None
    if topdf and utils.file_exists(html_file):
        out_file = "%s.pdf" % os.path.splitext(html_file)[0]
        if not utils.file_exists(out_file):
            cmd = ("sed 's/div.summary/div.summary-no/' %s | sed 's/div.main/div.main-no/' > %s"
                   % (html_file, html_fixed))
            do.run(cmd, "Fix fastqc CSS to be compatible with wkhtmltopdf")
            cmd = [topdf, html_fixed, out_file]
            do.run(cmd, "Convert QC HTML to PDF")
        return out_file

def _run_qc_tools(bam_file, data):
    """Run a set of third party quality control tools, returning QC directory and metrics.

        :param bam_file: alignments in bam format
        :param data: dict with all configuration information

        :returns: dict with output of different tools
    """
    metrics = {}
    to_run = []
    if "fastqc" not in tz.get_in(("config", "algorithm", "tools_off"), data, []):
        to_run.append(("fastqc", _run_fastqc))
    if data["analysis"].lower().startswith("rna-seq"):
        to_run += [("samtools", _run_samtools_stats)]
        if gtf.is_qualimap_compatible(dd.get_gtf_file(data)):
            to_run.append(("qualimap", _rnaseq_qualimap))
    elif data["analysis"].lower().startswith("chip-seq"):
        bam_file = bam_file.replace(".unique", "")
        to_run += [("samtools", _run_samtools_stats)]
    elif not data["analysis"].lower().startswith("smallrna-seq"):
        to_run += [("samtools", _run_samtools_stats), ("gemini", _run_gemini_stats)]

    if data["analysis"].lower().startswith(("standard", "variant", "variant2")):
        to_run.append(["qsignature", _run_qsignature_generator])
        if any([tool in tz.get_in(("config", "algorithm", "tools_on"), data, []) for tool in ["qualimap", "qualimap_full"]]):
            to_run.append(("qualimap", _run_qualimap))
    qc_dir = utils.safe_makedir(os.path.join(data["dirs"]["work"], "qc", data["description"]))
    metrics = {}
    qc_out = {}
    for program_name, qc_fn in to_run:
        cur_qc_dir = os.path.join(qc_dir, program_name)
        out = qc_fn(bam_file, data, cur_qc_dir)
        qc_files = None
        if out and isinstance(out, dict):
            metrics.update(out)
        elif out and isinstance(out, basestring) and os.path.exists(out):
            qc_files = {"base": out, "secondary": []}
        if not qc_files:
            qc_files = _organize_qc_files(program_name, cur_qc_dir)
        if qc_files:
            qc_out[program_name] = qc_files
    if data['config']["algorithm"].get("kraken", None):
        if data["analysis"].lower().startswith("smallrna-seq"):
            logger.info("Kraken is not compatible with srnaseq pipeline yet.")
        else:
            ratio = bam.get_aligned_reads(bam_file, data)
            cur_metrics = _run_kraken(data, ratio)
            metrics.update(cur_metrics)

    bam.remove("%s-downsample%s" % os.path.splitext(bam_file))

    metrics["Name"] = dd.get_sample_name(data)
    metrics["Quality format"] = dd.get_quality_format(data).lower()
    return {"qc": qc_out, "metrics": metrics}

def _organize_qc_files(program, qc_dir):
    """Organize outputs from quality control runs into a base file and secondary outputs.

    Provides compatibility with CWL output. Returns None if no files created during processing.
    """
    base_files = {"fastqc": "fastqc_report.html", "qualimap": "qualimapReport.html"}
    if os.path.exists(qc_dir):
        out_files = []
        for fname in [os.path.join(qc_dir, x) for x in os.listdir(qc_dir)]:
            if os.path.isfile(fname):
                out_files.append(fname)
            elif os.path.isdir(fname) and not fname.endswith("tx"):
                out_files.extend([os.path.join(fname, x) for x in os.listdir(fname)])
        if len(out_files) > 0:
            if len(out_files) == 1:
                base = out_files[0]
                secondary = []
            else:
                base = None
                if program in base_files:
                    base_choices = [x for x in out_files if x.endswith("/%s" % base_files[program])]
                    if len(base_choices) == 1:
                        base = base_choices[0]
                if not base:
                    base = out_files[0]
                secondary = [x for x in out_files if x != base]
            return {"base": base, "secondary": secondary}

# ## Generate project level QC summary for quickly assessing large projects

def write_project_summary(samples, qsign_info=None):
    """Write project summary information on the provided samples.
    write out dirs, genome resources,

    """
    work_dir = samples[0][0]["dirs"]["work"]
    out_file = os.path.join(work_dir, "project-summary.yaml")
    upload_dir = (os.path.join(work_dir, samples[0][0]["upload"]["dir"])
{'fraction_non_alphanumeric': 0.12310191839030182, 'fraction_numerical': 0.0029994376054489785, 'mean_word_length': 3.145040145040145, 'pattern_counts': {'":': 17, '<': 3, '<?xml version=': 0, '>': 7, 'https://': 2, 'lorem ipsum': 0, 'www.': 4, 'xml': 2}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4909026', 'n_tokens_mistral': 5199, 'n_tokens_neox': 5070, 'n_words': 1255}