%prog fastq fastqfile

    Convert reads formatted as FASTQ file, and convert to CA frg file.
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)
264762
def fastqc(self, file, output_dir):
        
        # You can find the fastqc help with fastqc --help
        try:
            pm = self.pm
        except AttributeError:
            # Do nothing, this is just for path construction.
            pass
        else:
            if not os.path.isabs(output_dir) and pm is not None:
                output_dir = os.path.join(pm.outfolder, output_dir)
        self.make_sure_path_exists(output_dir)
        return "{} --noextract --outdir {} {}".\
                format(self.tools.fastqc, output_dir, file)
122892
def dump(args):
    
    p = OptionParser(dump.__doc__)
    p.add_option("--dir",
                help="Working directory [default: %default]")
    p.add_option("--nosim", default=False, action="store_true",
                 help="Do not simulate qual to 50 [default: %default]")
    opts, args = p.parse_args(args)

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

    fastbfile, = args
    d = opts.dir
    if d:
        from jcvi.assembly.preprocess import export_fastq

        rc = "jump" in fastbfile
        export_fastq(d, fastbfile, rc=rc)
        return

    sim = not opts.nosim
    pf = "j" if "jump" in fastbfile else "f"

    statsfile = "{0}.lib_stats".format(pf)
    if op.exists(statsfile):
        os.remove(statsfile)

    cmd = "SplitReadsByLibrary READS_IN={0}".format(fastbfile)
    cmd += " READS_OUT={0} QUALS=True".format(pf)
    sh(cmd)

    libs = []
    fp = open(statsfile)
    next(fp); next(fp)  # skip two rows
    for row in fp:
        if row.strip() == "":
            continue

        libname = row.split()[0]
        if libname == "Unpaired":
            continue

        libs.append(libname)

    logging.debug("Found libraries: {0}".format(",".join(libs)))

    cmds = []
    for libname in libs:
        cmd = "FastbQualbToFastq"
        cmd += " HEAD_IN={0}.{1}.AB HEAD_OUT={1}".format(pf, libname)
        cmd += " PAIRED=True PHRED_OFFSET=33"
        if sim:
            cmd += " SIMULATE_QUALS=True"
        if pf == 'j':
            cmd += " FLIP=True"

        cmds.append((cmd, ))

    m = Jobs(target=sh, args=cmds)
    m.run()

    for libname in libs:
        cmd = "mv {0}.A.fastq {0}.1.fastq".format(libname)
        sh(cmd)
        cmd = "mv {0}.B.fastq {0}.2.fastq".format(libname)
        sh(cmd)
1058553
def process_fastq_minimal(fastq, **kwargs):
    """"""
    infastq = handle_compressed_input(fastq)
    try:
        df = pd.DataFrame(
            data=[rec for rec in fq_minimal(infastq) if rec],
            columns=["timestamp", "lengths"]
        )
    except IndexError:
        logging.error("Fatal: Incorrect file structure for fastq_minimal")
        sys.exit("Error: file does not match expected structure for fastq_minimal")
    return ut.reduce_memory_usage(df)
526958
def some(args):
    
    p = OptionParser(some.__doc__)
    p.add_option("--exclude", default=False, action="store_true",
            help="Output sequences not in the list file [default: %default]")
    p.add_option("--uniprot", default=False, action="store_true",
            help="Header is from uniprot [default: %default]")

    opts, args = p.parse_args(args)

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

    fastafile, listfile, outfastafile = args
    outfastahandle = must_open(outfastafile, "w")
    qualfile = get_qual(fastafile)

    names = set(x.strip() for x in open(listfile))
    if qualfile:
        outqualfile = outfastafile + ".qual"
        outqualhandle = open(outqualfile, "w")
        parser = iter_fasta_qual(fastafile, qualfile)
    else:
        parser = SeqIO.parse(fastafile, "fasta")

    num_records = 0
    for rec in parser:
        name = rec.id
        if opts.uniprot:
            name = name.split("|")[-1]

        if opts.exclude:
            if name in names:
                continue
        else:
            if name not in names:
                continue

        SeqIO.write([rec], outfastahandle, "fasta")
        if qualfile:
            SeqIO.write([rec], outqualhandle, "qual")

        num_records += 1

    logging.debug("A total of %d records written to `%s`" % \
            (num_records, outfastafile))
