%prog catread fastqfile1 fastqfile2

    Concatenate paired end reads into one. Useful for example to do single-end
    mapping and perform filtering on the whole read pair level.
198722
def catread(args):
    
    p = OptionParser(catread.__doc__)
    opts, args = p.parse_args(args)

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

    r1, r2 = args
    p1fp, p2fp = FastqPairedIterator(r1, r2)
    outfile = pairspf((r1, r2)) + ".cat.fastq"
    fw = must_open(outfile, "w")
    while True:
        a = list(islice(p1fp, 4))
        if not a:
            break
        atitle, aseq, _, aqual = a
        btitle, bseq, _, bqual = list(islice(p2fp, 4))
        print("\n".join((atitle.strip(), aseq.strip() + bseq.strip(), \
                                "+", aqual.strip() + bqual.strip())), file=fw)
420566
def shred(args):
    
    p = OptionParser(shred.__doc__)
    p.set_depth(depth=2)
    p.add_option("--readlen", default=1000, type="int",
            help="Desired length of the reads [default: %default]")
    p.add_option("--minctglen", default=0, type="int",
            help="Ignore contig sequence less than [default: %default]")
    p.add_option("--shift", default=50, type="int",
            help="Overlap between reads must be at least [default: %default]")
    p.add_option("--fasta", default=False, action="store_true",
            help="Output shredded reads as FASTA sequences [default: %default]")
    opts, args = p.parse_args(args)

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

    fastafile, = args
    libID = fastafile.split(".")[0]
    depth = opts.depth
    readlen = opts.readlen
    shift = opts.shift

    outfile = libID + ".depth{0}".format(depth)
    if opts.fasta:
        outfile += ".fasta"
    else:
        outfile += ".frg"
    f = Fasta(fastafile, lazy=True)

    fw = must_open(outfile, "w", checkexists=True)
    if not opts.fasta:
       print(headerTemplate.format(libID=libID), file=fw)

    """
    Taken from runCA:

                    |*********|
                    |###################|
    |--------------------------------------------------|
     ---------------1---------------
               ---------------2---------------
                         ---------------3---------------
    *** - center_increments
    ### - center_range_width
    """
    for ctgID, (name, rec) in enumerate(f.iteritems_ordered()):
        seq = rec.seq
        seqlen = len(seq)
        if seqlen < opts.minctglen:
            continue

        shredlen = min(seqlen - shift, readlen)
        numreads = max(seqlen * depth / shredlen, 1)
        center_range_width = seqlen - shredlen

        ranges = []
        if depth == 1:
            if seqlen < readlen:
                ranges.append((0, seqlen))
            else:
                for begin in xrange(0, seqlen, readlen - shift):
                    end = min(seqlen, begin + readlen)
                    ranges.append((begin, end))
        else:
            if numreads == 1:
                ranges.append((0, shredlen))
            else:
                prev_begin = -1
                center_increments = center_range_width * 1. / (numreads - 1)
                for i in xrange(numreads):
                    begin = center_increments * i
                    end = begin + shredlen
                    begin, end = int(begin), int(end)

                    if begin == prev_begin:
                        continue

                    ranges.append((begin, end))
                    prev_begin = begin

        for shredID, (begin, end) in enumerate(ranges):
            shredded_seq = seq[begin:end]
            fragID = "{0}.{1}.frag{2}.{3}-{4}".format(libID, ctgID, shredID, begin, end)
            emitFragment(fw, fragID, libID, shredded_seq, fasta=opts.fasta)

    fw.close()
    logging.debug("Shredded reads are written to `{0}`.".format(outfile))
    return outfile
314860
def catFiles(filesToCat, catFile):
    
    if len(filesToCat) == 0: #We must handle this case or the cat call will hang waiting for input
        open(catFile, 'w').close()
        return
    maxCat = 25
    system("cat %s > %s" % (" ".join(filesToCat[:maxCat]), catFile))
    filesToCat = filesToCat[maxCat:]
    while len(filesToCat) > 0:
        system("cat %s >> %s" % (" ".join(filesToCat[:maxCat]), catFile))
        filesToCat = filesToCat[maxCat:]
