Filter input file, removing items with low depth of support.

    cortex_var calls are tricky to filter by depth. Count information is in
    the COV FORMAT field grouped by alleles, so we need to sum up values and
    compare.
503211
def _select_final_variants(base_vcf, out_vcf, config):
    
    min_depth = int(config["algorithm"].get("min_depth", 4))
    with file_transaction(out_vcf) as tx_out_file:
        with open(base_vcf) as in_handle:
            with open(tx_out_file, "w") as out_handle:
                for line in in_handle:
                    if line.startswith("#"):
                        passes = True
                    else:
                        passes = _passes_cortex_depth(line, min_depth)
                    if passes:
                        out_handle.write(line)
    return out_vcf
175533
def _validate_fold_has_outputs_or_count_filter(fold_scope_location, fold_has_count_filter, outputs):
    """"""
    # This function makes sure that the @fold scope has an effect.
    # Folds either output data, or filter the data enclosing the fold based on the size of the fold.
    if fold_has_count_filter:
        # This fold has a filter on the "_x_count" property, so it is legal and has an effect.
        return True

    # At least one output in the outputs list must point to the fold_scope_location,
    # or the scope corresponding to fold_scope_location had no @outputs and is illegal.
    for output in six.itervalues(outputs):
        if output['fold'] == fold_scope_location:
            return True

    raise GraphQLCompilationError(u'Found a @fold scope that has no effect on the query. '
                                  u'Each @fold scope must either perform filtering, or contain at '
                                  u'least one field marked for output. Fold location: {}'
                                  .format(fold_scope_location))
228589
def _filter_seqs(fn):
    """"""
    out_file = op.splitext(fn)[0] + "_unique.fa"
    idx = 0
    if not file_exists(out_file):
        with open(out_file, 'w') as out_handle:
            with open(fn) as in_handle:
                for line in in_handle:
                    if line.startswith("@") or line.startswith(">"):
                        fixed_name = _make_unique(line.strip(), idx)
                        seq = in_handle.next().strip()
                        counts = _get_freq(fixed_name)
                        if len(seq) < 26 and (counts > 1 or counts == 0):
                            idx += 1
                            print(fixed_name, file=out_handle, end="\n")
                            print(seq, file=out_handle, end="\n")
                        if line.startswith("@"):
                            in_handle.next()
                            in_handle.next()
    return out_file
219294
def filter_input_sequences
      # Ensure the sequences are in plain text.
      begin
        ActiveSupport::Multibyte::Unicode.u_unpack(self.sequence)
      rescue Exception => e
        logger.error e.message
        errors.add(
          :sequence,
          "Please enter your sequence(s) in Plain Text as FASTA."
        )
        self.sequence = ""
        return
      end

      hash = create_hash(self.sequence)
      tmp  = File.join(::Rails.root.to_s, 'tmp')

      # Ensure the sequences are FASTA via #write_input_sequence_to_file.
      begin
        write_input_sequence_to_file(tmp, hash, self.sequence)
      rescue Exception => e
        errors.add(:sequence, e.message)
        return
      else
        fasta = File.read(File.join(tmp, hash + ".fa"))
      ensure
        File.delete(File.join(tmp, hash + ".fa"))
        File.delete(File.join(tmp, hash + ".seq"))
      end

      self.na_sequence = ""
      self.aa_sequence = ""

      # Split the sequences on >, check the type (AA or NA) and separate.
      seqs = fasta.split('>')
      seqs.delete_if { |s| s.empty? }
      seqs.each do |s|
        type = discover_input_sequence_type(s)
        if type == "nucleic_acid"
          self.na_sequence << ">" << s
        end
        if type == "amino_acid"
          self.aa_sequence << ">" << s
        end
      end

      self.na_sequence = nil if self.na_sequence.empty?
      self.aa_sequence = nil if self.aa_sequence.empty?
    end
965861
def _cnv_prioritize(data):
    
    supported = {"cnvkit": {"inputs": ["call_file", "segmetrics"], "fn": _cnvkit_prioritize}}
    pcall = None
    priority_files = None
    for call in data.get("sv", []):
        if call["variantcaller"] in supported:
            priority_files = [call.get(x) for x in supported[call["variantcaller"]]["inputs"]]
            priority_files = [x for x in priority_files if x is not None and utils.file_exists(x)]
            if len(priority_files) == len(supported[call["variantcaller"]]["inputs"]):
                pcall = call
                break
    prioritize_by = tz.get_in(["config", "algorithm", "svprioritize"], data)
    if pcall and prioritize_by:
        out_file = "%s-prioritize.tsv" % utils.splitext_plus(priority_files[0])[0]
        gene_list = _find_gene_list_from_bed(prioritize_by, out_file, data)
        if gene_list:
            with open(gene_list) as in_handle:
                genes = [x.strip() for x in in_handle]
            args = [dd.get_sample_name(data), genes] + priority_files
            df = supported[pcall["variantcaller"]]["fn"](*args)
            with file_transaction(data, out_file) as tx_out_file:
                df.to_csv(tx_out_file, sep="\t", index=False)
            pcall["priority"] = out_file
    return data