875469
def clr(args):
    
    p = OptionParser(clr.__doc__)
    opts, args = p.parse_args(args)

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

    blastfile = args[0]
    fastafiles = args[1:]

    sizes = {}
    for fa in fastafiles:
        f = Fasta(fa)
        sizes.update(f.itersizes())

    b = Blast(blastfile)
    for query, hits in b.iter_hits():

        qsize = sizes[query]
        vectors = list((x.qstart, x.qstop) for x in hits)
        vmin, vmax = range_minmax(vectors)

        left_size = vmin - 1
        right_size = qsize - vmax

        if left_size > right_size:
            clr_start, clr_end = 0, vmin
        else:
            clr_start, clr_end = vmax, qsize

        print("\t".join(str(x) for x in (query, clr_start, clr_end)))
        del sizes[query]

    for q, size in sorted(sizes.items()):
        print("\t".join(str(x) for x in (q, 0, size)))
1286087
def is_transformed(fastq):
    

    with open_fastq(fastq) as in_handle:
        for line in islice(in_handle, 400):
            if "UMI_" in line:
                return True
    return False
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
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()
644764
def suffix(args):
    
    p = OptionParser(suffix.__doc__)
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    fastqfile, sf = args
    fw = must_open(opts.outfile, "w")
    nreads = nselected = 0
    for rec in iter_fastq(fastqfile):
        nreads += 1
        if rec is None:
            break
        if rec.seq.endswith(sf):
            print(rec, file=fw)
            nselected += 1
    logging.debug("Selected reads with suffix {0}: {1}".\
                  format(sf, percentage(nselected, nreads)))
46355
def fpkm(args):
    
    p = OptionParser(fpkm.__doc__)
    opts, args = p.parse_args(args)

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

    fastafile = args[0]
    bamfiles = args[1:]
    # Create a DUMMY gff file for cuffdiff
    gffile = fastafile.rsplit(".", 1)[0] + ".gff"
    if need_update(fastafile, gffile):
        fw = open(gffile, "w")
        f = Fasta(fastafile, lazy=True)
        for key, size in f.itersizes_ordered():
            print("\t".join(str(x) for x in (key, "dummy", "transcript",\
                1, size, ".", ".", ".", "ID=" + key)), file=fw)
        fw.close()
        logging.debug("Dummy GFF created: {0}".format(gffile))

    cmd = "cuffdiff {0} {1}".format(gffile, " ".join(bamfiles))
    sh(cmd)
242949
def to_fastq_str(self):
    
    return "@" + self.name + "\n" + self.sequenceData +\
           "\n" + "+" + self.name + "\n" + self.seq_qual
