Demultiplex a fastqtransformed FASTQ file into a FASTQ file for
    each cell.
356599
def demultiplex_cells(fastq, out_dir, readnumber, prefix, cb_histogram, 
                      cb_cutoff):
    ''' 
    '''
    annotations = detect_fastq_annotations(fastq)
    re_string = construct_transformed_regex(annotations)
    parser_re = re.compile(re_string)
    readstring = "" if not readnumber else "_R{}".format(readnumber)
    filestring = "{prefix}{sample}{readstring}.fq"
    cb_set = set()
    if cb_histogram:
        cb_set = get_cb_depth_set(cb_histogram, cb_cutoff)
    sample_set = set()
    batch = collections.defaultdict(list)
    parsed = 0
    safe_makedir(out_dir)
    for read in read_fastq(fastq):
        parsed += 1
        match = parser_re.search(read).groupdict()
        sample = match['CB']
        if cb_set and sample not in cb_set:
            continue
        sample_set.add(sample)
        batch[sample].append(read)
        # write in batches to avoid opening up file handles repeatedly
        if not parsed % 10000000:
            for sample, reads in batch.items():
                out_file = os.path.join(out_dir, filestring.format(**locals()))
                with open(out_file, "a") as out_handle:
                    for read in reads:
                        out_handle.write(read)
            batch = collections.defaultdict(list)

    for sample, reads in batch.items():
        out_file = os.path.join(out_dir, filestring.format(**locals()))
        with open(out_file, "a") as out_handle:
            for read in reads:
                out_handle.write(read)
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())
732217
def _seqtk_fastq_prep_cl(data, in_file=None, read_num=0):
    
    needs_convert = dd.get_quality_format(data).lower() == "illumina"
    trim_ends = dd.get_trim_ends(data)
    seqtk = config_utils.get_program("seqtk", data["config"])
    if in_file:
        in_file = objectstore.cl_input(in_file)
    else:
        in_file = "/dev/stdin"
    cmd = ""
    if needs_convert:
        cmd += "{seqtk} seq -Q64 -V {in_file}".format(**locals())
    if trim_ends:
        left_trim, right_trim = trim_ends[0:2] if data.get("read_num", read_num) == 0 else trim_ends[2:4]
        if left_trim or right_trim:
            trim_infile = "/dev/stdin" if needs_convert else in_file
            pipe = " | " if needs_convert else ""
            cmd += "{pipe}{seqtk} trimfq -b {left_trim} -e {right_trim} {trim_infile}".format(**locals())
    return cmd
152376
def fastqWrite(fileHandleOrFile, name, seq, qualValues, mode="w"):
    
    fileHandle = _getFileHandle(fileHandleOrFile, mode)
    assert seq.__class__ == "".__class__
    for i in seq:
        if not ((i >= 'A' and i <= 'Z') or (i >= 'a' and i <= 'z') or i == '-'): #For safety and sanity I only allows roman alphabet characters in fasta sequences.
            raise RuntimeError("Invalid FASTQ character, ASCII code = \'%d\', char = '%s' found in input sequence %s" % (ord(i), i, name))
    if qualValues != None and qualValues != '*':
        if len(seq) != len(qualValues):
            raise RuntimeError("Got a mismatch between the number of sequence characters (%s) and number of qual values (%s) for sequence: %s " % (len(seq), len(qualValues), name))
        for i in qualValues:
            if i < 33 or i > 126:
                raise RuntimeError("Got a qual value out of range %s (range is 33 to 126)" % i)
        fileHandle.write("@%s\n%s\n+\n%s\n" % (name, seq, "".join([ chr(i) for i in qualValues ])))
    else:
        fileHandle.write("@%s\n%s\n+\n*\n" % (name, seq))
    if isinstance(fileHandleOrFile, "".__class__):
        fileHandle.close()
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
720725
def demultiplex_samples(data):
    
    work_dir = os.path.join(dd.get_work_dir(data), "umis")
    sample_dir = os.path.join(work_dir, dd.get_sample_name(data))
    demulti_dir = os.path.join(sample_dir, "demultiplexed")

    files = data["files"]
    if len(files) == 2:
        logger.error("Sample demultiplexing doesn't handle paired-end reads, but "
            "we can add it. Open an issue here https://github.com/bcbio/bcbio-nextgen/issues if you need this and we'll add it.")
        sys.exit(1)
    else:
        fq1 = files[0]
    # check if samples need to be demultiplexed
    with open_fastq(fq1) as in_handle:
        read = next(in_handle)
        if "SAMPLE_" not in read:
            return [[data]]

    bcfile = get_sample_barcodes(dd.get_sample_barcodes(data), sample_dir)
    demultiplexed = glob.glob(os.path.join(demulti_dir, "*.fq*"))
    if demultiplexed:
        return [split_demultiplexed_sampledata(data, demultiplexed)]
    umis = _umis_cmd(data)
    cmd = ("{umis} demultiplex_samples --nedit 1 --barcodes {bcfile} "
           "--out_dir {tx_dir} {fq1}")
    msg = "Demultiplexing {fq1}."
    with file_transaction(data, demulti_dir) as tx_dir:
        do.run(cmd.format(**locals()), msg.format(**locals()))
    demultiplexed = glob.glob(os.path.join(demulti_dir, "*.fq*"))
    return [split_demultiplexed_sampledata(data, demultiplexed)]
