%prog needle nw.pairs a.pep.fasta b.pep.fasta

    Take protein pairs and needle them
    Automatically writes output file `nw.scores`
167965
def needle(args):
    
    from jcvi.formats.fasta import Fasta, SeqIO

    p = OptionParser(needle.__doc__)

    opts, args = p.parse_args(args)

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

    manager = mp.Manager()
    results = manager.list()
    needle_pool = mp.Pool(processes=mp.cpu_count())

    pairsfile, apep, bpep = args
    afasta, bfasta = Fasta(apep), Fasta(bpep)
    fp = must_open(pairsfile)
    for i, row in enumerate(fp):
        a, b = row.split()
        a, b = afasta[a], bfasta[b]
        fa, fb = must_open("{0}_{1}_a.fasta".format(pairsfile, i), "w"), \
            must_open("{0}_{1}_b.fasta".format(pairsfile, i), "w")
        SeqIO.write([a], fa, "fasta")
        SeqIO.write([b], fb, "fasta")
        fa.close()
        fb.close()

        needlefile = "{0}_{1}_ab.needle".format(pairsfile, i)
        needle_pool.apply_async(_needle, \
            (fa.name, fb.name, needlefile, a.id, b.id, results))

    needle_pool.close()
    needle_pool.join()

    fp.close()

    scoresfile = "{0}.scores".format(pairsfile.rsplit(".")[0])
    fw = must_open(scoresfile, "w")
    for result in results:
        print(result, file=fw)
    fw.close()
1341677
def write_seq_as_temp_fasta(seq):
    
    sr = ssbio.protein.sequence.utils.cast_to_seq_record(seq, id='tempfasta')
    return write_fasta_file(seq_records=sr, outname='temp', outdir=tempfile.gettempdir(), force_rerun=True)
58612
def writeFasta(sequence, sequence_name, output_file):
    
    i=0
    f=open(output_file,'w')
    f.write(">"+str(sequence_name)+"\n")
    while i <= len(sequence):
        f.write(sequence[i:i+60]+"\n")
        i=i+60
    f.close()
188179
def create_psm_lookup(fn, fastafn, mapfn, header, pgdb, unroll=False,
                      specfncol=None, decoy=False,
                      fastadelim=None, genefield=None):
    
    proteins = store_proteins_descriptions(pgdb, fastafn, fn, mapfn, header,
                                           decoy, fastadelim, genefield)
    mzmlmap = pgdb.get_mzmlfile_map()
    sequences = {}
    for psm in tsvreader.generate_tsv_psms(fn, header):
        seq = tsvreader.get_psm_sequence(psm, unroll)
        sequences[seq] = 1
    pgdb.store_pepseqs(((seq,) for seq in sequences))
    pepseqmap = pgdb.get_peptide_seq_map()
    psms = []
    for row, psm in enumerate(tsvreader.generate_tsv_psms(fn, header)):
        specfn, psm_id, scan, seq, score = tsvreader.get_psm(psm, unroll,
                                                             specfncol)
        if len(psms) % DB_STORE_CHUNK == 0:
            pgdb.store_psms(psms)
            psms = []
        psms.append({'rownr': row,
                     'psm_id': psm_id,
                     'seq': pepseqmap[seq],
                     'score': score,
                     'specfn': mzmlmap[specfn],
                     'scannr': scan,
                     'spec_id': '{}_{}'.format(mzmlmap[specfn], scan),
                     })
    pgdb.store_psms(psms)
    pgdb.index_psms()
    store_psm_protein_relations(fn, header, pgdb, proteins)
