Convert a fastq file into a CSV of phred quality scores.
548696
def fastq_to_csv(in_file, fastq_format, work_dir):
    
    out_file = "%s.csv" % (os.path.splitext(os.path.basename(in_file))[0])
    out_file = os.path.join(work_dir, out_file)
    if not (os.path.exists(out_file) and os.path.getsize(out_file) > 0):
        with open(in_file) as in_handle:
            with open(out_file, "w") as out_handle:
                writer = csv.writer(out_handle)
                for rec in SeqIO.parse(in_handle, fastq_format):
                    writer.writerow([rec.id] + rec.letter_annotations["phred_quality"])
    return out_file
900989
def _gzip_fastq(in_file, out_dir=None):
    
    if fastq.is_fastq(in_file) and not objectstore.is_remote(in_file):
        if utils.is_bzipped(in_file):
            return _bzip_gzip(in_file, out_dir)
        elif not utils.is_gzipped(in_file):
            if out_dir:
                gzipped_file = os.path.join(out_dir, os.path.basename(in_file) + ".gz")
            else:
                gzipped_file = in_file + ".gz"
            if file_exists(gzipped_file):
                return gzipped_file
            message = "gzipping {in_file} to {gzipped_file}.".format(
                in_file=in_file, gzipped_file=gzipped_file)
            with file_transaction(gzipped_file) as tx_gzipped_file:
                do.run("gzip -c {in_file} > {tx_gzipped_file}".format(**locals()),
                       message)
            return gzipped_file
    return in_file
676770
def write(self):
        """"""
        outfn = self.create_outfilepath(self.fn, self.outsuffix)
        command = ['qvality']
        command.extend(self.qvalityoptions)
        command.extend([self.scores['target']['fn'], self.scores['decoy']['fn'],
                        '-o', outfn])
        subprocess.call(command)
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())
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
345429
def write_prepareref_fasta_file(outfile, gene_coords, genes_need_upstream, genes_non_upstream, upstream_before=100, upstream_after=100):
    ''''''
    tmp_dict = {}
    fasta_in = os.path.join(data_dir, 'NC_000962.3.fa.gz')
    pyfastaq.tasks.file_to_dict(fasta_in, tmp_dict)
    ref_seq = tmp_dict['NC_000962.3']

    with open(outfile, 'w') as f:
        for gene in genes_non_upstream:
            start = gene_coords[gene]['start']
            end = gene_coords[gene]['end']
            if start < end:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[start:end+1])
            else:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[end:start+1])
                gene_fa.revcomp()

            print(gene_fa, file=f)

        for gene in genes_need_upstream:
            start = gene_coords[gene]['start']
            end = gene_coords[gene]['end']
            if start < end:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[start - upstream_before:start + upstream_after])
            else:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[start - upstream_after + 1:start + upstream_before + 1])
                gene_fa.revcomp()

            gene_fa.id += '_upstream'
            print(gene_fa, file=f)
713935
def kallisto(fastq, out_dir, cb_histogram, cb_cutoff):
    ''' 
    '''
    parser_re = re.compile('(.*):CELL_(?<CB>.*):UMI_(?P<UMI>.*)\\n(.*)\\n\\+\\n(.*)\\n')
    if fastq.endswith('gz'):
        fastq_fh = gzip.GzipFile(fileobj=open(fastq))
    elif fastq == "-":
        fastq_fh = sys.stdin
    else:
        fastq_fh = open(fastq)

    cb_depth_set = get_cb_depth_set(cb_histogram, cb_cutoff)

    cb_set = set()
    cb_batch = collections.defaultdict(list)
    parsed = 0
    for read in stream_fastq(fastq_fh):
        match = parser_re.search(read).groupdict()
        umi = match['UMI']
        cb = match['CB']
        if cb_depth_set and cb not in cb_depth_set:
            continue
        parsed += 1
        cb_set.add(cb)
        cb_batch[cb].append((read, umi))
        # write in batches to avoid opening up file handles repeatedly
        if not parsed % 10000000:
            for cb, chunk in cb_batch.items():
                write_kallisto_chunk(out_dir, cb, chunk)
            cb_batch = collections.defaultdict(list)

    for cb, chunk in cb_batch.items():
        write_kallisto_chunk(out_dir, cb, chunk)

    with open(os.path.join(out_dir, "barcodes.batch"), "w") as out_handle:
        out_handle.write("#id umi-file file-1\n")
        batchformat = "{cb} {cb}.umi {cb}.fq\n"
        for cb in cb_set:
            out_handle.write(batchformat.format(**locals()))