1084862
def _cutadapt_se_cmd(fastq_files, out_files, base_cmd, data):
    
    min_length = dd.get_min_read_length(data)
    cmd = base_cmd + " --minimum-length={min_length} ".format(**locals())
    fq1 = objectstore.cl_input(fastq_files[0])
    of1 = out_files[0]
    cmd += " -o {of1_tx} " + str(fq1)
    cmd = "%s | tee > {log_tx}" % cmd
    return cmd
3635
def movefastq(self):
        """"""
        logging.info('Moving FASTQ files')
        # Iterate through each sample
        for sample in self.metadata.runmetadata.samples:
            # Retrieve the output directory
            outputdir = os.path.join(self.path, sample.name)
            # Find any fastq files with the sample name
            fastqfiles = sorted(glob(os.path.join(self.path, '{}_*.fastq*'.format(sample.name)))) \
                if sorted(glob(os.path.join(self.path, '{}_*.fastq*'.format(sample.name)))) \
                else sorted(glob(os.path.join(self.path, '{}.fastq*'.format(sample.name)))) \
                if sorted(glob(os.path.join(self.path, '{}.fastq*'.format(sample.name)))) \
                else sorted(glob(os.path.join(self.path, '{}*.fastq*'.format(sample.name))))
            # Only try and move the files if the files exist
            if fastqfiles:
                make_path(outputdir)
                # Symlink the fastq files to the directory
                try:
                    list(map(lambda x: os.symlink(os.path.join('..', os.path.basename(x)),
                                                  os.path.join(outputdir, os.path.basename(x))), fastqfiles))
                except OSError:
                    pass
                # Find any fastq files with the sample name
                fastqfiles = [fastq for fastq in sorted(glob(os.path.join(outputdir, '{}*.fastq*'.format(sample.name))))
                              if 'trimmed' not in fastq and 'normalised' not in fastq and 'corrected' not in fastq
                              and 'paired' not in fastq and 'unpaired' not in fastq]
            else:
                if outputdir:
                    # Find any fastq files with the sample name
                    fastqfiles = [fastq for fastq in sorted(glob(os.path.join(
                        outputdir, '{}*.fastq*'.format(outputdir, sample.name))))
                                  if 'trimmed' not in fastq and 'normalised' not in fastq and 'corrected' not in fastq
                                  and 'paired' not in fastq and 'unpaired' not in fastq]
            sample.general.fastqfiles = fastqfiles
961002
def split_seq(iterable, size):
    """"""
    elm = iter(iterable)
    item = list(itertools.islice(elm, size))
    while item:
        yield item
        item = list(itertools.islice(elm, size))
466499
def fastqc(self):
        """"""
        while True:  # while daemon
            threadlock = threading.Lock()
            # Unpack the variables from the queue
            (sample, systemcall, outputdir, fastqcreads) = self.qcqueue.get()
            # Check to see if the output HTML file already exists
            try:
                _ = glob(os.path.join(outputdir, '*.html'))[0]
            except IndexError:
                # Make the output directory
                make_path(outputdir)
                # Run the system calls
                outstr = str()
                errstr = str()
                out, err = run_subprocess(systemcall)
                outstr += out
                errstr += err
                out, err = run_subprocess(fastqcreads)
                outstr += out
                errstr += err
                # Acquire thread lock, and write the logs to file
                threadlock.acquire()
                write_to_logfile(systemcall, systemcall, self.logfile, sample.general.logout, sample.general.logerr,
                                 None, None)
                write_to_logfile(fastqcreads, fastqcreads, self.logfile, sample.general.logout, sample.general.logerr,
                                 None, None)
                write_to_logfile(outstr, errstr, self.logfile, sample.general.logout, sample.general.logerr, None, None)
                threadlock.release()
                # Rename the outputs
                try:
                    shutil.move(os.path.join(outputdir, 'stdin_fastqc.html'),
                                os.path.join(outputdir, '{}_fastqc.html'.format(sample.name)))
                    shutil.move(os.path.join(outputdir, 'stdin_fastqc.zip'),
                                os.path.join(outputdir, '{}_fastqc.zip'.format(sample.name)))
                except IOError:
                    pass
            # Signal to qcqueue that job is done
            self.qcqueue.task_done()