322514
def ld(args):
    
    import numpy as np
    from random import sample

    from jcvi.algorithms.matrix import symmetrize

    p = OptionParser(ld.__doc__)
    p.add_option("--subsample", default=1000, type="int",
                 help="Subsample markers to speed up [default: %default]")
    opts, args, iopts = p.set_image_options(args, figsize="8x8")

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

    mstmap, = args
    subsample = opts.subsample
    data = MSTMap(mstmap)

    markerbedfile = mstmap + ".subsample.bed"
    ldmatrix = mstmap + ".subsample.matrix"
    # Take random subsample while keeping marker order
    if subsample < data.nmarkers:
        data = [data[x] for x in \
                sorted(sample(xrange(len(data)), subsample))]
    else:
        logging.debug("Use all markers, --subsample ignored")

    nmarkers = len(data)
    if need_update(mstmap, (ldmatrix, markerbedfile)):
        fw = open(markerbedfile, "w")
        print("\n".join(x.bedline for x in data), file=fw)
        logging.debug("Write marker set of size {0} to file `{1}`."\
                        .format(nmarkers, markerbedfile))
        fw.close()

        M = np.zeros((nmarkers, nmarkers), dtype=float)
        for i, j in combinations(range(nmarkers), 2):
            a = data[i]
            b = data[j]
            M[i, j] = calc_ldscore(a.genotype, b.genotype)

        M = symmetrize(M)

        logging.debug("Write LD matrix to file `{0}`.".format(ldmatrix))
        M.tofile(ldmatrix)
    else:
        nmarkers = len(Bed(markerbedfile))
        M = np.fromfile(ldmatrix, dtype="float").reshape(nmarkers, nmarkers)
        logging.debug("LD matrix `{0}` exists ({1}x{1})."\
                        .format(ldmatrix, nmarkers))

    from jcvi.graphics.base import plt, savefig, Rectangle, draw_cmap

    plt.rcParams["axes.linewidth"] = 0

    fig = plt.figure(1, (iopts.w, iopts.h))
    root = fig.add_axes([0, 0, 1, 1])
    ax = fig.add_axes([.1, .1, .8, .8])  # the heatmap

    ax.matshow(M, cmap=iopts.cmap)

    # Plot chromosomes breaks
    bed = Bed(markerbedfile)
    xsize = len(bed)
    extent = (0, nmarkers)
    chr_labels = []
    ignore_size = 20

    for (seqid, beg, end) in bed.get_breaks():
        ignore = abs(end - beg) < ignore_size
        pos = (beg + end) / 2
        chr_labels.append((seqid, pos, ignore))
        if ignore:
            continue
        ax.plot((end, end), extent, "w-", lw=1)
        ax.plot(extent, (end, end), "w-", lw=1)

    # Plot chromosome labels
    for label, pos, ignore in chr_labels:
        pos = .1 + pos * .8 / xsize
        if not ignore:
            root.text(pos, .91, label,
                ha="center", va="bottom", rotation=45, color="grey")
            root.text(.09, pos, label,
                ha="right", va="center", color="grey")

    ax.set_xlim(extent)
    ax.set_ylim(extent)
    ax.set_axis_off()

    draw_cmap(root, "Pairwise LD (r2)", 0, 1, cmap=iopts.cmap)

    root.add_patch(Rectangle((.1, .1), .8, .8, fill=False, ec="k", lw=2))
    m = mstmap.split(".")[0]
    root.text(.5, .06, "Linkage Disequilibrium between {0} markers".format(m), ha="center")

    root.set_xlim(0, 1)
    root.set_ylim(0, 1)
    root.set_axis_off()

    image_name = m + ".subsample" + "." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)
697094
def needle_msa(reference, results, gap_open=-15, gap_extend=0,
               matrix=submat.DNA_SIMPLE):
    '''

    '''
    gap = '-'
    # Convert alignments to list of strings
    alignments = []
    for result in results:
        ref_dna, res_dna, score = needle(reference, result, gap_open=gap_open,
                                         gap_extend=gap_extend,
                                         matrix=matrix)
        alignments.append([str(ref_dna), str(res_dna), score])

    def insert_gap(sequence, position):
        return sequence[:position] + gap + sequence[position:]

    i = 0
    while True:
        # Iterate over 'columns' in every reference
        refs = [alignment[0][i] for alignment in alignments]

        # If there's a non-unanimous gap, insert gap into alignments
        gaps = [ref == gap for ref in refs]
        if any(gaps) and not all(gaps):
            for alignment in alignments:
                if alignment[0][i] != gap:
                    alignment[0] = insert_gap(alignment[0], i)
                    alignment[1] = insert_gap(alignment[1], i)

        # If all references match, we're all done
        alignment_set = set(alignment[0] for alignment in alignments)
        if len(alignment_set) == 1:
            break

        # If we've reach the end of some, but not all sequences, add end gap
        lens = [len(alignment[0]) for alignment in alignments]
        if i + 1 in lens:
            for alignment in alignments:
                if len(alignment[0]) == i + 1:
                    alignment[0] = alignment[0] + gap
                    alignment[1] = alignment[1] + gap

        i += 1

        if i > 20:
            break

    # Convert into MSA format
    output_alignment = [cr.DNA(alignments[0][0])]
    for alignment in alignments:
        output_alignment.append(cr.DNA(alignment[1]))

    return output_alignment