1239141
def merge_pairs(self):
        
        logging.info('Merging paired reads')
        for sample in self.metadata:
            # Can only merge paired-end
            if len(sample.general.fastqfiles) == 2:
                # Set the name of the merged, and unmerged files
                sample.general.mergedreads = \
                    os.path.join(sample.general.outputdirectory, '{}_paired.fastq.gz'.format(sample.name))
                sample.general.unmergedforward = \
                    os.path.join(sample.general.outputdirectory, '{}_unpaired_R1.fastq.gz'.format(sample.name))
                sample.general.unmergedreverse = \
                    os.path.join(sample.general.outputdirectory, '{}_unpaired_R2.fastq.gz'.format(sample.name))
                try:
                    # Run the merging command - forward_in=sample.general.normalisedreads[0],
                    out, err, cmd = bbtools.bbmerge(forward_in=sorted(sample.general.trimmedcorrectedfastqfiles)[0],
                                                    merged_reads=sample.general.mergedreads,
                                                    returncmd=True,
                                                    outu1=sample.general.unmergedforward,
                                                    outu2=sample.general.unmergedreverse,
                                                    threads=self.cpus)
                    sample[self.analysistype].bbmergecmd = cmd
                    write_to_logfile(out, err, self.logfile, sample.general.logout, sample.general.logerr, None, None)
                except CalledProcessError:
                    delattr(sample.general, 'mergedreads')
                    delattr(sample.general, 'unmergedforward')
                    delattr(sample.general, 'unmergedreverse')
                except IndexError:
                    delattr(sample.general, 'mergedreads')
                    delattr(sample.general, 'unmergedforward')
                    delattr(sample.general, 'unmergedreverse')
            else:
                sample.general.mergedreads = sorted(sample.general.trimmedcorrectedfastqfiles)[0]
909489
def concat_main(args):
    

    # Get files directly
    if args.input_filepaths is not None:
        files = args.input_filepaths

    # Or find them
    else:
        files = get_file_list(args.file_wildcard)
        
        # No files found
        if len(files) == 0:
            msg = "No files were found. args.file_wildcard: {}".format(args.file_wildcard)
            logger.error(msg)
            raise Exception(msg)

    # Only 1 file found
    if len(files) == 1:
        logger.warning("Only 1 file found. No concatenation needs to be done, exiting")
        return

    # More than 1 file found
    else:
        # Parse each file and append to a list
        gctoos = []
        for f in files:
            gctoos.append(parse.parse(f))

        # Create concatenated gctoo object
        if args.concat_direction == "horiz":
            out_gctoo = hstack(gctoos, args.remove_all_metadata_fields, args.error_report_output_file,
                               args.fields_to_remove, args.reset_ids)

        elif args.concat_direction == "vert":
            out_gctoo = vstack(gctoos, args.remove_all_metadata_fields, args.error_report_output_file,
                               args.fields_to_remove, args.reset_ids)

    # Write out_gctoo to file
    logger.info("Writing to output file args.out_name:  {}".format(args.out_name))

    if args.out_type == "gctx":
        write_gctx.write(out_gctoo, args.out_name)

    elif args.out_type == "gct":
        write_gct.write(out_gctoo, args.out_name,
                          filler_null=args.filler_null,
                          metadata_null=args.metadata_null,
                          data_null=args.data_null)
