%prog filter *.consensus.fasta

    Filter consensus sequence with min cluster size.
4519
def filter(args):
    
    from jcvi.formats.fasta import Fasta, SeqIO

    p = OptionParser(filter.__doc__)
    p.add_option("--minsize", default=2, type="int",
                 help="Minimum cluster size")
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    fastafiles = args
    minsize = opts.minsize
    totalreads = totalassembled = 0
    fw = must_open(opts.outfile, "w")
    for i, fastafile in enumerate(fastafiles):
        f = Fasta(fastafile, lazy=True)
        pf = "s{0:03d}".format(i)
        nreads = nsingletons = nclusters = 0
        for desc, rec in f.iterdescriptions_ordered():
            nclusters += 1
            if desc.startswith("singleton"):
                nsingletons += 1
                nreads += 1
                continue
            # consensus_for_cluster_0 with 63 sequences
            name, w, size, seqs = desc.split()
            assert w == "with"
            size = int(size)
            nreads += size
            if size < minsize:
                continue
            rec.description = rec.description.split(None, 1)[-1]
            rec.id = pf + "_" + rec.id
            SeqIO.write(rec, fw, "fasta")
        logging.debug("Scanned {0} clusters with {1} reads ..".\
                       format(nclusters, nreads))
        cclusters, creads = nclusters - nsingletons, nreads - nsingletons
        logging.debug("Saved {0} clusters (min={1}) with {2} reads (avg:{3}) [{4}]".\
                       format(cclusters, minsize, creads, creads / cclusters, pf))
        totalreads += nreads
        totalassembled += nreads - nsingletons
    logging.debug("Total assembled: {0}".\
                  format(percentage(totalassembled, totalreads)))
262000
def filtervcf(args):
    
    p = OptionParser(filtervcf.__doc__)
    p.set_home("lobstr", default="/mnt/software/lobSTR")
    p.set_aws_opts(store="hli-mv-data-science/htang/str")
    p.set_cpus()
    opts, args = p.parse_args(args)

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

    samples, = args
    lhome = opts.lobstr_home
    store = opts.output_path

    if samples.endswith((".vcf", ".vcf.gz")):
        vcffiles = [samples]
    else:
        vcffiles = [x.strip() for x in must_open(samples)]

    vcffiles = [x for x in vcffiles if ".filtered." not in x]

    run_args = [(x, lhome, x.startswith("s3://") and store) for x in vcffiles]
    cpus = min(opts.cpus, len(run_args))
    p = Pool(processes=cpus)
    for res in p.map_async(run_filter, run_args).get():
        continue
1027106
def _filter_by_moys_slow(self, moys):
        """"""
        _filt_values = []
        _filt_datetimes = []
        for i, d in enumerate(self.datetimes):
            if d.moy in moys:
                _filt_datetimes.append(d)
                _filt_values.append(self._values[i])
        return _filt_values, _filt_datetimes
284926
def mb_filter(fastq, cores):
    ''' 
    '''
    filter_mb = partial(umi_filter)
    p = multiprocessing.Pool(cores)

    chunks = tz.partition_all(10000, read_fastq(fastq))
    bigchunks = tz.partition_all(cores, chunks)
    for bigchunk in bigchunks:
        for chunk in p.map(filter_mb, list(bigchunk)):
            for read in chunk:
                sys.stdout.write(read)