139759
def compile_output(self, input_sample, num_samples, num_params,
                       maximum_combo, num_groups=None):
        

        if num_groups is None:
            num_groups = num_params

        self.check_input_sample(input_sample, num_groups, num_samples)

        index_list = self._make_index_list(num_samples, num_params, num_groups)

        output = np.zeros(
            (np.size(maximum_combo) * (num_groups + 1), num_params))
        for counter, combo in enumerate(maximum_combo):
            output[index_list[counter]] = np.array(
                input_sample[index_list[combo]])
        return output
954950
def generate_proteins(pepfn, proteins, pepheader, scorecol, minlog,
                      higherbetter=True, protcol=False):
    """"""
    protein_peptides = {}
    if minlog:
        higherbetter = False
    if not protcol:
        protcol = peptabledata.HEADER_MASTERPROTEINS
    for psm in reader.generate_tsv_psms(pepfn, pepheader):
        p_acc = psm[protcol]
        if ';' in p_acc:
            continue
        protein_peptides = evaluate_peptide(protein_peptides, psm, p_acc,
                                            higherbetter, scorecol,
                                            fncol=False)
    if minlog:
        try:
            nextbestscore = min([pep['score'] for pep in
                                 protein_peptides.values()
                                 if pep['score'] > 0])
        except ValueError:
            import sys
            sys.stderr.write('Cannot find score of type {} which is above 0. '
                             'Only scores above zero can have a -log value. '
                             'Exiting.'.format(scorecol))
            sys.exit(1)
        nextbestscore = -log(nextbestscore, 10)
    for protein in proteins:
        try:
            peptide = protein_peptides[protein[prottabledata.HEADER_PROTEIN]]
        except KeyError:
            print('WARNING - protein {} not found in peptide '
                  'table'.format(protein[prottabledata.HEADER_PROTEIN]))
            peptide = {'score': 'NA'}
        if minlog and peptide['score'] != 'NA':
            peptide['score'] = log_score(peptide['score'], nextbestscore)
        protein[prottabledata.HEADER_QSCORE] = str(
            peptide['score'])
        yield protein
1082530
def write(self):
        

        rabs = [x.seqrecord for x in self.gproteins.values() if x.is_rab()]
        return SeqIO.write(rabs, self.tmpfname + '.phase2', 'fasta')
471534
def fromimpute2(args):
    
    p = OptionParser(fromimpute2.__doc__)
    opts, args = p.parse_args(args)

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

    impute2file, fastafile, chr = args
    fasta = Fasta(fastafile)
    print(get_vcfstanza(fastafile, fasta))
    fp = open(impute2file)
    seen = set()
    for row in fp:
        snp_id, rsid, pos, ref, alt, aa, ab, bb = row.split()
        pos = int(pos)
        if pos in seen:
            continue
        seen.add(pos)
        code = max((float(aa), "0/0"), (float(ab), "0/1"), (float(bb), "1/1"))[-1]
        tag = "PR" if snp_id == chr else "IM"
        print("\t".join(str(x) for x in \
                (chr, pos, rsid, ref, alt, ".", ".", tag, \
                "GT:GP", code + ":" + ",".join((aa, ab, bb)))))
669003
def print_genome_matrix(hits, fastas, id2desc, file_name):
    
    out = open(file_name, 'w')
    fastas = sorted(fastas)
    print('## percent identity between genomes', file=out)
    print('# - \t %s' % ('\t'.join(fastas)), file=out)
    for fasta in fastas:
        line = [fasta]
        for other in fastas:
            if other == fasta:
                average = '-'
            else:
                average = numpy.average([hits[fasta][other][i][3] for i in hits[fasta][other]])
            line.append(str(average))
        print('\t'.join(line), file=out)
    print('', file=out)
    print('## percent of orfs that are orthologous between genomes', file=out)
    print('# - \t %s' % ('\t'.join(fastas)), file=out)
    for fasta in fastas:
        line = [fasta]
        for other in fastas:
            if other == fasta:
                percent = '-'
            else:
                orthologs = float(len(hits[fasta][other]))
                orfs = float(len([i for i in id2desc if id2desc[i][0] == fasta]))
                percent = float(orthologs / orfs) * 100
            line.append(str(percent))
        print('\t'.join(line), file=out)
