Convert a fastq or BAM input into a SDF indexed file.
356819
def to_sdf(files, data):
    
    # BAM
    if len(files) == 1 and files[0].endswith(".bam"):
        qual = []
        format = ["-f", "sam-pe" if bam.is_paired(files[0]) else "sam-se"]
        inputs = [files[0]]
    # fastq
    else:
        qual = ["-q", "illumina" if dd.get_quality_format(data).lower() == "illumina" else "sanger"]
        format = ["-f", "fastq"]
        if len(files) == 2:
            inputs = ["-l", files[0], "-r", files[1]]
        else:
            assert len(files) == 1
            inputs = [files[0]]
    work_dir = utils.safe_makedir(os.path.join(data["dirs"]["work"], "align_prep"))
    out_file = os.path.join(work_dir,
                            "%s.sdf" % utils.splitext_plus(os.path.basename(os.path.commonprefix(files)))[0])
    if not utils.file_exists(out_file):
        with file_transaction(data, out_file) as tx_out_file:
            cmd = _rtg_cmd(["rtg", "format", "-o", tx_out_file] + format + qual + inputs)
            do.run(cmd, "Format inputs to indexed SDF")
    return out_file
1300715
def _bgzip_from_bam(bam_file, dirs, data, is_retry=False, output_infix=''):
    
    # tools
    config = data["config"]
    bamtofastq = config_utils.get_program("bamtofastq", config)
    resources = config_utils.get_resources("bamtofastq", config)
    cores = config["algorithm"].get("num_cores", 1)
    max_mem = config_utils.convert_to_bytes(resources.get("memory", "1G")) * cores
    bgzip = tools.get_bgzip_cmd(config, is_retry)
    # files
    work_dir = utils.safe_makedir(os.path.join(dirs["work"], "align_prep"))
    out_file_1 = os.path.join(work_dir, "%s%s-1.fq.gz" % (os.path.splitext(os.path.basename(bam_file))[0], output_infix))
    out_file_2 = out_file_1.replace("-1.fq.gz", "-2.fq.gz")
    needs_retry = False
    if is_retry or not utils.file_exists(out_file_1):
        if not bam.is_paired(bam_file):
            out_file_2 = None
        with file_transaction(config, out_file_1) as tx_out_file:
            for f in [tx_out_file, out_file_1, out_file_2]:
                if f and os.path.exists(f):
                    os.remove(f)
            fq1_bgzip_cmd = "%s -c /dev/stdin > %s" % (bgzip, tx_out_file)
            prep_cmd = _seqtk_fastq_prep_cl(data, read_num=0)
            if prep_cmd:
                fq1_bgzip_cmd = prep_cmd + " | " + fq1_bgzip_cmd
            sortprefix = "%s-sort" % os.path.splitext(tx_out_file)[0]
            if bam.is_paired(bam_file):
                prep_cmd = _seqtk_fastq_prep_cl(data, read_num=1)
                fq2_bgzip_cmd = "%s -c /dev/stdin > %s" % (bgzip, out_file_2)
                if prep_cmd:
                    fq2_bgzip_cmd = prep_cmd + " | " + fq2_bgzip_cmd
                out_str = ("F=>({fq1_bgzip_cmd}) F2=>({fq2_bgzip_cmd}) S=/dev/null O=/dev/null "
                           "O2=/dev/null collate=1 colsbs={max_mem}")
            else:
                out_str = "S=>({fq1_bgzip_cmd})"
            bam_file = objectstore.cl_input(bam_file)
            extra_opts = " ".join([str(x) for x in resources.get("options", [])])
            cmd = "{bamtofastq} filename={bam_file} T={sortprefix} {extra_opts} " + out_str
            try:
                do.run(cmd.format(**locals()), "BAM to bgzipped fastq",
                       checks=[do.file_reasonable_size(tx_out_file, bam_file)],
                       log_error=False)
            except subprocess.CalledProcessError as msg:
                if not is_retry and "deflate failed" in str(msg):
                    logger.info("bamtofastq deflate IO failure preparing %s. Retrying with single core."
                                % (bam_file))
                    needs_retry = True
                else:
                    logger.exception()
                    raise
    if needs_retry:
        return _bgzip_from_bam(bam_file, dirs, data, is_retry=True)
    else:
        return [x for x in [out_file_1, out_file_2] if x is not None and utils.file_exists(x)]
