Determine if the pipeline can proceed with a BAM file, or needs fastq conversion.
340860
def _pipeline_needs_fastq(config, data):
    
    aligner = config["algorithm"].get("aligner")
    support_bam = aligner in alignment.metadata.get("support_bam", [])
    return aligner and not support_bam
239016
def has_fast_mps_support():
    
    lib = _load_tcmps_lib()
    if lib is None:
        return False

    c_bool = _ctypes.c_bool()
    ret = lib.TCMPSHasHighPowerMetalDevice(_ctypes.byref(c_bool))
    return ret == 0 and c_bool.value
307631
def validate_fastq(self):
        
        logging.info('Validating FASTQ files')
        with progressbar(self.metadata) as bar:
            for sample in bar:
                # Tiny files can pass the validation tests - ensure that they don't
                size = os.path.getsize(sample.general.fastqfiles[0])
                if size >= 1000000:
                    # Try to run reformat.sh on the reads - on any errors try to run repair.sh
                    try:
                        out, err, cmd = bbtools.validate_reads(forward_in=sample.general.fastqfiles[0],
                                                               returncmd=True)
                        write_to_logfile(out, err, self.logfile, sample.general.logout, sample.general.logerr,
                                         None, None)
                    except CalledProcessError:
                        # Set the file names for the reformatted and repaired files
                        outputfile1 = os.path.join(sample.general.outputdirectory, '{}_reformatted_R1.fastq.gz'
                                                   .format(sample.name))
                        repair_file1 = os.path.join(sample.general.outputdirectory, '{}_repaired_R1.fastq.gz'
                                                    .format(sample.name))
                        if len(sample.general.fastqfiles) == 2:
                            outputfile2 = os.path.join(sample.general.outputdirectory, '{}_reformatted_R2.fastq.gz'
                                                       .format(sample.name))
                            repair_file2 = os.path.join(sample.general.outputdirectory, '{}_repaired_R2.fastq.gz'
                                                        .format(sample.name))
                        else:
                            outputfile2 = str()
                            repair_file2 = str()
                        # Use reformat.sh to repair the reads. If this fails, discard the sample from the analyses
                        try:
                            logging.warning('Errors detected in FASTQ files for sample {sample}. '
                                            'Please check the following files for details {log} {logout} {logerr}. '
                                            'The pipeline will use reformat.sh to attempt to repair issues'
                                            .format(sample=sample.name,
                                                    log=self.logfile,
                                                    logout=sample.general.logout,
                                                    logerr=sample.general.logerr))
                            if not os.path.isfile(outputfile1):
                                # Run reformat.sh
                                out, err, cmd = bbtools.reformat_reads(forward_in=sample.general.fastqfiles[0],
                                                                       forward_out=outputfile1,
                                                                       returncmd=True)
                                write_to_logfile(out, err, self.logfile, sample.general.logout, sample.general.logerr,
                                                 None, None)
                                # Run repair.sh (if necessary)
                                if outputfile2:
                                    out, err, cmd = bbtools.repair_reads(forward_in=outputfile1,
                                                                         forward_out=repair_file1,
                                                                         returncmd=True)
                                    write_to_logfile(out, err, self.logfile, sample.general.logout,
                                                     sample.general.logerr, None, None)
                            # Ensure that the output file(s) exist before declaring this a success
                            if os.path.isfile(outputfile1):
                                # Update the fastqfiles attribute to point to the repaired files
                                sample.general.fastqfiles = [repair_file1, repair_file2] if repair_file2 \
                                    else [outputfile1]
                        except CalledProcessError:
                            # The file(s) can be created even if there is STDERR from reformat.sh
                            if os.path.isfile(outputfile1) and outputfile2:
                                try:
                                    out, err, cmd = bbtools.repair_reads(forward_in=outputfile1,
                                                                         forward_out=repair_file1,
                                                                         returncmd=True)
                                    write_to_logfile(out, err, self.logfile, sample.general.logout,
                                                     sample.general.logerr, None, None)
                                    # Update the fastqfiles attribute to point to the repaired files
                                    sample.general.fastqfiles = [repair_file1, repair_file2] if repair_file2 else \
                                        [repair_file1]
                                except CalledProcessError:
                                    # Write in the logs that there was an error detected in the FASTQ files
                                    write_to_logfile('An error was detected in the FASTQ files for sample {}. '
                                                     'These files will not be processed further'.format(sample.name),
                                                     'An error was detected in the FASTQ files for sample {}. '
                                                     'These files will not be processed further'.format(sample.name),
                                                     self.logfile,
                                                     sample.general.logout,
                                                     sample.general.logerr,
                                                     None,
                                                     None)
                                    # Update metadata objects with error
                                    self.error(sample, 'fastq_error')
                            else:
                                # Write in the logs that there was an error detected in the FASTQ files
                                write_to_logfile('An error was detected in the FASTQ files for sample {}. '
                                                 'These files will not be processed further'.format(sample.name),
                                                 'An error was detected in the FASTQ files for sample {}. '
                                                 'These files will not be processed further'.format(sample.name),
                                                 self.logfile,
                                                 sample.general.logout,
                                                 sample.general.logerr,
                                                 None,
                                                 None)

                                # Update metadata objects with error
                                self.error(sample, 'fastq_error')
                else:
                    # Update metadata objects with error
                    self.error(sample, 'files_too_small')
        # Print the metadata to file
        metadataprinter.MetadataPrinter(self)