39303
def fastqmover(self):
        """"""
        # Create the project path variable
        self.projectpath = self.fastqdestination + "/Project_" + self.projectname
        # Iterate through all the sample names
        for sample in self.metadata.samples:
            # Glob all the .gz files in the subfolders - projectpath/Sample_:sample.name/*.gz
            for fastq in sorted(glob('{}/Sample_{}/*.gz'.format(self.projectpath, sample.name))):
                # Try/except loop link .gz files to self.path
                try:
                    # Move fastq file to the path, but renames them first using the sample number.
                    move(
                        fastq, '{}{}'.format(self.path, os.path.basename(
                            sub('\w{8}-\w{8}', 'S{}'.format(
                                sample.run.SampleNumber), fastq))))
                # Except os errors
                except OSError as exception:
                    # If there is an exception other than the file exists, raise it
                    if exception.errno != errno.EEXIST:
                        raise
            # Repopulate .strainfastqfiles with the freshly-linked files
            fastqfiles = glob('{}/{}*.fastq*'.format(self.fastqdestination, sample.name))
            fastqfiles = [fastq for fastq in fastqfiles if 'trimmed' not in fastq]
            # Populate the metadata object with the name/path of the fastq files
            sample.general.fastqfiles = fastqfiles
            # Save the outputdir to the metadata object
            sample.run.outputdirectory = self.fastqdestination
        # Copy the sample sheet and the run info files to the path
        copyfile(self.assertions.samplesheet, os.path.join(self.path, 'SampleSheet.csv'))
        copyfile(os.path.join(self.miseqfolder, 'RunInfo.xml'), os.path.join(self.path, 'RunInfo.xml'))
678433
def aligned_read_to_read(read, revcomp=True, qual=None, ignore_quality=False):
    ''''''
    if read.qual is None or ignore_quality:
        if qual is None or ignore_quality:
            seq = pyfastaq.sequences.Fasta(read.qname, common.decode(read.seq))
        else:
            seq = pyfastaq.sequences.Fastq(read.qname, common.decode(read.seq), qual * read.query_length)
    else:
        if qual is None:
            seq = pyfastaq.sequences.Fastq(read.qname, common.decode(read.seq), common.decode(read.qual))
        else:
            seq = pyfastaq.sequences.Fastq(read.qname, common.decode(read.seq), qual * read.query_length)

    if read.is_reverse and revcomp:
        seq.revcomp()

    return seq
732951
def split_sequence_file_on_sample_ids_to_files(seqs,
                                               outdir):
    
    logger = logging.getLogger(__name__)
    logger.info('split_sequence_file_on_sample_ids_to_files'
                ' for file %s into dir %s' % (seqs, outdir))

    outputs = {}

    for bits in sequence_generator(seqs):
        sample = sample_id_from_read_id(bits[0])

        if sample not in outputs:
            outputs[sample] = open(join(outdir, sample + '.fasta'), 'w')
        outputs[sample].write(">%s\n%s\n" % (bits[0], bits[1]))
    for sample in outputs:
        outputs[sample].close()
    logger.info('split to %d files' % len(outputs))