1202506
function filterFunc(inputValue, optionValue) {
  if (inputValue && optionValue) {
    return optionValue.toLowerCase().indexOf(inputValue.toLowerCase()) > -1;
  }

  return false;
}
389100
def tagcount(sam, out, genemap, output_evidence_table, positional, minevidence,
             cb_histogram, cb_cutoff, no_scale_evidence, subsample, sparse,
             parse_tags, gene_tags):
    ''' 
    '''
    from pysam import AlignmentFile

    from io import StringIO
    import pandas as pd

    from utils import weigh_evidence

    logger.info('Reading optional files')

    gene_map = None
    if genemap:
        with open(genemap) as fh:
            try:
                gene_map = dict(p.strip().split() for p in fh)
            except ValueError:
                logger.error('Incorrectly formatted gene_map, need to be tsv.')
                sys.exit()

    if positional:
        tuple_template = '{0},{1},{2},{3}'
    else:
        tuple_template = '{0},{1},{3}'

    if not cb_cutoff:
        cb_cutoff = 0

    if cb_histogram and cb_cutoff == "auto":
        cb_cutoff = guess_depth_cutoff(cb_histogram)

    cb_cutoff = int(cb_cutoff)

    cb_hist = None
    filter_cb = False
    if cb_histogram:
        cb_hist = pd.read_csv(cb_histogram, index_col=0, header=-1, squeeze=True, sep="\t")
        total_num_cbs = cb_hist.shape[0]
        cb_hist = cb_hist[cb_hist > cb_cutoff]
        logger.info('Keeping {} out of {} cellular barcodes.'.format(cb_hist.shape[0], total_num_cbs))
        filter_cb = True

    parser_re = re.compile('.*:CELL_(?P<CB>.*):UMI_(?P<MB>.*)')

    if subsample:
        logger.info('Creating reservoir of subsampled reads ({} per cell)'.format(subsample))
        start_sampling  = time.time()

        reservoir = collections.defaultdict(list)
        cb_hist_sampled = 0 * cb_hist
        cb_obs = 0 * cb_hist

        track = stream_bamfile(sam)
        current_read = 'none_observed_yet'
        for i, aln in enumerate(track):
            if aln.qname == current_read:
                continue

            current_read = aln.qname

            if parse_tags:
                CB = aln.get_tag('CR')
            else:
                match = parser_re.match(aln.qname)
                CB = match.group('CB')

            if CB not in cb_hist.index:
                continue

            cb_obs[CB] += 1
            if len(reservoir[CB]) < subsample:
                reservoir[CB].append(i)
                cb_hist_sampled[CB] += 1
            else:
                s = pd.np.random.randint(0, cb_obs[CB])
                if s < subsample:
                    reservoir[CB][s] = i

        index_filter = set(itertools.chain.from_iterable(reservoir.values()))
        sam_file.close()
        sampling_time = time.time() - start_sampling
        logger.info('Sampling done - {:.3}s'.format(sampling_time))

    evidence = collections.defaultdict(int)

    logger.info('Tallying evidence')
    start_tally = time.time()

    sam_mode = 'r' if sam.endswith(".sam") else 'rb'
    sam_file = AlignmentFile(sam, mode=sam_mode)
    targets = [x["SN"] for x in sam_file.header["SQ"]]
    track = sam_file.fetch(until_eof=True)
    count = 0
    unmapped = 0
    kept = 0
    nomatchcb = 0
    current_read = 'none_observed_yet'
    count_this_read = True
    missing_transcripts = set()
    for i, aln in enumerate(track):
        if count and not count % 1000000:
            logger.info("Processed %d alignments, kept %d." % (count, kept))
            logger.info("%d were filtered for being unmapped." % unmapped)
            if filter_cb:
                logger.info("%d were filtered for not matching known barcodes."
                            % nomatchcb)
        count += 1

        if aln.is_unmapped:
            unmapped += 1
            continue

        if gene_tags and not aln.has_tag('GX'):
            unmapped += 1
            continue

        if aln.qname != current_read:
            current_read = aln.qname
            if subsample and i not in index_filter:
                count_this_read = False
                continue
            else:
                count_this_read = True
        else:
            if not count_this_read:
                continue

        if parse_tags:
            CB = aln.get_tag('CR')
        else:
            match = parser_re.match(aln.qname)
            CB = match.group('CB')

        if filter_cb:
            if CB not in cb_hist.index:
                nomatchcb += 1
                continue

        if parse_tags:
            MB = aln.get_tag('UM')
        else:
            MB = match.group('MB')

        if gene_tags:
            target_name = aln.get_tag('GX').split(',')[0]
        else:
            txid = sam_file.getrname(aln.reference_id)
            if gene_map:
                if txid in gene_map:
                    target_name = gene_map[txid]
                else:
                    missing_transcripts.add(txid)
                    target_name = txid
            else:
                target_name = txid

        e_tuple = tuple_template.format(CB, target_name, aln.pos, MB)

        # Scale evidence by number of hits
        if no_scale_evidence:
            evidence[e_tuple] += 1.0
        else:
            evidence[e_tuple] += weigh_evidence(aln.tags)

        kept += 1

    tally_time = time.time() - start_tally
    if missing_transcripts:
        logger.warn('The following transcripts were missing gene_ids, so we added them as the transcript ids: %s' % str(missing_transcripts))
    logger.info('Tally done - {:.3}s, {:,} alns/min'.format(tally_time, int(60. * count / tally_time)))
    logger.info('Collapsing evidence')

    logger.info('Writing evidence')
    with tempfile.NamedTemporaryFile('w+t') as out_handle:
        for key in evidence:
            line = '{},{}\n'.format(key, evidence[key])
            out_handle.write(line)

        out_handle.flush()
        out_handle.seek(0)

        evidence_table = pd.read_csv(out_handle, header=None)

    del evidence

    evidence_query = 'evidence >= %f' % minevidence
    if positional:
        evidence_table.columns=['cell', 'gene', 'umi', 'pos', 'evidence']
        collapsed = evidence_table.query(evidence_query).groupby(['cell', 'gene'])['umi', 'pos'].size()

    else:
        evidence_table.columns=['cell', 'gene', 'umi', 'evidence']
        collapsed = evidence_table.query(evidence_query).groupby(['cell', 'gene'])['umi'].size()

    expanded = collapsed.unstack().T

    if gene_map:
        # This Series is just for sorting the index
        genes = pd.Series(index=set(gene_map.values()))
        genes = genes.sort_index()
        # Now genes is assigned to a DataFrame
        genes = expanded.ix[genes.index]

    elif gene_tags:
        expanded.sort_index()
        genes = expanded

    else:
        # make data frame have a complete accounting of transcripts
        targets = pd.Series(index=set(targets))
        targets = targets.sort_index()
        expanded = expanded.reindex(targets.index.values, fill_value=0)
        genes = expanded

    genes.fillna(0, inplace=True)
    genes = genes.astype(int)
    genes.index.name = "gene"

    logger.info('Output results')

    if subsample:
        cb_hist_sampled.to_csv('ss_{}_'.format(subsample) + os.path.basename(cb_histogram), sep='\t')

    if output_evidence_table:
        import shutil
        buf.seek(0)
        with open(output_evidence_table, 'w') as etab_fh:
            shutil.copyfileobj(buf, etab_fh)

    if sparse:
        pd.Series(genes.index).to_csv(out + ".rownames", index=False, header=False)
        pd.Series(genes.columns.values).to_csv(out + ".colnames", index=False, header=False)
        with open(out, "w+b") as out_handle:
            scipy.io.mmwrite(out_handle, scipy.sparse.csr_matrix(genes))
    else:
        genes.to_csv(out)