251426
def _cram_to_fastq_region(cram_file, work_dir, base_name, region, data):
    
    ref_file = tz.get_in(["reference", "fasta", "base"], data)
    resources = config_utils.get_resources("bamtofastq", data["config"])
    cores = tz.get_in(["config", "algorithm", "num_cores"], data, 1)
    max_mem = config_utils.convert_to_bytes(resources.get("memory", "1G")) * cores
    rext = "-%s" % region.replace(":", "_").replace("-", "_") if region else "full"
    out_s, out_p1, out_p2, out_o1, out_o2 = [os.path.join(work_dir, "%s%s-%s.fq.gz" %
                                                          (base_name, rext, fext))
                                             for fext in ["s1", "p1", "p2", "o1", "o2"]]
    if not utils.file_exists(out_p1):
        with file_transaction(data, out_s, out_p1, out_p2, out_o1, out_o2) as \
             (tx_out_s, tx_out_p1, tx_out_p2, tx_out_o1, tx_out_o2):
            cram_file = objectstore.cl_input(cram_file)
            sortprefix = "%s-sort" % utils.splitext_plus(tx_out_s)[0]
            cmd = ("bamtofastq filename={cram_file} inputformat=cram T={sortprefix} "
                   "gz=1 collate=1 colsbs={max_mem} exclude=SECONDARY,SUPPLEMENTARY "
                   "F={tx_out_p1} F2={tx_out_p2} S={tx_out_s} O={tx_out_o1} O2={tx_out_o2} "
                   "reference={ref_file}")
            if region:
                cmd += " ranges='{region}'"
            do.run(cmd.format(**locals()), "CRAM to fastq %s" % region if region else "")
    return [[out_p1, out_p2, out_s]]
1240585
def prep_bam_inputs(out_dir, sample, call_file, bam_file):
    
    base = utils.splitext_plus(os.path.basename(bam_file))[0]
    with open(call_file) as in_handle:
        for cur_hla in (x.strip() for x in in_handle):
            out_file = os.path.join(utils.safe_makedir(os.path.join(out_dir, base)),
                                    "%s.type.%s.filtered.bam" % (base, cur_hla))
            if not os.path.exists(out_file):
                cmd = ["samtools", "view", "-b","-o", out_file, bam_file, cur_hla]
                subprocess.check_call(cmd)
1068509
def fastqIterate(infile):
    ''''''

    def convert2string(b):
        if type(b) == str:
            return b
        else:
            return b.decode("utf-8")

    while 1:
        line1 = convert2string(infile.readline())
        if not line1:
            break
        if not line1.startswith('@'):
            U.error("parsing error: expected '@' in line %s" % line1)
        line2 = convert2string(infile.readline())
        line3 = convert2string(infile.readline())
        if not line3.startswith('+'):
            U.error("parsing error: expected '+' in line %s" % line3)
        line4 = convert2string(infile.readline())
        # incomplete entry
        if not line4:
            U.error("incomplete entry for %s" % line1)

        yield Record(line1[1:-1], line2[:-1], line4[:-1])
518751
def filter_reads(self, input_bam, output_bam, metrics_file, paired=False, cpus=16, Q=30):
        
        nodups = re.sub("\.bam$", "", output_bam) + ".nodups.nofilter.bam"
        cmd1 = self.tools.sambamba + " markdup -t {0} -r --compression-level=0 {1} {2} 2> {3}".format(cpus, input_bam, nodups, metrics_file)
        cmd2 = self.tools.sambamba + ' view -t {0} -f bam --valid'.format(cpus)
        if paired:
            cmd2 += ' -F "not (unmapped or mate_is_unmapped) and proper_pair'
        else:
            cmd2 += ' -F "not unmapped'
        cmd2 += ' and not (secondary_alignment or supplementary) and mapping_quality >= {0}"'.format(Q)
        cmd2 += ' {0} |'.format(nodups)
        cmd2 += self.tools.sambamba + " sort -t {0} /dev/stdin -o {1}".format(cpus, output_bam)
        cmd3 = "if [[ -s {0} ]]; then rm {0}; fi".format(nodups)
        cmd4 = "if [[ -s {0} ]]; then rm {0}; fi".format(nodups + ".bai")
        return [cmd1, cmd2, cmd3, cmd4]