513323
def cluster(args):
    

    args = _check_args(args)
    read_stats_file = op.join(args.dir_out, "read_stats.tsv")
    if file_exists(read_stats_file):
        os.remove(read_stats_file)

    bam_file, seq_obj = _clean_alignment(args)

    logger.info("Parsing matrix file")
    seqL, y, l = parse_ma_file(seq_obj, args.ffile)
    # y, l = _total_counts(seqL.keys(), seqL)
    logger.info("counts after: %s" % sum(y.values()))
    logger.info("# sequences after: %s" % l)
    dt = pd.DataFrame({'sample': y.keys(), 'counts': y.values()})
    dt['step'] = 'aligned'
    dt.to_csv(read_stats_file, sep="\t", index=False, header=False, mode='a')

    if len(seqL.keys()) < 10:
        logger.error("It seems you have low coverage. Please check your fastq files have enough sequences.")
        raise ValueError("So few sequences.")

    logger.info("Cleaning bam file")
    y, l = _total_counts(seqL.keys(), seqL)
    logger.info("counts after: %s" % sum(y.values()))
    logger.info("# sequences after: %s" % l)
    dt = pd.DataFrame({'sample': y.keys(), 'counts': y.values()})
    dt['step'] = 'cleaned'
    dt.to_csv(read_stats_file, sep="\t", index=False, header=False, mode='a')

    clusL = _create_clusters(seqL, bam_file, args)
    y, l = _total_counts(clusL.seq.keys(), clusL.seq, aligned=True)
    logger.info("counts after: %s" % sum(y.values()))
    logger.info("# sequences after: %s" % l)
    dt = pd.DataFrame({'sample': y.keys(), 'counts': y.values()})
    dt['step'] = 'clusters'
    dt.to_csv(read_stats_file, sep="\t", index=False, header=False, mode='a')

    logger.info("Solving multi-mapping events in the network of clusters")
    clusLred = _cleaning(clusL, args.dir_out)
    y, l = _total_counts(clusLred.clus, seqL)
    logger.info("counts after: %s" % sum(y.values()))
    logger.info("# sequences after: %s" % l)
    dt = pd.DataFrame({'sample': y.keys(), 'counts': y.values()})
    dt['step'] = 'meta-cluster'
    dt.to_csv(read_stats_file, sep="\t", index=False, header=False, mode='a')
    logger.info("Clusters up to %s" % (len(clusLred.clus.keys())))

    if args.show:
        logger.info("Creating sequences alignment to precursor")
        clusLred = show_seq(clusLred, args.index)

    clusLred = peak_calling(clusLred)

    clusLred = _annotate(args, clusLred)

    logger.info("Creating json and count matrix")
    json_file = _create_json(clusLred, args)

    logger.info("Output file in: %s" % args.dir_out)

    if args.db:
        name = args.db + ".db"
        logger.info("Create database: database/" + name)
        data = load_data(json_file)
        out_dir = op.join(args.dir_out, "database")
        make_database(data, name, out_dir)
    logger.info("Finished")
497162
def blat(args):
    
    p = OptionParser(blat.__doc__)
    p.add_option("--minscore", default=100, type="int",
                 help="Matches minus mismatches gap penalty [default: %default]")
    p.add_option("--minid", default=98, type="int",
                 help="Minimum sequence identity [default: %default]")
    p.set_cpus()
    opts, args = p.parse_args(args)

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

    oldfasta, newfasta = args
    twobitfiles = []
    for fastafile in args:
        tbfile = faToTwoBit(fastafile)
        twobitfiles.append(tbfile)

    oldtwobit, newtwobit = twobitfiles
    cmd = "pblat -threads={0}".format(opts.cpus) if which("pblat") else "blat"
    cmd += " {0} {1}".format(oldtwobit, newfasta)
    cmd += " -tileSize=12 -minScore={0} -minIdentity={1} ".\
                format(opts.minscore, opts.minid)
    pslfile = "{0}.{1}.psl".format(*(op.basename(x).split('.')[0] \
                for x in (newfasta, oldfasta)))
    cmd += pslfile
    sh(cmd)