789033
def count(self, cls, filters, quick=True, sort_by=None, select=None):
        
        query = "select count(*) from `%s` %s" % (self.domain.name, self._build_filter_part(cls, filters, sort_by, select))
        count = 0
        for row in self.domain.select(query):
            count += int(row['Count'])
            if quick:
                return count
        return count
458921
def count(self, val=True):
        
        return sum((elem.count(val) for elem in self._iter_components()))
966188
def forall(self, vars_list: List[str]) -> 'TensorFluent':
        '''
        '''
        return self._aggregation_op(tf.reduce_all, self, vars_list)
833882
def filter_single_reads_by_length(in_file, quality_format, min_length=20,
                                  out_file=None):
    
    logger.info("Removing reads in %s thare are less than %d bases."
                % (in_file, min_length))
    in_iterator = SeqIO.parse(in_file, quality_format)
    out_iterator = (record for record in in_iterator if
                    len(record.seq) > min_length)
    with file_transaction(out_file) as tmp_out_file:
        with open(tmp_out_file, "w") as out_handle:
            SeqIO.write(out_iterator, out_handle, quality_format)
    return out_file
23619
function map(cov) {
  var ret = {
      instrumentation: 'node-jscoverage'
    , sloc: 0
    , hits: 0
    , misses: 0
    , coverage: 0
    , files: []
  };

  for (var filename in cov) {
    var data = coverage(filename, cov[filename]);
    ret.files.push(data);
    ret.hits += data.hits;
    ret.misses += data.misses;
    ret.sloc += data.sloc;
  }

  ret.files.sort(function(a, b) {
    return a.filename.localeCompare(b.filename);
  });

  if (ret.sloc > 0) {
    ret.coverage = (ret.hits / ret.sloc) * 100;
  }

  return ret;
}
629432
def _delly_count_evidence_filter(in_file, data):
    
    filtname = "DVSupport"
    filtdoc = "FMT/DV < 4 || (SR < 1 && PE < 5) || (SR < 5 && PE < 1)"
    out_file = "%s-filter%s" % utils.splitext_plus(in_file)
    cur_out_file = out_file.replace(".vcf.gz", ".vcf")
    if not utils.file_exists(out_file):
        with file_transaction(data, cur_out_file) as tx_out_file:
            with utils.open_gzipsafe(in_file) as in_handle:
                with open(tx_out_file, "w") as out_handle:
                    inp = vcf.Reader(in_handle, in_file)
                    inp.filters["DVSupport"] = vcf.parser._Filter(filtname, filtdoc)
                    outp = vcf.Writer(out_handle, inp)
                    for rec in inp:
                        sr = rec.INFO.get("SR", 0)
                        pe = rec.INFO.get("PE", 0)
                        call = rec.samples[0].data
                        dv = call.DV if hasattr(call, "DV") else 0
                        if dv < 4 or (sr < 1 and pe < 5) or (sr < 5 and pe < 1):
                            rec.add_filter(filtname)
                        outp.write_record(rec)
    if out_file.endswith(".vcf.gz"):
        out_file = vcfutils.bgzip_and_index(cur_out_file, data["config"])
    return out_file