70351
def fromagp(args):
    
    from jcvi.formats.agp import AGP
    from jcvi.formats.sizes import Sizes

    p = OptionParser(fromagp.__doc__)
    p.add_option("--novalidate", default=False, action="store_true",
                 help="Do not validate AGP")
    opts, args = p.parse_args(args)

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

    agpfile, componentfasta, objectfasta = args
    chainfile = agpfile.rsplit(".", 1)[0] + ".chain"
    fw = open(chainfile, "w")
    agp = AGP(agpfile, validate=(not opts.novalidate))
    componentsizes = Sizes(componentfasta).mapping
    objectsizes = Sizes(objectfasta).mapping
    chain = "chain"
    score = 1000
    tStrand = "+"
    id = 0
    for a in agp:
        if a.is_gap:
            continue

        tName = a.component_id
        tSize = componentsizes[tName]
        tStart = a.component_beg
        tEnd = a.component_end
        tStart -= 1

        qName = a.object
        qSize = objectsizes[qName]
        qStrand = "-" if a.orientation == "-" else "+"
        qStart = a.object_beg
        qEnd = a.object_end
        if qStrand == '-':
            _qStart = qSize - qEnd + 1
            _qEnd = qSize - qStart + 1
            qStart, qEnd = _qStart, _qEnd
        qStart -= 1

        id += 1
        size = a.object_span
        headerline = "\t".join(str(x) for x in (
             chain, score, tName, tSize, tStrand, tStart,
             tEnd, qName, qSize, qStrand, qStart, qEnd, id
        ))
        alignmentline = size
        print(headerline, file=fw)
        print(alignmentline, file=fw)
        print(file=fw)

    fw.close()
    logging.debug("File written to `{0}`.".format(chainfile))
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())
255055
def handle_program_options():
    
    parser = argparse.ArgumentParser(description="Convert Sanger-sequencing \
                                     derived data files for use with the \
                                     metagenomics analysis program QIIME, by \
                                     extracting Sample ID information, adding\
                                     barcodes and primers to the sequence \
                                     data, and outputting a mapping file and\
                                     single FASTA-formatted sequence file \
                                     formed by concatenating all input data.")
    parser.add_argument('-i', '--input_dir', required=True,
                        help="The directory containing sequence data files. \
                              Assumes all data files are placed in this \
                              directory. For files organized within folders by\
                              sample, use -s in addition.")
    parser.add_argument('-m', '--map_file', default='map.txt',
                        help="QIIME-formatted mapping file linking Sample IDs \
                              with barcodes and primers.")
    parser.add_argument('-o', '--output', default='output.fasta',
                        metavar='OUTPUT_FILE',
                        help="Single file containing all sequence data found \
                              in input_dir, FASTA-formatted with barcode and \
                              primer preprended to sequence. If the -q option \
                              is passed, any quality data will also be output \
                              to a single file of the same name with a .qual \
                              extension.")
    parser.add_argument('-b', '--barcode_length', type=int, default=12,
                        help="Length of the generated barcode sequences. \
                              Default is 12 (QIIME default), minimum is 8.")

    parser.add_argument('-q', '--qual', action='store_true', default=False,
                        help="Instruct the program to look for quality \
                              input files")
    parser.add_argument('-u', '--utf16', action='store_true', default=False,
                        help="UTF-16 encoded input files")

    parser.add_argument('-t', '--treatment',
                        help="Inserts an additional column into the mapping \
                              file specifying some treatment or other variable\
                              that separates the current set of sequences \
                              from any other set of seqeunces. For example:\
                              -t DiseaseState=healthy")

    # data input options
    sidGroup = parser.add_mutually_exclusive_group(required=True)
    sidGroup.add_argument('-d', '--identifier_pattern',
                          action=ValidateIDPattern,
                          nargs=2, metavar=('SEPARATOR', 'FIELD_NUMBER'),
                          help="Indicates how to extract the Sample ID from \
                               the description line. Specify two things: \
                               1. Field separator, 2. Field number of Sample \
                               ID (1 or greater). If the separator is a space \
                               or tab, use \s or \\t respectively. \
                               Example: >ka-SampleID-2091, use -i - 2, \
                               indicating - is the separator and the Sample ID\
                               is field #2.")
    sidGroup.add_argument('-f', '--filename_sample_id', action='store_true',
                          default=False, help='Specify that the program should\
                          the name of each fasta file as the Sample ID for use\
                          in the mapping file. This is meant to be used when \
                          all sequence data for a sample is stored in a single\
                          file.')

    return parser.parse_args()
593965
def run_bcl2fastq(run_folder, ss_csv, config):
    
    bc_dir = os.path.join(run_folder, "Data", "Intensities", "BaseCalls")
    output_dir = os.path.join(run_folder, "fastq")

    if not os.path.exists(os.path.join(output_dir, "Makefile")):
        subprocess.check_call(["configureBclToFastq.pl", "--no-eamss",
                               "--input-dir", bc_dir, "--output-dir", output_dir,
                               "--sample-sheet", ss_csv])
    with utils.chdir(output_dir):
        cores = str(utils.get_in(config, ("algorithm", "num_cores"), 1))
        cmd = ["make", "-j", cores]
        if "submit_cmd" in config["process"] and "bcl2fastq_batch" in config["process"]:
            _submit_and_wait(cmd, cores, config, output_dir)
        else:
            subprocess.check_call(cmd)
    return output_dir
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()