1003084
def createfastq(self):
        """"""
        # Initialise samplecount
        samplecount = 0
        # If the fastq destination folder is not provided, make the default value of :path/:miseqfoldername
        self.fastqdestination = self.fastqdestination if self.fastqdestination else self.path + self.miseqfoldername
        # Make the path
        make_path(self.fastqdestination)
        # Initialise variables for storing index information
        index = ''
        indexlength = int()
        # bcl2fastq requires an older version of the sample sheet, this recreates the required version
        # Create the new sample sheet
        with open('{}/SampleSheet_modified.csv'.format(self.fastqdestination), "w") as modifiedsamplesheet:
            # Write the required headings to the file
            modifiedsamplesheet.write(
                "FCID,Lane,SampleID,SampleRef,Index,Description,Control,Recipe,Operator,SampleProject\n")
            for strain in self.samples:
                # Create a combined index of index1-index2
                try:
                    strain.run.modifiedindex = '{}-{}'.format(strain.run.index, strain.run.index2)
                    indexlength = 16
                    index = 'I8,I8'
                except KeyError:
                    strain.run.modifiedindex = strain.run.index
                    indexlength = 6
                    index = 'I6'
                # The list of items to print to each line of the modified sample sheet
                printlist = [self.flowcell, '1', strain.name, str(strain.run.SampleNumber), strain.run.modifiedindex,
                             strain.run.Description, 'N', 'NA',
                             strain.run.InvestigatorName, self.projectname]
                modifiedsamplesheet.write('{}\n'.format(",".join(printlist)))
                samplecount += 1
        # Set :forward/reverse length to :header.forward/reverse length if the argument is not provided, or it's 'full',
        # otherwise  use the supplied argument
        self.forwardlength = self.metadata.header.forwardlength if self.forwardlength.lower()\
            == 'full' else self.forwardlength
        # Set :reverselength to :header.reverselength
        self.reverselength = self.metadata.header.reverselength if self.reverselength.lower() \
            == 'full' else self.reverselength
        # As the number of cycles required is the number of forward reads + the index(8) + the second index(8)
        # Also set the basemask variable as required
        if self.reverselength != '0':
            self.readsneeded = int(self.forwardlength) + int(self.reverselength) + indexlength
            basemask = "Y{}n*,{},Y{}n*".format(self.forwardlength, index, self.reverselength)
            nohup = "nohup make -j 16 > nohup.out"
        else:
            #  + 1
            self.readsneeded = int(self.forwardlength) + indexlength
            basemask = "Y{}n*,{},n*".format(self.forwardlength, index)
            nohup = "nohup make -j 16 r1 > nohup.out"
        # Handle plurality appropriately
        samples = 'samples' if samplecount > 1 else 'sample'
        number = 'are' if samplecount > 1 else 'is'
        printtime('There {} {} {} in this run. '
                  'Running fastq creating module with the following parameters:\n'
                  'MiSeqPath: {},\n'
                  'MiSeqFolder: {},\n'
                  'Fastq destination: {},\n'
                  'SampleSheet: {}'
                  .format(number, samplecount, samples, self.miseqpath, self.miseqfolder,
                          self.fastqdestination, '{}/SampleSheet_modified.csv'.format(self.fastqdestination)),
                  self.start)
        # Count the number of completed cycles in the run of interest
        cycles = glob('{}Data/Intensities/BaseCalls/L001/C*'.format(self.miseqfolder))
        while len(cycles) < self.readsneeded:
            printtime('Currently at {} cycles. Waiting until the MiSeq reaches cycle {}'.format(len(cycles),
                      self.readsneeded), self.start)
            sleep(1800)
            cycles = glob('{}Data/Intensities/BaseCalls/L001/C*'.format(self.miseqfolder))
        # configureBClToFastq requires :self.miseqfolder//Data/Intensities/BaseCalls/config.xml in order to work
        # When you download runs from BaseSpace, this file is not provided. There is an empty config.xml file that
        # can be populated with run-specific values and moved to the appropriate folder
        if not os.path.isfile('{}Data/Intensities/BaseCalls/config.xml'.format(self.miseqfolder)):
            self.configfilepopulator()
        # Define the bcl2fastq system call
        bclcall = "configureBclToFastq.pl --input-dir {}Data/Intensities/BaseCalls " \
                  "--output-dir {} --force --sample-sheet {}/SampleSheet_modified.csv " \
                  "--mismatches 1 --no-eamss --fastq-cluster-count 0 --compression none --use-bases-mask {}"\
            .format(self.miseqfolder, self.fastqdestination, self.fastqdestination, basemask)
        # Define the nohup system call
        nohupcall = "cd {} && {}".format(self.fastqdestination, nohup)
        # fnull = open(os.devnull, 'wb')
        if not os.path.isdir("{}/Project_{}".format(self.fastqdestination, self.projectname)):
            # Call configureBclToFastq.pl
            printtime('Running bcl2fastq', self.start)
            # Run the commands
            threadlock = threading.Lock()
            outstr = ''
            outerr = ''
            out, err = run_subprocess(bclcall)
            outstr += out
            outerr += out
            out, err = run_subprocess(nohupcall)
            outstr += out
            outerr += out
            # call(bclcall, shell=True, stdout=fnull, stderr=fnull)
            # call(nohupcall, shell=True, stdout=fnull, stderr=fnull)
            threadlock.acquire()
            write_to_logfile(bclcall, bclcall, self.logfile)
            write_to_logfile(nohupcall, nohupcall, self.logfile)
            write_to_logfile(outstr, outerr, self.logfile)
            threadlock.release()
        # Populate the metadata
        for sample in self.metadata.samples:
            sample.commands = GenObject()
            sample.commands.nohup = nohupcall
            sample.commands.bcl = bclcall
            sample.run.forwardlength = self.forwardlength
            sample.run.reverselength = self.reverselength
        # Copy the fastq files to a central folder so they can be processed
        self.fastqmover()