769588
def cb_filter(fastq, bc1, bc2, bc3, cores, nedit):
    ''' 
    '''

    with open_gzipsafe(bc1) as bc1_fh:
        bc1 = set(cb.strip() for cb in bc1_fh)
    if bc2:
        with open_gzipsafe(bc2) as bc2_fh:
            bc2 = set(cb.strip() for cb in bc2_fh)
    if bc3:
        with open_gzipsafe(bc3) as bc3_fh:
            bc3 = set(cb.strip() for cb in bc3_fh)

    annotations = detect_fastq_annotations(fastq)
    re_string = construct_transformed_regex(annotations)

    if nedit == 0:
        filter_cb = partial(exact_barcode_filter, bc1=bc1, bc2=bc2, bc3=bc3,
                            re_string=re_string)
    else:
        bc1hash = MutationHash(bc1, nedit)
        bc2hash = None
        bc3hash = None
        if bc2:
            bc2hash = MutationHash(bc2, nedit)
        if bc3:
            bc3hash = MutationHash(bc3, nedit)
        filter_cb = partial(correcting_barcode_filter, bc1hash=bc1hash,
                            bc2hash=bc2hash, bc3hash=bc3hash, re_string=re_string)
    p = multiprocessing.Pool(cores)

    chunks = tz.partition_all(10000, read_fastq(fastq))
    bigchunks = tz.partition_all(cores, chunks)
    for bigchunk in bigchunks:
        for chunk in p.map(filter_cb, list(bigchunk)):
            for read in chunk:
                sys.stdout.write(read)
579905
def filter(args):
    
    p = OptionParser(filter.__doc__)
    p.add_option("--less", default=False, action="store_true",
                 help="filter the sizes < certain cutoff [default: >=]")
    p.set_outfile()

    opts, args = p.parse_args(args)

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

    fastafile, cutoff = args
    try:
        cutoff = int(cutoff)
    except ValueError:
        sys.exit(not p.print_help())

    f = Fasta(fastafile, lazy=True)

    fw = must_open(opts.outfile, "w")
    for name, rec in f.iteritems_ordered():

        if opts.less and len(rec) >= cutoff:
            continue

        if (not opts.less) and len(rec) < cutoff:
            continue

        SeqIO.write([rec], fw, "fasta")
        fw.flush()

    return fw.name
345984
def filter_genes_fano_deprecated(X, Ecutoff, Vcutoff):
    
    if issparse(X):
        raise ValueError('Not defined for sparse input. See `filter_genes_dispersion`.')
    mean_filter = np.mean(X, axis=0) > Ecutoff
    var_filter = np.var(X, axis=0) / (np.mean(X, axis=0) + .0001) > Vcutoff
    gene_subset = np.nonzero(np.all([mean_filter, var_filter], axis=0))[0]
    return gene_subset
512101
def partition(predicate, iterable):
    

    t1, t2 = tee(iterable)

    return filterfalse(predicate, t1), filter(predicate, t2)
278027
def get_consensus_hashes(cls, impl, working_dir, start_block_height=None, end_block_height=None, db_con=None, completeness_check=True):
        
        if (start_block_height is None and end_block_height is not None) or (start_block_height is not None and end_block_height is None):
            raise ValueError("Need either both or neither start/end block height")

        we_opened = False
        if db_con is not None:
            con = db_con

        else:
            assert impl and working_dir, 'Need impl and working_dir if db_con is not given'
            con = cls.db_open(impl, working_dir)
            we_opened = True

        range_query = ''
        range_args = ()

        if start_block_height and end_block_height:
            range_query += ' WHERE block_id >= ? AND block_id < ?'
            range_args += (start_block_height,end_block_height)

        query = 'SELECT block_id,consensus_hash FROM snapshots' + range_query + ';'
        args = range_args

        rows = cls.db_query_execute(con, query, range_args, verbose=False)
        ret = {}
        block_min = None
        block_max = None

        for r in rows:
            block_id = int(r['block_id'])
            ret[block_id] = r['consensus_hash']

            if block_min is None or block_min > block_id:
                block_min = block_id

            if block_max is None or block_max < block_id:
                block_max = block_id

        if we_opened:
            con.close()
       
        if completeness_check:
            # sanity check
            for i in range(block_min,block_max+1):
                ch = ret.get(i, None)

                assert ch is not None, 'Missing consensus hash for {}'.format(i)
                assert isinstance(ret[i], (str,unicode)), 'consensus hash for {} is type {}'.format(i, type(ret[i]))

        return ret