888641
def pasteprepare(args):
    
    p = OptionParser(pasteprepare.__doc__)
    p.add_option("--flank", default=5000, type="int",
                 help="Get the seq of size on two ends [default: %default]")
    opts, args = p.parse_args(args)

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

    goodfasta, = args
    flank = opts.flank
    pf = goodfasta.rsplit(".", 1)[0]
    extbed = pf + ".ext.bed"

    sizes = Sizes(goodfasta)
    fw = open(extbed, "w")
    for bac, size in sizes.iter_sizes():
        print("\t".join(str(x) for x in \
                               (bac, 0, min(flank, size), bac + "L")), file=fw)
        print("\t".join(str(x) for x in \
                               (bac, max(size - flank, 0), size, bac + "R")), file=fw)
    fw.close()

    fastaFromBed(extbed, goodfasta, name=True)
61893
def main(arguments):
    
    global verbose
    global veryVerbose
    global iteration_num
    global single_score
    global pr_flag
    global match_triple_dict
    # set the iteration number
    # total iteration number = restart number + 1
    iteration_num = arguments.r + 1
    if arguments.ms:
        single_score = False
    if arguments.v:
        verbose = True
    if arguments.vv:
        veryVerbose = True
    if arguments.pr:
        pr_flag = True
    # significant digits to print out
    floatdisplay = "%%.%df" % arguments.significant
    for (precision, recall, best_f_score) in score_amr_pairs(args.f[0], args.f[1],
                                                             justinstance=arguments.justinstance,
                                                             justattribute=arguments.justattribute,
                                                             justrelation=arguments.justrelation):
        # print("Sentence", sent_num)
        if pr_flag:
            print("Precision: " + floatdisplay % precision)
            print("Recall: " + floatdisplay % recall)
        print("F-score: " + floatdisplay % best_f_score)
    args.f[0].close()
    args.f[1].close()
193971
def _needle(fa, fb, needlefile, a, b, results):
    
    from Bio.Emboss.Applications import NeedleCommandline

    needle_cline = NeedleCommandline(asequence=fa, bsequence=fb,
                        gapopen=10, gapextend=0.5, outfile=needlefile)
    stdout, stderr = needle_cline()
    nh = NeedleHeader(needlefile)
    FileShredder([fa, fb, needlefile], verbose=False)
    r = ["\t".join((a, b, nh.identity, nh.score))]

    results.extend(r)