351511
def piped_bamprep(data, region=None, out_file=None):
    
    data["region"] = region
    if not _need_prep(data):
        return [data]
    else:
        utils.safe_makedir(os.path.dirname(out_file))
        if region[0] == "nochrom":
            prep_bam = shared.write_nochr_reads(data["work_bam"], out_file, data["config"])
        elif region[0] == "noanalysis":
            prep_bam = shared.write_noanalysis_reads(data["work_bam"], region[1], out_file,
                                                     data["config"])
        else:
            if not utils.file_exists(out_file):
                with tx_tmpdir(data) as tmp_dir:
                    _piped_bamprep_region(data, region, out_file, tmp_dir)
            prep_bam = out_file
        bam.index(prep_bam, data["config"])
        data["work_bam"] = prep_bam
        return [data]
1354632
def run(bam_file, sample, out_dir):
    """"""
    out = {}
    # if "rchipqc" in dd.get_tools_on(sample):
    #    out = chipqc(bam_file, sample, out_dir)

    peaks = sample.get("peaks_files", {}).get("main")
    if peaks:
        out.update(_reads_in_peaks(bam_file, peaks, sample))
    return out
831980
def check_required_arguments(ribo_file, transcriptome_fasta, transcript_name=None):
    """"""
    # Is this a valid BAM file? i.e., can pysam read it?
    try:
        is_bam_valid(ribo_file)
    except ValueError:
        log.error('The given RiboSeq BAM file is not valid')
        raise

    # Does the BAM file have an index? If not, create it.
    if not bam_has_index(ribo_file):
        log.info('Creating an index for the BAM file...')
        create_bam_index(ribo_file)

        if not bam_has_index(ribo_file):
            msg = ('Could not create an index for this BAM file. Is this a valid BAM file '
                   'and/or is the BAM file sorted by chromosomal coordinates?')
            log.error(msg)
            raise BamFileError(msg)

    # Is FASTA file valid?
    fasta_valid = False
    try:
        fasta_valid = is_fasta_valid(transcriptome_fasta)
    except IOError:
        log.error('Transcriptome FASTA file is not valid')
        raise

    if fasta_valid:
        if transcript_name:
            try:
                get_fasta_records(transcriptome_fasta, [transcript_name])
            except IOError:
                log.error('Could not get FASTA sequence of "{}" from transcriptome FASTA file'.format(transcript_name))
                raise
        else:
            # ribocount doesn't have a transcript option so we get the first
            # sequence name from the fasta file
            transcript_name = get_first_transcript_name(transcriptome_fasta)

        # check if transcript also exists in BAM
        with pysam.AlignmentFile(ribo_file, 'rb') as bam_file:
            if transcript_name not in bam_file.references:
                msg = 'Transcript "{}" does not exist in BAM file'.format(transcript_name)
                log.error(msg)
                raise ArgumentError(msg)
1062248
def bam_to_fastq(bam_file, is_paired):
    
    out_files, out_handles = _get_fastq_handles(bam_file,
            is_paired)
    if len(out_handles) > 0:
        in_bam = pysam.Samfile(bam_file, mode='rb')
        for read in in_bam:
            num = 1 if (not read.is_paired or read.is_read1) else 2
            # reverse the sequence and quality if mapped to opposite strand
            if read.is_reverse:
                seq = str(Seq.reverse_complement(Seq.Seq(read.seq)))
                qual = "".join(reversed(read.qual))
            else:
                seq = read.seq
                qual = read.qual
            out_handles[num].write("@%s\n%s\n+\n%s\n" % (read.qname,
                seq, qual))
    [h.close() for h in out_handles.values()]
    return out_files