1085785
def ispcr(args):
    
    from jcvi.utils.iter import grouper

    p = OptionParser(ispcr.__doc__)
    p.add_option("-r", dest="rclip", default=1, type="int",
            help="pair ID is derived from rstrip N chars [default: %default]")
    opts, args = p.parse_args(args)

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

    fastafile, = args
    ispcrfile = fastafile + ".isPcr"
    fw = open(ispcrfile, "w")

    N = opts.rclip
    strip_name = lambda x: x[:-N] if N else str

    npairs = 0
    fastaiter = SeqIO.parse(fastafile, "fasta")
    for a, b in grouper(fastaiter, 2):

        aid, bid = [strip_name(x) for x in (a.id, b.id)]
        assert aid == bid, "Name mismatch {0}".format((aid, bid))

        print("\t".join((aid, str(a.seq), str(b.seq))), file=fw)
        npairs += 1

    fw.close()
    logging.debug("A total of {0} pairs written to `{1}`.".\
                  format(npairs, ispcrfile))
1363440
def _make_reads_for_assembly(number_of_wanted_reads, total_reads, reads_in1, reads_in2, reads_out1, reads_out2, random_seed=None):
        ''''''
        random.seed(random_seed)

        if number_of_wanted_reads < total_reads:
            reads_written = 0
            percent_wanted = 100 * number_of_wanted_reads / total_reads
            file_reader1 = pyfastaq.sequences.file_reader(reads_in1)
            file_reader2 = pyfastaq.sequences.file_reader(reads_in2)
            out1 = pyfastaq.utils.open_file_write(reads_out1)
            out2 = pyfastaq.utils.open_file_write(reads_out2)

            for read1 in file_reader1:
                try:
                    read2 = next(file_reader2)
                except StopIteration:
                    pyfastaq.utils.close(out1)
                    pyfastaq.utils.close(out2)
                    raise Error('Error subsetting reads. No mate found for read ' + read1.id)

                if random.randint(0, 100) <= percent_wanted:
                    print(read1, file=out1)
                    print(read2, file=out2)
                    reads_written += 2

            pyfastaq.utils.close(out1)
            pyfastaq.utils.close(out2)
            return reads_written
        else:
            os.symlink(reads_in1, reads_out1)
            os.symlink(reads_in2, reads_out2)
            return total_reads
463453
def read_fastq(filename):
    
    if not filename:
        return itertools.cycle((None,))
    if filename == "-":
        filename_fh = sys.stdin
    elif filename.endswith('gz'):
        if is_python3:
            filename_fh = gzip.open(filename, mode='rt')
        else:
            filename_fh = BufferedReader(gzip.open(filename, mode='rt'))
    else:
        filename_fh = open(filename)
    return stream_fastq(filename_fh)
896104
def merge(args):
    
    p = OptionParser(merge.__doc__)
    p.set_outfile(outfile="merged_results.delta")
    opts, args = p.parse_args(args)

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

    ref, query = args[:2]
    deltafiles = args[2:]
    outfile = opts.outfile

    ref = get_abs_path(ref)
    query = get_abs_path(query)
    fw = must_open(outfile, "w")
    print(" ".join((ref, query)), file=fw)
    print("NUCMER", file=fw)
    fw.close()

    for d in deltafiles:
        cmd = "awk 'NR > 2 {{print $0}}' {0}".format(d)
        sh(cmd, outfile=outfile, append=True)
771506
def __read_single_fasta_query_lines(f):
    '''
    
    '''

    def readline():
        l = f.readline()
        if l == '':
            raise EOFError()
        return l

    rec = None
    try:
        l = readline()
        assert l.startswith('>')

        rec = [l]
        while True:
            pos = f.tell()
            l = readline()
            if l.startswith('>'):
                f.seek(pos, 0)
                break
            rec += [l]
    except EOFError:
        pass

    return rec
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])
944603
def cat_files(infiles, outfile):
    ''''''
    f_out = pyfastaq.utils.open_file_write(outfile)

    for filename in infiles:
        if os.path.exists(filename):
            f_in = pyfastaq.utils.open_file_read(filename)
            for line in f_in:
                print(line, end='', file=f_out)
            pyfastaq.utils.close(f_in)

    pyfastaq.utils.close(f_out)