254349
def sam_conversions(self, sam_file, depth=True):
        
        cmd = self.tools.samtools + " view -bS " + sam_file + " > " + sam_file.replace(".sam", ".bam") + "\n"
        cmd += self.tools.samtools + " sort " + sam_file.replace(".sam", ".bam") + " -o " + sam_file.replace(".sam", "_sorted.bam") + "\n"
        cmd += self.tools.samtools + " index " + sam_file.replace(".sam", "_sorted.bam") + "\n"
        if depth:
            cmd += self.tools.samtools + " depth " + sam_file.replace(".sam", "_sorted.bam") + " > " + sam_file.replace(".sam", "_sorted.depth") + "\n"
        return cmd
420415
def _stdio_stream_as(src_fd, dst_fd, dst_sys_attribute, mode):
  """"""
  if src_fd == -1:
    src = open('/dev/null', mode)
    src_fd = src.fileno()

  # Capture the python and os level file handles.
  old_dst = getattr(sys, dst_sys_attribute)
  old_dst_fd = os.dup(dst_fd)
  if src_fd != dst_fd:
    os.dup2(src_fd, dst_fd)

  # Open up a new file handle to temporarily replace the python-level io object, then yield.
  new_dst = os.fdopen(dst_fd, mode)
  setattr(sys, dst_sys_attribute, new_dst)
  try:
    yield
  finally:
    new_dst.close()

    # Restore the python and os level file handles.
    os.dup2(old_dst_fd, dst_fd)
    setattr(sys, dst_sys_attribute, old_dst)
574342
def _convert_fastq(srafn, outdir, single=False):
    ""
    cmd = "fastq-dump --split-files --gzip {srafn}"
    cmd = "%s %s" % (utils.local_path_export(), cmd)
    sraid = os.path.basename(utils.splitext_plus(srafn)[0])
    if not srafn:
        return None
    if not single:
        out_file = [os.path.join(outdir, "%s_1.fastq.gz" % sraid),
                    os.path.join(outdir, "%s_2.fastq.gz" % sraid)]
        if not utils.file_exists(out_file[0]):
            with utils.chdir(outdir):
                do.run(cmd.format(**locals()), "Covert to fastq %s" % sraid)
        if not utils.file_exists(out_file[0]):
            raise IOError("SRA %s didn't convert, something happened." % srafn)
        return [out for out in out_file if utils.file_exists(out)]
    else:
        raise ValueError("Not supported single-end sra samples for now.")
736957
def _merge_list_fastqs(files, out_file, config):
    """"""
    if not all(map(fastq.is_fastq, files)):
        raise ValueError("Not all of the files to merge are fastq files: %s " % (files))
    assert all(map(utils.file_exists, files)), ("Not all of the files to merge "
                                                "exist: %s" % (files))
    if not file_exists(out_file):
        files = [_gzip_fastq(fn) for fn in files]
        if len(files) == 1:
            if "remove_source" in config and config["remove_source"]:
                shutil.move(files[0], out_file)
            else:
                os.symlink(files[0], out_file)
            return out_file
        with file_transaction(out_file) as file_txt_out:
            files_str = " ".join(list(files))
            cmd = "cat {files_str} > {file_txt_out}".format(**locals())
            do.run(cmd, "merge fastq files %s" % files)
    return out_file
1004328
def _align_from_fastq(fastq1, fastq2, aligner, align_ref, sam_ref, names,
                      align_dir, data):
    
    config = data["config"]
    align_fn = TOOLS[aligner].align_fn
    out = align_fn(fastq1, fastq2, align_ref, names, align_dir, data)
    # handle align functions that update the main data dictionary in place
    if isinstance(out, dict):
        assert out.get("work_bam"), (dd.get_sample_name(data), out.get("work_bam"))
        return out
    # handle output of raw SAM files that need to be converted to BAM
    else:
        work_bam = bam.sam_to_bam(out, config)
        data["work_bam"] = bam.sort(work_bam, config)
        return data