1146061
def _cutadapt_trim(fastq_files, quality_format, adapters, out_files, log_file, data):
    
    if all([utils.file_exists(x) for x in out_files]):
        return out_files
    cmd = _cutadapt_trim_cmd(fastq_files, quality_format, adapters, out_files, data)
    if len(fastq_files) == 1:
        of = [out_files[0], log_file]
        message = "Trimming %s in single end mode with cutadapt." % (fastq_files[0])
        with file_transaction(data, of) as of_tx:
            of1_tx, log_tx = of_tx
            do.run(cmd.format(**locals()), message)
    else:
        of = out_files + [log_file]
        with file_transaction(data, of) as tx_out_files:
            of1_tx, of2_tx, log_tx = tx_out_files
            tmp_fq1 = utils.append_stem(of1_tx, ".tmp")
            tmp_fq2 = utils.append_stem(of2_tx, ".tmp")
            singles_file = of1_tx + ".single"
            message = "Trimming %s and %s in paired end mode with cutadapt." % (fastq_files[0],
                                                                                fastq_files[1])
            do.run(cmd.format(**locals()), message)
    return out_files
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'))
661519
def format(file_metrics, build_metrics):
    """"""
    # TODO maybe we need different output for build_metrics in csv format, too?
    # filter out positions metric
    def report_header(file_metrics):
        values = list(file_metrics.values())[0]
        print(values)
        values.pop('block_positions', None)
        return 'filename,' + ','.join(values) + '\n'

    def report_metrics(file_metrics):
        report = ''
        for key, values in file_metrics.items():
            report += key + ','
            report += ','.join([str(v) for k, v in values.items() if k not in ['block_positions']])
            report += '\n'
        return report

    report = report_header(file_metrics)
    report += report_metrics(file_metrics)
    return report