234561
def _filter_by_background(base_name, back_samples, gt_vcfs, data):
    
    filtname = "InBackground"
    filtdoc = "Variant also present in background samples with same genotype"
    orig_vcf = gt_vcfs[base_name]
    out_file = "%s-backfilter.vcf" % (utils.splitext_plus(orig_vcf)[0])
    if not utils.file_exists(out_file) and not utils.file_exists(out_file + ".gz"):
        with file_transaction(data, out_file) as tx_out_file:
            with utils.open_gzipsafe(orig_vcf) as in_handle:
                inp = vcf.Reader(in_handle, orig_vcf)
                inp.filters[filtname] = vcf.parser._Filter(filtname, filtdoc)
                with open(tx_out_file, "w") as out_handle:
                    outp = vcf.Writer(out_handle, inp)
                    for rec in inp:
                        if _genotype_in_background(rec, base_name, back_samples):
                            rec.add_filter(filtname)
                        outp.write_record(rec)
    if utils.file_exists(out_file + ".gz"):
        out_file = out_file + ".gz"
    gt_vcfs[base_name] = vcfutils.bgzip_and_index(out_file, data["config"])
    return gt_vcfs
956727
function filterByThreshold(nodeModel, totalArea, sum, orderBy, orderedChildren) {

    // visibleMin is not supported yet when no option.sort.
    if (!orderBy) {
        return sum;
    }

    var visibleMin = nodeModel.get('visibleMin');
    var len = orderedChildren.length;
    var deletePoint = len;

    // Always travel from little value to big value.
    for (var i = len - 1; i >= 0; i--) {
        var value = orderedChildren[
            orderBy === 'asc' ? len - i - 1 : i
        ].getValue();

        if (value / sum * totalArea < visibleMin) {
            deletePoint = i;
            sum -= value;
        }
    }

    orderBy === 'asc'
        ? orderedChildren.splice(0, len - deletePoint)
        : orderedChildren.splice(deletePoint, len - deletePoint);

    return sum;
}
1271625
def _remove_prioritization(in_file, data, out_dir=None):
    
    out_file = "%s-germline.vcf" % utils.splitext_plus(in_file)[0]
    if out_dir:
        out_file = os.path.join(out_dir, os.path.basename(out_file))
    if not utils.file_uptodate(out_file, in_file) and not utils.file_uptodate(out_file + ".gz", in_file):
        with file_transaction(data, out_file) as tx_out_file:
            reader = cyvcf2.VCF(str(in_file))
            reader.add_filter_to_header({'ID': 'Somatic', 'Description': 'Variant called as Somatic'})
            # with open(tx_out_file, "w") as out_handle:
            #     out_handle.write(reader.raw_header)
            with contextlib.closing(cyvcf2.Writer(tx_out_file, reader)) as writer:
                for rec in reader:
                    rec = _update_prioritization_filters(rec)
                    # out_handle.write(str(rec))
                    writer.write_record(rec)
    return out_file
828640
function filterFiles(files, recursive) {
    var fileOrDir = function fileOrDir(lstat) {
        return lstat.isFile() || lstat.isDirectory();
    };
    return files.filter(function (_ref) {
        var lstat = _ref.lstat;

        return recursive ? fileOrDir(lstat) : lstat.isFile();
    });
}