233451
def load_data(self, idx):
        """"""
        for name in self:
            ndim = getattr(self, '_%s_ndim' % name)
            diskflag = getattr(self, '_%s_diskflag' % name)
            ramflag = getattr(self, '_%s_ramflag' % name)
            if diskflag:
                file_ = getattr(self, '_%s_file' % name)
                length_tot = 1
                shape = []
                for jdx in range(ndim):
                    length = getattr(self, '_%s_length_%s' % (name, jdx))
                    length_tot *= length
                    shape.append(length)
                raw = file_.read(length_tot*8)
                values = struct.unpack(length_tot*'d', raw)
                if ndim:
                    values = numpy.array(values).reshape(shape)
                else:
                    values = values[0]
            elif ramflag:
                array = getattr(self, '_%s_array' % name)
                values = array[idx]
            if diskflag or ramflag:
                if ndim == 0:
                    setattr(self, name, values)
                else:
                    getattr(self, name)[:] = values
849285
def tohdf5(input_files, output_file, n_events, conv_times_to_jte, **kwargs):
    """"""
    if len(input_files) > 1:
        cprint(
            "Preparing to convert {} files to HDF5.".format(len(input_files))
        )

    from km3pipe import Pipeline    # noqa
    from km3pipe.io import GenericPump, HDF5Sink, HDF5MetaData    # noqa

    for input_file in input_files:
        cprint("Converting '{}'...".format(input_file))
        if len(input_files) > 1:
            output_file = input_file + '.h5'

        meta_data = kwargs.copy()
        meta_data['origin'] = input_file

        pipe = Pipeline()
        pipe.attach(HDF5MetaData, data=meta_data)
        pipe.attach(GenericPump, filenames=input_file, **kwargs)
        pipe.attach(StatusBar, every=250)
        if conv_times_to_jte:
            from km3modules.mc import MCTimeCorrector
            pipe.attach(MCTimeCorrector)
        pipe.attach(HDF5Sink, filename=output_file, **kwargs)
        pipe.drain(n_events)
        cprint("File '{}' was converted.".format(input_file))
814889
def get_fastq_files(data):
    
    assert "files" in data, "Did not find `files` in input; nothing to process"
    ready_files = []
    should_gzip = True

    # Bowtie does not accept gzipped fastq
    if 'bowtie' in data['reference'].keys():
        should_gzip = False
    for fname in data["files"]:
        if fname.endswith(".bam"):
            if _pipeline_needs_fastq(data["config"], data):
                ready_files = convert_bam_to_fastq(fname, data["dirs"]["work"],
                                                   data, data["dirs"], data["config"])
            else:
                ready_files = [fname]
        elif objectstore.is_remote(fname):
            ready_files.append(fname)
        # Trimming does quality conversion, so if not doing that, do an explicit conversion
        elif not(dd.get_trim_reads(data)) and dd.get_quality_format(data) != "standard":
            out_dir = utils.safe_makedir(os.path.join(dd.get_work_dir(data), "fastq_convert"))
            ready_files.append(fastq.groom(fname, data, out_dir=out_dir))
        else:
            ready_files.append(fname)
    ready_files = [x for x in ready_files if x is not None]
    if should_gzip:
        out_dir = utils.safe_makedir(os.path.join(dd.get_work_dir(data), "fastq"))
        ready_files = [_gzip_fastq(x, out_dir) for x in ready_files]
    for in_file in ready_files:
        if not objectstore.is_remote(in_file):
            assert os.path.exists(in_file), "%s does not exist." % in_file
    return ready_files