242949
def to_fastq_str(self):
    
    return "@" + self.name + "\n" + self.sequenceData +\
           "\n" + "+" + self.name + "\n" + self.seq_qual
423875
def fastq(args):
    
    from jcvi.formats.fastq import guessoffset

    p = OptionParser(fastq.__doc__)
    p.add_option("--outtie", dest="outtie", default=False, action="store_true",
            help="Are these outie reads? [default: %default]")
    p.set_phred()
    p.set_size()

    opts, args = p.parse_args(args)

    if len(args) < 1:
        sys.exit(p.print_help())

    fastqfiles = [get_abs_path(x) for x in args]
    size = opts.size
    outtie = opts.outtie
    if size > 1000 and (not outtie):
        logging.debug("[warn] long insert size {0} but not outtie".format(size))

    mated = (size != 0)
    libname = op.basename(args[0]).split(".")[0]
    libname = libname.replace("_1_sequence", "")

    frgfile = libname + ".frg"
    mean, sv = get_mean_sv(opts.size)

    cmd = "fastqToCA"
    cmd += " -libraryname {0} ".format(libname)
    fastqs = " ".join("-reads {0}".format(x) for x in fastqfiles)
    if mated:
        assert len(args) in (1, 2), "you need one or two fastq files for mated library"
        fastqs = "-mates {0}".format(",".join(fastqfiles))
        cmd += "-insertsize {0} {1} ".format(mean, sv)
    cmd += fastqs

    offset = int(opts.phred) if opts.phred else guessoffset([fastqfiles[0]])
    illumina = (offset == 64)
    if illumina:
        cmd += " -type illumina"
    if outtie:
        cmd += " -outtie"

    sh(cmd, outfile=frgfile)
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
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
604051
def scrobble_data_dir(dataDir, sampleMap, outF, qualF=None, idopt=None,
                      utf16=False):
    
    seqcount = 0
    outfiles = [osp.split(outF.name)[1]]
    if qualF:
        outfiles.append(osp.split(qualF.name)[1])

    for item in os.listdir(dataDir):
        if item in outfiles or not osp.isfile(os.path.join(dataDir, item)):
            continue
        # FASTA files
        if osp.splitext(item)[1] in file_types['fasta']:
            fh = open_enc(os.path.join(dataDir, item), utf16)
            records = SeqIO.parse(fh, 'fasta')
            for record in records:
                if isinstance(idopt, tuple):
                    sep, field = idopt
                    sampleID = record.id.split(sep)[field - 1]
                else:
                    sampleID = osp.splitext(item)[0]
                record.seq = (sampleMap[sampleID].barcode +
                              sampleMap[sampleID].primer +
                              record.seq)
                SeqIO.write(record, outF, 'fasta')
                seqcount += 1
            fh.close()
        # QUAL files
        elif qualF and osp.splitext(item)[1] in file_types['qual']:
            fh = open_enc(os.path.join(dataDir, item), utf16)
            records = SeqIO.parse(fh, 'qual')
            for record in records:
                mi = sampleMap[sampleMap.keys()[0]]
                quals = [40 for _ in range(len(mi.barcode) + len(mi.primer))]
                record.letter_annotations['phred_quality'][0:0] = quals
                SeqIO.write(record, qualF, 'qual')
            fh.close()
    return seqcount
756761
def get_seqs_fasta(seqs, names, out_fa):
    """"""
    with open(out_fa, 'w') as fa_handle:
        for s, n in itertools.izip(seqs, names):
            print(">cx{1}-{0}\n{0}".format(s, n), file=fa_handle)
    return out_fa