511061
def fasta(args):
    
    from jcvi.formats.sizes import Sizes
    from jcvi.formats.agp import OO, build

    p = OptionParser(fasta.__doc__)
    opts, args = p.parse_args(args)

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

    bedfile, scffasta, pmolfasta = args
    pf = bedfile.rsplit(".", 1)[0]
    bed = Bed(bedfile)
    selected = select_bed(bed)
    oo = OO()
    seen = set()
    sizes = Sizes(scffasta).mapping
    agpfile = pf + ".agp"
    agp = open(agpfile, "w")
    for b in selected:
        scf = range_parse(b.accn).seqid
        chr = b.seqid
        cs = (chr, scf)
        if cs not in seen:
            oo.add(chr, scf, sizes[scf], b.strand)
            seen.add(cs)
        else:
            logging.debug("Seen {0}, ignored.".format(cs))

    oo.write_AGP(agp, gaptype="contig")
    agp.close()
    build([agpfile, scffasta, pmolfasta])
444353
def filter(polylines, min_len=20):
    
    filtered = []
    for n in range(len(polylines) - 1, -1, -1):
        if lengths(polylines[n]).sum() < min_len:
            filtered.append(polylines.pop(n))
    return filtered
932696
def fastFilter(arr, ksize=30, every=None, resize=True, fn='median',
                  interpolation=cv2.INTER_LANCZOS4,
                  smoothksize=0,
                  borderMode=cv2.BORDER_REFLECT):
    '''
    
    '''
    if every is None:
        every = max(ksize//3, 1)
    else:
        assert ksize >= 3*every
    s0,s1 = arr.shape[:2]
    
    ss0 = s0//every
    every = s0//ss0
    ss1 = s1//every
    
    out = np.full((ss0+1,ss1+1), np.nan)
    
    c = {'median':_calcMedian,
         'nanmedian':_calcNanMedian,
         'nanmean':_calcNanMean,
         'mean':_calcMean,
         }[fn]
    ss0,ss1 = c(arr, out, ksize, every)
    out = out[:ss0,:ss1]
    
    if smoothksize:
        out = gaussian_filter(out, smoothksize)
        
    
    if not resize:
        return out
    return cv2.resize(out, arr.shape[:2][::-1],
               interpolation=interpolation)
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())
92186
def parallel_fn(f):
    

    def simple_parallel(func, sequence, **args):
        """ f takes an element of sequence as input and the keyword args in **args"""
        if 'processes' in args:
            processes = args.get('processes')
            del args['processes']
        else:
            processes = multiprocessing.cpu_count()

        pool = multiprocessing.Pool(processes)  # depends on available cores

        result = pool.map_async(universal_worker, pool_args(func, sequence, args))
        pool.close()
        pool.join()
        cleaned = [x for x in result.get() if x is not None]  # getting results
        cleaned = asarray(cleaned)
        return cleaned
    return partial(simple_parallel, f)
229585
def filter_seq(seq):
     ''''''
     if seq.res:
          return None
     n = nt.Factors(seq.factors)
     guide, s, t = aq.canonical_form(n)
     seq.guide = guide
     # The target_tau for the composite is at most the class minus extant prime factor count
     cls = aq.get_class(guide=guide)
     num_larges = seq.factors.count('P')
     upper_bound_tau = cls - num_larges - len(t)

     if cls < 2 or upper_bound_tau < 2: # Cheap tests to eliminate almost all sequences
          return None

     # Next we ignore sequences whose guide is solely a power of 2 greater than 3
     v = nt.Factors({p: a for p, a in guide.items() if p != 2 and a > 0})
     if int(v) == 1 and cls > 3:
          return None
     # This condition greatly reduces fdb load, but excludes a lot of sequences
     if not aq.is_driver(guide=guide):
          return None

     return n, guide