1044511
def abi_to_fasta(input, output):
    '''
    
    '''
    direcs = [input, ]
    # unzip any zip archives
    zip_files = list_files(input, ['zip'])
    if zip_files:
        direcs.extend(_process_zip_files(zip_files))
    # convert files
    for d in direcs:
        files = list_files(d, ['ab1', 'abi'])
        seqs = [SeqIO.read(open(f, 'rb'), 'abi') for f in files]
        # seqs = list(chain.from_iterable(seqs))
        fastas = ['>{}\n{}'.format(s.id, str(s.seq)) for s in seqs]
        ofile = os.path.basename(os.path.normpath(d)) + '.fasta'
        opath = os.path.join(output, ofile)
        open(opath, 'w').write('\n'.join(fastas))
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]]
109282
def trimquality(self):
        """"""
        logging.info("Trimming fastq files")
        # Iterate through strains with fastq files
        with progressbar(self.metadata) as bar:
            for sample in bar:
                # As the metadata can be populated with 'NA' (string) if there are no fastq files, only process if
                # :fastqfiles is a list
                if type(sample.general.fastqfiles) is list:
                    # Check to see if the fastq files exist
                    fastqfiles = sorted(sample.general.fastqfiles)
                    # Define the output directory
                    outputdir = sample.general.outputdirectory
                    # Define the name of the trimmed fastq files
                    cleanforward = os.path.join(outputdir, '{}_R1_trimmed.fastq.gz'.format(sample.name))
                    cleanreverse = os.path.join(outputdir, '{}_R2_trimmed.fastq.gz'.format(sample.name))
                    # Incorporate read length into the minlength parameter - set it to 50 unless one or more of the
                    # reads has a lower calculated length than 50
                    try:
                        lesser_length = min(int(sample.run.forwardlength), int(sample.run.reverselength))
                    except ValueError:
                        lesser_length = int(sample.run.forwardlength)
                    min_len = 50 if lesser_length >= 50 else lesser_length
                    # Initialise a variable to store the number of bases to automatically trim from the beginning of
                    # each read, as these bases tend to have lower quality scores. If trimming the reads will cause
                    trim_left = 0
                    # If, for some reason, only the reverse reads are present, use the appropriate output file name
                    try:
                        if 'R2' in fastqfiles[0]:
                            if not os.path.isfile(cleanreverse):
                                out, \
                                    err, \
                                    bbdukcall = bbtools.bbduk_trim(forward_in=fastqfiles[0],
                                                                   reverse_in=None,
                                                                   forward_out=cleanreverse,
                                                                   trimq=10,
                                                                   minlength=min_len,
                                                                   forcetrimleft=trim_left,
                                                                   returncmd=True)
                            else:
                                bbdukcall = str()
                                out = str()
                                err = str()
                        else:
                            if not os.path.isfile(cleanforward):
                                    out, \
                                        err, \
                                        bbdukcall = bbtools.bbduk_trim(forward_in=fastqfiles[0],
                                                                       forward_out=cleanforward,
                                                                       trimq=10,
                                                                       minlength=min_len,
                                                                       forcetrimleft=trim_left,
                                                                       returncmd=True)
                            else:
                                bbdukcall = str()
                                out = str()
                                err = str()
                    except (IndexError, CalledProcessError):
                        bbdukcall = str()
                        out = str()
                        err = str()
                    # Write the command, stdout, and stderr to the logfile
                    write_to_logfile(bbdukcall, bbdukcall, self.logfile, sample.general.logout, sample.general.logerr,
                                     None, None)
                    write_to_logfile(out, err, self.logfile, sample.general.logout, sample.general.logerr, None, None)
                    # Add the trimmed fastq files to a list
                    trimmedfastqfiles = sorted(glob(os.path.join(sample.general.outputdirectory, '*trimmed.fastq.gz')))
                    # Populate the metadata if the files exist
                    sample.general.trimmedfastqfiles = trimmedfastqfiles if trimmedfastqfiles else list()
        # Add all the trimmed files to the metadata
        logging.info('Fastq files trimmed')
630803
def kallisto_table(kallisto_dir, index):
    
    quant_dir = os.path.join(kallisto_dir, "quant")
    out_file = os.path.join(quant_dir, "matrix.csv")
    if file_exists(out_file):
        return out_file
    tsvfile = os.path.join(quant_dir, "matrix.tsv")
    ecfile = os.path.join(quant_dir, "matrix.ec")
    cellsfile = os.path.join(quant_dir, "matrix.cells")
    fastafile = os.path.splitext(index)[0] + ".fa"
    fasta_names = fasta.sequence_names(fastafile)
    ec_names = get_ec_names(ecfile, fasta_names)
    df = pd.read_table(tsvfile, header=None, names=["ec", "cell", "count"])
    df["ec"] = [ec_names[x] for x in df["ec"]]
    df = df.pivot(index='ec', columns='cell', values='count')
    cellnames = get_cell_names(cellsfile)
    colnames = [cellnames[x] for x in df.columns]
    df.columns = colnames
    df.to_csv(out_file)
    return out_file