921743
def build_parser(parser):
    
    parser.add_argument(
        'sequence_file',
        type=FileType('r'),
        help="""Input fastq file. A fasta-format file may also be provided
            if --input-qual is also specified.""")
    parser.add_argument(
        '--input-qual',
        type=FileType('r'),
        help="""The quality scores associated with the input file. Only
            used if input file is fasta.""")
    parser.add_argument(
        'output_file',
        type=FileType('w'),
        help="""Output file. Format determined from extension.""")

    output_group = parser.add_argument_group("Output")
    output_group.add_argument(
        '--report-out',
        type=FileType('w'),
        default=sys.stdout,
        help="""Output file for report [default:
            stdout]""")
    output_group.add_argument(
        '--details-out',
        type=FileType('w'),
        help="""Output file to report fate of each sequence""")
    output_group.add_argument(
        '--no-details-comment',
        action='store_false',
        default=True,
        dest='details_comment',
        help="""Do not write comment
            lines with version and call to start --details-out""")

    parser.add_argument(
        '--min-mean-quality',
        metavar='QUALITY',
        type=float,
        default=DEFAULT_MEAN_SCORE,
        help="""Minimum mean quality score for
            each read [default: %(default)s]""")
    parser.add_argument(
        '--min-length',
        metavar='LENGTH',
        type=int,
        default=200,
        help="""Minimum length to keep sequence [default:
            %(default)s]""")
    parser.add_argument(
        '--max-length',
        metavar='LENGTH',
        type=int,
        default=1000,
        help="""Maximum length to keep before truncating
            [default: %(default)s]. This operation occurs before
            --max-ambiguous""")

    window_group = parser.add_argument_group('Quality window options')
    window_group.add_argument(
        '--quality-window-mean-qual',
        type=float,
        help="""Minimum quality score within the window defined by
            --quality-window. [default: same as --min-mean-quality]""")
    window_group.add_argument(
        '--quality-window-prop',
        help="""Proportion of
            reads within quality window to that must pass filter. Floats are [default:
            %(default).1f]""",
        default=1.0,
        type=typed_range(float, 0.0, 1.0))
    window_group.add_argument(
        '--quality-window',
        type=int,
        metavar='WINDOW_SIZE',
        default=0,
        help="""Window size for truncating sequences.  When set
            to a non-zero value, sequences are truncated where the mean mean
            quality within the window drops below --min-mean-quality.
            [default: %(default)s]""")

    parser.add_argument(
        '--ambiguous-action',
        choices=('truncate', 'drop'),
        help="""Action to take on ambiguous base in sequence (N's).
            [default: no action]""")
    parser.add_argument(
        '--max-ambiguous',
        default=None,
        help="""Maximum number
            of ambiguous bases in a sequence. Sequences exceeding this count
            will be removed.""",
        type=int)
    parser.add_argument(
        '--pct-ambiguous',
        help="""Maximun percent of
            ambiguous bases in a sequence.  Sequences exceeding this percent
            will be removed.""",
        type=float)

    barcode_group = parser.add_argument_group('Barcode/Primer')
    primer_group = barcode_group.add_mutually_exclusive_group()
    primer_group.add_argument(
        '--primer', help="""IUPAC ambiguous primer to
            require""")
    primer_group.add_argument(
        '--no-primer',
        help="""Do not use a primer.""",
        action='store_const',
        const='',
        dest='primer')
    barcode_group.add_argument(
        '--barcode-file',
        help="""CSV file containing
            sample_id,barcode[,primer] in the rows. A single primer for all
            sequences may be specified with `--primer`, or `--no-primer` may be
            used to indicate barcodes should be used without a primer
            check.""",
        type=FileType('r'))
    barcode_group.add_argument(
        '--barcode-header',
        action='store_true',
        default=False,
        help="""Barcodes have a header row [default:
            %(default)s]""")
    barcode_group.add_argument(
        '--map-out',
        help="""Path to write
            sequence_id,sample_id pairs""",
        type=FileType('w'),
        metavar='SAMPLE_MAP')
    barcode_group.add_argument(
        '--quoting',
        help="""A string naming an
            attribute of the csv module defining the quoting behavior for
            `SAMPLE_MAP`.  [default: %(default)s]""",
        default='QUOTE_MINIMAL',
        choices=[s for s in dir(csv) if s.startswith('QUOTE_')])