881969
def overlapbatch(args):
    
    p = OptionParser(overlap.__doc__)
    opts, args = p.parse_args(args)
    if len(args) != 2:
        sys.exit(not p.print_help())

    ctgfasta, poolfasta = args
    f = Fasta(ctgfasta)
    for k, rec in f.iteritems_ordered():
        fastafile = k + ".fasta"
        fw = open(fastafile, "w")
        SeqIO.write([rec], fw, "fasta")
        fw.close()

        overlap([fastafile, poolfasta])
199809
def shuffle(args):
    
    p = OptionParser(shuffle.__doc__)
    p.set_tag()
    opts, args = p.parse_args(args)

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

    p1, p2 = args
    pairsfastq = pairspf((p1, p2)) + ".fastq"
    tag = opts.tag

    p1fp = must_open(p1)
    p2fp = must_open(p2)
    pairsfw = must_open(pairsfastq, "w")
    nreads = 0
    while True:
        a = list(islice(p1fp, 4))
        if not a:
            break

        b = list(islice(p2fp, 4))
        if tag:
            name = a[0].rstrip()
            a[0] = name + "/1\n"
            b[0] = name + "/2\n"

        pairsfw.writelines(a)
        pairsfw.writelines(b)
        nreads += 2

    pairsfw.close()
    extra = nreads * 2 if tag else 0
    checkShuffleSizes(p1, p2, pairsfastq, extra=extra)

    logging.debug("File `{0}` verified after writing {1} reads.".\
                     format(pairsfastq, nreads))
    return pairsfastq
371787
def merge_fastqs(job, job_vars):
    
    input_args, ids = job_vars
    work_dir = job.fileStore.getLocalTempDir()
    cores = input_args['cpu_count']
    single_end_reads = input_args['single_end_reads']
    # I/O
    sample = return_input_paths(job, work_dir, ids, 'sample.tar')
    # Untar File
    # subprocess.check_call(['unzip', sample, '-d', work_dir])
    subprocess.check_call(['tar', '-xvf', sample, '-C', work_dir])
    # Remove large files before creating concat versions.
    os.remove(os.path.join(work_dir, 'sample.tar'))
    # Zcat files in parallel
    if single_end_reads:
        files = sorted(glob.glob(os.path.join(work_dir, '*')))
        with open(os.path.join(work_dir, 'R1.fastq'), 'w') as f1:
            subprocess.check_call(['zcat'] + files, stdout=f1)
        # FileStore
        ids['R1.fastq'] = job.fileStore.writeGlobalFile(os.path.join(work_dir, 'R1.fastq'))
    else:
        r1_files = sorted(glob.glob(os.path.join(work_dir, '*R1*')))
        r2_files = sorted(glob.glob(os.path.join(work_dir, '*R2*')))
        with open(os.path.join(work_dir, 'R1.fastq'), 'w') as f1:
            p1 = subprocess.Popen(['zcat'] + r1_files, stdout=f1)
        with open(os.path.join(work_dir, 'R2.fastq'), 'w') as f2:
            p2 = subprocess.Popen(['zcat'] + r2_files, stdout=f2)
        p1.wait()
        p2.wait()
        # FileStore
        ids['R1.fastq'] = job.fileStore.writeGlobalFile(os.path.join(work_dir, 'R1.fastq'))
        ids['R2.fastq'] = job.fileStore.writeGlobalFile(os.path.join(work_dir, 'R2.fastq'))
    job.fileStore.deleteGlobalFile(ids['sample.tar'])
    # Spawn child job
    return job.addChildJobFn(mapsplice, job_vars, cores=cores, disk='130 G').rv()
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
216287
def _all_reads_from_contig(self, contig, fout):
        ''''''
        sam_reader = pysam.Samfile(self.bam, "rb")
        for read in sam_reader.fetch(contig):
            print(mapping.aligned_read_to_read(read, ignore_quality=not self.fastq_out), file=fout)