632148
def _piped_bamprep_region_gatk(data, region, prep_params, out_file, tmp_dir):
    
    broad_runner = broad.runner_from_config(data["config"])
    cur_bam, cl = _piped_input_cl(data, region, tmp_dir, out_file, prep_params)
    if not prep_params["realign"]:
        prerecal_bam = None
    elif prep_params["realign"] == "gatk":
        prerecal_bam, cl = _piped_realign_gatk(data, region, cl, out_file, tmp_dir,
                                               prep_params)
    else:
        raise NotImplementedError("Realignment method: %s" % prep_params["realign"])
    with file_transaction(data, out_file) as tx_out_file:
        out_flag = ("-o" if (prep_params["realign"] == "gatk"
                             or not prep_params["realign"])
                    else ">")
        cmd = "{cl} {out_flag} {tx_out_file}".format(**locals())
        do.run(cmd, "GATK: realign {0}".format(region), data)
        _cleanup_tempfiles(data, [cur_bam, prerecal_bam])
451199
def process_fastq_plain(fastq, **kwargs):
    """"""
    logging.info("Nanoget: Starting to collect statistics from plain fastq file.")
    inputfastq = handle_compressed_input(fastq)
    return ut.reduce_memory_usage(pd.DataFrame(
        data=[res for res in extract_from_fastq(inputfastq) if res],
        columns=["quals", "lengths"]
    ).dropna())
666804
def fastq(args):
    
    p = OptionParser(fastq.__doc__)
    opts, args = p.parse_args(args)

    if len(args) != 2:
        sys.exit(not p.print_help())

    bamfile, pf = args
    singletons = pf + ".se.fastq"
    a = pf + ".read1.fastq"
    b = pf + ".read2.fastq"

    cmd  = "samtools collate -uOn 128 {} tmp-prefix".format(bamfile)
    cmd += " | samtools fastq -s {} -1 {} -2 {} -"\
                .format(singletons, a, b)
    sh(cmd)

    if os.stat(singletons).st_size == 0:  # singleton file is empty
        os.remove(singletons)
    return a, b
814112
def collapse_fastq(args):
    
    try:
        umi_fn = args.fastq
        if _is_umi(args.fastq):
            umis = collapse(args.fastq)
            umi_fn = os.path.join(args.out, splitext_plus(os.path.basename(args.fastq))[0] + "_umi_trimmed.fastq")
            write_output(umi_fn, umis, args.minimum)
        seqs = collapse(umi_fn)
        out_file = splitext_plus(os.path.basename(args.fastq))[0] + "_trimmed.fastq"
    except IOError as e:
        logger.error("I/O error({0}): {1}".format(e.errno, e.strerror))
        raise "Can not read file"
    out_file = os.path.join(args.out, out_file)
    write_output(out_file, seqs, args.minimum)
    return out_file
861267
def _run_qc_tools(bam_file, data):
    
    from bcbio.qc import (atropos, contamination, coverage, damage, fastqc, kraken,
                          qsignature, qualimap, samtools, picard, srna, umi, variant,
                          viral, preseq, chipseq)
    tools = {"fastqc": fastqc.run,
             "atropos": atropos.run,
             "small-rna": srna.run,
             "samtools": samtools.run,
             "qualimap": qualimap.run,
             "qualimap_rnaseq": qualimap.run_rnaseq,
             "qsignature": qsignature.run,
             "contamination": contamination.run,
             "coverage": coverage.run,
             "damage": damage.run,
             "variants": variant.run,
             "peddy": peddy.run_qc,
             "kraken": kraken.run,
             "picard": picard.run,
             "umi": umi.run,
             "viral": viral.run,
             "preseq": preseq.run,
             "chipqc": chipseq.run
             }
    qc_dir = utils.safe_makedir(os.path.join(data["dirs"]["work"], "qc", data["description"]))
    metrics = {}
    qc_out = utils.deepish_copy(dd.get_summary_qc(data))
    for program_name in dd.get_algorithm_qc(data):
        if not bam_file and program_name != "kraken":  # kraken doesn't need bam
            continue
        if dd.get_phenotype(data) == "germline" and program_name != "variants":
            continue
        qc_fn = tools[program_name]
        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):
            # Check for metrics output, two cases:
            # 1. output with {"metrics"} and files ("base")
            if "metrics" in out:
                metrics.update(out.pop("metrics"))
            # 2. a dictionary of metrics
            elif "base" not in out:
                metrics.update(out)
            # Check for files only output
            if "base" in out:
                qc_files = out
        elif out and isinstance(out, six.string_types) 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

    metrics["Name"] = dd.get_sample_name(data)
    metrics["Quality format"] = dd.get_quality_format(data).lower()
    return {"qc": qc_out, "metrics": metrics}
1307327
def is_gzipped_fastq(file_name):
    
    _, ext = os.path.splitext(file_name)
    return file_name.endswith(".fastq.gz") or file_name.endswith(".fq.gz")