580919
def report_pairs(data, cutoff=0, mateorientation=None,
        pairsfile=None, insertsfile=None, rclip=1, ascii=False, bins=20,
        distmode="ss", mpcutoff=1000):
    
    allowed_mateorientations = ("++", "--", "+-", "-+")

    if mateorientation:
        assert mateorientation in allowed_mateorientations

    num_fragments, num_pairs = 0, 0

    all_dist = []
    linked_dist = []
    # +- (forward-backward) is `innie`, -+ (backward-forward) is `outie`
    orientations = defaultdict(int)

    # clip how many chars from end of the read name to get pair name
    key = (lambda x: x.accn[:-rclip]) if rclip else (lambda x: x.accn)
    data.sort(key=key)

    if pairsfile:
        pairsfw = open(pairsfile, "w")
    if insertsfile:
        insertsfw = open(insertsfile, "w")

    for pe, lines in groupby(data, key=key):
        lines = list(lines)
        if len(lines) != 2:
            num_fragments += len(lines)
            continue

        num_pairs += 1
        a, b = lines

        asubject, astart, astop = a.seqid, a.start, a.end
        bsubject, bstart, bstop = b.seqid, b.start, b.end

        aquery, bquery = a.accn, b.accn
        astrand, bstrand = a.strand, b.strand

        dist, orientation = range_distance(\
                (asubject, astart, astop, astrand),
                (bsubject, bstart, bstop, bstrand),
                distmode=distmode)

        if dist >= 0:
            all_dist.append((dist, orientation, aquery, bquery))

    # select only pairs with certain orientations - e.g. innies, outies, etc.
    if mateorientation:
        all_dist = [x for x in all_dist if x[1] == mateorientation]

    # try to infer cutoff as twice the median until convergence
    if cutoff <= 0:
        dists = np.array([x[0] for x in all_dist], dtype="int")
        p0 = analyze_dists(dists, cutoff=mpcutoff)
        cutoff = int(2 * p0)  # initial estimate
        cutoff = int(math.ceil(cutoff / bins)) * bins
        logging.debug("Insert size cutoff set to {0}, ".format(cutoff) +
            "use '--cutoff' to override")

    for dist, orientation, aquery, bquery in all_dist:
        if dist > cutoff:
            continue
        if cutoff > 2 * mpcutoff and dist < mpcutoff:
            continue

        linked_dist.append(dist)
        if pairsfile:
            print("{0}\t{1}\t{2}".format(aquery, bquery, dist), file=pairsfw)
        orientations[orientation] += 1

    print("{0} fragments, {1} pairs ({2} total)".\
                format(num_fragments, num_pairs, num_fragments + num_pairs * 2), file=sys.stderr)

    s = SummaryStats(linked_dist, dtype="int")
    num_links = s.size

    meandist, stdev = s.mean, s.sd
    p0, p1, p2 = s.median, s.p1, s.p2

    print("%d pairs (%.1f%%) are linked (cutoff=%d)" % \
            (num_links, num_links * 100. / num_pairs, cutoff), file=sys.stderr)
    print("mean distance between mates: {0} +/- {1}".\
            format(meandist, stdev), file=sys.stderr)
    print("median distance between mates: {0}".format(p0), file=sys.stderr)
    print("95% distance range: {0} - {1}".format(p1, p2), file=sys.stderr)
    print("\nOrientations:", file=sys.stderr)

    orientation_summary = []
    for orientation, count in sorted(orientations.items()):
        o = "{0}:{1}".format(orientation, \
                percentage(count, num_links, mode=1))
        orientation_summary.append(o.split()[0])
        print(o, file=sys.stderr)

    if insertsfile:
        from jcvi.graphics.histogram import histogram

        print("\n".join(str(x) for x in linked_dist), file=insertsfw)
        insertsfw.close()
        prefix = insertsfile.rsplit(".", 1)[0]
        if prefix > 10:
            prefix = prefix.split("-")[0]
        osummary = " ".join(orientation_summary)
        title="{0} ({1}; median:{2} bp)".format(prefix, osummary, p0)
        histogram(insertsfile, vmin=0, vmax=cutoff, bins=bins,
                xlabel="Insertsize", title=title, ascii=ascii)
        if op.exists(insertsfile):
            os.remove(insertsfile)

    return s
133334
def extract_pairs(fastqfile, p1fw, p2fw, fragsfw, p, suffix=False):
    
    fp = open(fastqfile)
    currentID = 0
    npairs = nfrags = 0
    for x, lib in izip(p.r1, p.libs):
        while currentID != x:
            fragsfw.writelines(islice(fp, 4))  # Exhaust the iterator
            currentID += 1
            nfrags += 1
        a = list(islice(fp, 4))
        b = list(islice(fp, 4))
        if suffix:
            name = a[0].rstrip()
            a[0] = name + "/1\n"
            b[0] = name + "/2\n"
        else:
            b[0] = a[0]  # Keep same read ID for pairs

        p1fw[lib].writelines(a)
        p2fw[lib].writelines(b)
        currentID += 2
        npairs += 2

    # Write the remaining single reads
    while True:
        contents = list(islice(fp, 4))
        if not contents:
            break
        fragsfw.writelines(contents)
        nfrags += 1

    logging.debug("A total of {0} paired reads written to `{1}`.".\
                  format(npairs, ",".join(x.name for x in p1fw + p2fw)))
    logging.debug("A total of {0} single reads written to `{1}`.".\
                  format(nfrags, fragsfw.name))

    # Validate the numbers
    expected_pairs = 2 * p.npairs
    expected_frags = p.nreads - 2 * p.npairs
    assert npairs == expected_pairs, "Expect {0} paired reads, got {1} instead".\
              format(expected_pairs, npairs)
    assert nfrags == expected_frags, "Expect {0} single reads, got {1} instead".\
              format(expected_frags, nfrags)