1038266
def align(fastq_file, pair_file, index_dir, names, align_dir, data):
    
    umi_ext = "-cumi" if "umi_bam" in data else ""
    out_file = os.path.join(align_dir, "{0}-sort{1}.bam".format(dd.get_sample_name(data), umi_ext))
    num_cores = data["config"]["algorithm"].get("num_cores", 1)
    rg_info = "rgid={rg} rgpl={pl} rgpu={pu} rgsm={sample}".format(**names)
    pair_file = pair_file if pair_file else ""
    final_file = None
    if data.get("align_split"):
        # BBMap does not accept input fastq streams
        raise ValueError("bbmap is not compatible with alignment splitting, set `align_split: false`")
    pair_arg = "in2=%s" % pair_file if pair_file else ""
    if not utils.file_exists(out_file) and (final_file is None or not utils.file_exists(final_file)):
        with postalign.tobam_cl(data, out_file, pair_file != "") as (tobam_cl, tx_out_file):
            if index_dir.endswith(("/ref", "/ref/")):
                index_dir = os.path.dirname(index_dir)
            # sam=1.3 required for compatibility with strelka2
            cmd = ("bbmap.sh sam=1.3 mdtag=t {rg_info} path={index_dir} in1={fastq_file} "
                   "{pair_arg} out=stdout.sam | ")
            do.run(cmd.format(**locals()) + tobam_cl, "bbmap alignment: %s" % dd.get_sample_name(data))
    data["work_bam"] = out_file
    return data
706977
def _prepare_inputs(ma_fn, bam_file, out_dir):
    
    fixed_fa = os.path.join(out_dir, "file_reads.fa")
    count_name =dict()
    with file_transaction(fixed_fa) as out_tx:
        with open(out_tx, 'w') as out_handle:
            with open(ma_fn) as in_handle:
                h = next(in_handle)
                for line in in_handle:
                    cols = line.split("\t")
                    name_with_counts = "%s_x%s" % (cols[0], sum(map(int, cols[2:])))
                    count_name[cols[0]] = name_with_counts
                    out_handle.write(">%s\n%s\n" % (name_with_counts, cols[1]))
    fixed_bam = os.path.join(out_dir, "align.bam")
    bam_handle = pysam.AlignmentFile(bam_file, "rb")
    with pysam.AlignmentFile(fixed_bam, "wb", template=bam_handle) as out_handle:
        for read in bam_handle.fetch():
            read.query_name = count_name[read.query_name]
            out_handle.write(read)

    return fixed_fa, fixed_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())
634899
def _merge_fastqc(samples):
    
    fastqc_list = collections.defaultdict(list)
    seen = set()
    for data in samples:
        name = dd.get_sample_name(data)
        if name in seen:
            continue
        seen.add(name)
        fns = glob.glob(os.path.join(dd.get_work_dir(data), "qc", dd.get_sample_name(data), "fastqc") + "/*")
        for fn in fns:
            if fn.endswith("tsv"):
                metric = os.path.basename(fn)
                fastqc_list[metric].append([name, fn])

    for metric in fastqc_list:
        dt_by_sample = []
        for fn in fastqc_list[metric]:
            dt = pd.read_csv(fn[1], sep="\t")
            dt['sample'] = fn[0]
            dt_by_sample.append(dt)
        dt = utils.rbind(dt_by_sample)
        dt.to_csv(metric, sep="\t", index=False, mode ='w')
    return samples
1294184
def process_fasta(fasta, **kwargs):
    """"""
    logging.info("Nanoget: Starting to collect statistics from a fasta file.")
    inputfasta = handle_compressed_input(fasta, file_type="fasta")
    return ut.reduce_memory_usage(pd.DataFrame(
        data=[len(rec) for rec in SeqIO.parse(inputfasta, "fasta")],
        columns=["lengths"]
    ).dropna())
1020906
private static byte[] extractBytes(
      String str, 
      int startIndex, 
      int endIndex,
      char[] charArray, 
      boolean canFastConvert) throws UnsupportedEncodingException {
    if (canFastConvert) {
      // fast conversion, just copy the raw bytes
      final int len = endIndex - startIndex;
      byte[] strBytes = new byte[len];
      for (int i = 0; i < len; i++) {
        strBytes[i] = (byte) charArray[startIndex + i];
      }
      return strBytes;
    } 
    // otherwise, do expensive conversion
    return str.substring(startIndex, endIndex).getBytes(utf8charsetName);
  }
1068733
def from_json(buffer):
        """"""
        buffer = to_bytes(buffer)
        return Index._from_ptr(rustcall(
            _lib.lsm_index_from_json,
            buffer, len(buffer)))