Subsample 1000 reads from the baited files
429608
def subsample(self):
        
        # Create the threads for the analysis
        logging.info('Subsampling FASTQ reads')
        for _ in range(self.cpus):
            threads = Thread(target=self.subsamplethreads, args=())
            threads.setDaemon(True)
            threads.start()
        with progressbar(self.runmetadata.samples) as bar:
            for sample in bar:
                if sample.general.bestassemblyfile != 'NA':
                    # Set the name of the subsampled FASTQ file
                    sample[self.analysistype].subsampledfastq = \
                        os.path.splitext(sample[self.analysistype].baitedfastq)[0] + '_subsampled.fastq'
                    # Set the system call
                    sample[self.analysistype].seqtkcall = 'reformat.sh in={} out={} samplereadstarget=1000'\
                        .format(sample[self.analysistype].baitedfastq,
                                sample[self.analysistype].subsampledfastq)
                    # Add the sample to the queue
                    self.samplequeue.put(sample)
        self.samplequeue.join()
1204834
def _fast_read(self, infile):
        """"""
        infile.seek(0)
        return(int(infile.read().decode().strip()))
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)
1025315
static long readUnary(
      final int[] compressedWords,
      final long[] ptrArr) {

    int nextWordIndex = (int) ptrArr[NEXT_WORD_IDX]; //must be int
    assert nextWordIndex == ptrArr[NEXT_WORD_IDX]; //catch truncation error
    long bitBuf = ptrArr[BIT_BUF];
    int bufBits = (int) ptrArr[BUF_BITS];

    assert compressedWords != null;
    assert nextWordIndex >= 0;
    assert bitBuf >= 0;
    assert bufBits >= 0;

    long subTotal = 0;
    int trailingZeros;

    //readUnaryLoop:
    while (true) {
      //MAYBE_FILL_BITBUF(compressedWords,nextWordIndex,8); // ensure 8 bits in bit buffer
      if (bufBits < 8) {  // Prepare for an 8-bit peek into the bitstream.
        bitBuf |= ((compressedWords[nextWordIndex++] & 0XFFFF_FFFFL) << bufBits);
        bufBits += 32;
      }

      // These 8 bits include either all or part of the Unary codeword.
      final int peek8 = (int) (bitBuf & 0XFFL);
      trailingZeros = Math.min(8, Integer.numberOfTrailingZeros(peek8));

      assert ((trailingZeros >= 0) && (trailingZeros <= 8)) : "TZ+ " + trailingZeros;

      if (trailingZeros == 8) { // The codeword was partial, so read some more.
        subTotal += 8;
        bufBits -= 8;
        bitBuf >>>= 8;
        continue;
      }
      break;
    }

    bufBits -= (1 + trailingZeros);
    bitBuf >>>= (1 + trailingZeros);

    ptrArr[NEXT_WORD_IDX] = nextWordIndex;
    assert nextWordIndex == ptrArr[NEXT_WORD_IDX]; //catch sign extension error
    ptrArr[BIT_BUF] = bitBuf;
    ptrArr[BUF_BITS] = bufBits;
    return subTotal + trailingZeros;
  }
897065
def starts_by_depth(bam_file, data, sample_size=10000000):
    
    binsize = (sample_size / 100) + 1
    seen_starts = set()
    counted = 0
    num_reads = []
    starts = []
    buffer = []
    downsampled = bam.downsample(bam_file, data, sample_size)
    with bam.open_samfile(downsampled) as samfile:
        for read in samfile:
            if read.is_unmapped:
                continue
            counted += 1
            buffer.append(str(read.tid) + ":" + str(read.pos))
            if counted % binsize == 0:
                seen_starts.update(buffer)
                buffer = []
                num_reads.append(counted)
                starts.append(len(seen_starts))
        seen_starts.update(buffer)
        num_reads.append(counted)
        starts.append(len(seen_starts))
    return pd.DataFrame({"reads": num_reads, "starts": starts})
1212024
def downsample(work_dir, sample_name, fastq_left_fpath, fastq_right_fpath, downsample_to, num_pairs=None):
    
    sample_name = sample_name or splitext(''.join(lc if lc == rc else '' for lc, rc in zip(fastq_left_fpath, fastq_right_fpath)))[0]

    l_out_fpath = make_downsampled_fpath(work_dir, fastq_left_fpath)
    r_out_fpath = make_downsampled_fpath(work_dir, fastq_right_fpath)
    if can_reuse(l_out_fpath, [fastq_left_fpath, fastq_right_fpath]):
        return l_out_fpath, r_out_fpath

    info('Processing ' + sample_name)
    if num_pairs is None:
        info(sample_name + ': counting number of reads in fastq...')
        num_pairs = _count_records_in_fastq(fastq_left_fpath)
    if num_pairs > LIMIT:
        info(sample_name + ' the number of reads is higher than ' + str(LIMIT) +
             ', sampling from only first ' + str(LIMIT))
        num_pairs = LIMIT
    info(sample_name + ': ' + str(num_pairs) + ' reads')
    num_downsample_pairs = int(downsample_to * num_pairs) if isinstance(downsample_to, float) else downsample_to
    if num_pairs <= num_downsample_pairs:
        info(sample_name + ': and it is less than ' + str(num_downsample_pairs) + ', so no downsampling.')
        return fastq_left_fpath, fastq_right_fpath
    else:
        info(sample_name + ': downsampling to ' + str(num_downsample_pairs))
        rand_records = sorted(random.sample(range(num_pairs), num_downsample_pairs))

    info('Opening ' + fastq_left_fpath)
    fh1 = open_gzipsafe(fastq_left_fpath)
    info('Opening ' + fastq_right_fpath)
    fh2 = open_gzipsafe(fastq_right_fpath) if fastq_right_fpath else None

    out_files = (l_out_fpath, r_out_fpath) if r_out_fpath else (l_out_fpath,)

    written_records = 0
    with file_transaction(work_dir, out_files) as tx_out_files:
        if isinstance(tx_out_files, six.string_types):
            tx_out_f1 = tx_out_files
        else:
            tx_out_f1, tx_out_f2 = tx_out_files
        info('Opening ' + str(tx_out_f1) + ' to write')
        sub1 = open_gzipsafe(tx_out_f1, "w")
        info('Opening ' + str(tx_out_f2) + ' to write')
        sub2 = open_gzipsafe(tx_out_f2, "w") if r_out_fpath else None
        rec_no = -1
        for rr in rand_records:
            while rec_no < rr:
                rec_no += 1
                for i in range(4): fh1.readline()
                if fh2:
                    for i in range(4): fh2.readline()
            for i in range(4):
                sub1.write(fh1.readline())
                if sub2:
                    sub2.write(fh2.readline())
            written_records += 1
            if written_records % 10000 == 0:
                info(sample_name + ': written ' + str(written_records) + ', rec_no ' + str(rec_no + 1))
            if rec_no > num_pairs:
                info(sample_name + ' reached the limit of ' + str(num_pairs), ' read lines, stopping.')
                break
        info(sample_name + ': done, written ' + str(written_records) + ', rec_no ' + str(rec_no))
        fh1.close()
        sub1.close()
        if fastq_right_fpath:
            fh2.close()
            sub2.close()

    info(sample_name + ': done downsampling, saved to ' + l_out_fpath + ' and ' + r_out_fpath + ', total ' + str(written_records) + ' paired reads written')
    return l_out_fpath, r_out_fpath
817484
private void readMbean() {
		if (clazz == null) {
			return;
		}
		try {
			if (isUnixOS) {
				long openFiles = getValue("OpenFileDescriptorCount");
				long maxOpenFiles = getValue("MaxFileDescriptorCount");

				long freePhysicalMemorySize = getValue("FreePhysicalMemorySize");
				long totalPhysicalMemorySize = getValue("TotalPhysicalMemorySize");

				long processTime = getValue("ProcessCpuTime");

				double processCPULoad = getDoubleValue("ProcessCpuLoad");
				double systemCPULoad = getDoubleValue("SystemCpuLoad");

				long processors = getValue("AvailableProcessors");

				stats.update((int) openFiles, (int) maxOpenFiles, freePhysicalMemorySize, totalPhysicalMemorySize, processTime, (int) processors, processCPULoad, systemCPULoad);

			} else if (isWindowsOS) {
				long processors = getValue("AvailableProcessors");
				double systemCPULoad = readWindowsCPULoad();

				long totalPhysicalMemorySize = readWindowsTotalMemory();
				long freePhysicalMemorySize = readWindowsFreeMemory();

				stats.update(-1, -1, freePhysicalMemorySize, totalPhysicalMemorySize, -1, (int) processors, -1, systemCPULoad);

			} else {
				long processors = getValue("AvailableProcessors");

				stats.update(-1, -1, -1, -1, -1, (int) processors, -1, -1);
			}

		} catch (Exception e) {
			log.warn("Couldn't read value due to", e);
		}
	}
605427
def corrupt_bytes(s, p=0.01, n=None):
    """"""
    s = bytes(s)
    s = array.array("B",s)
    l = len(s)
    if n is None:
        n = max(1,int(l*p))
    for i in random.sample(range(l), n):
        s[i] = (s[i]+random.randint(1,255))%256
    return s.tobytes()
92729
def run(_, data, out_dir):
    
    # logger.info("Number of aligned reads < than 0.60 in %s: %s" % (dd.get_sample_name(data), ratio))
    logger.info("Running kraken to determine contaminant: %s" % dd.get_sample_name(data))
    # ratio = bam.get_aligned_reads(bam_file, data)
    out = out_stats = None
    db = tz.get_in(["config", "algorithm", "kraken"], data)
    if db and isinstance(db, (list, tuple)):
        db = db[0]
    kraken_cmd = config_utils.get_program("kraken", data["config"])
    if db == "minikraken":
        db = os.path.join(install._get_data_dir(), "genomes", "kraken", "minikraken")

    if not os.path.exists(db):
        logger.info("kraken: no database found %s, skipping" % db)
        return {"kraken_report": "null"}

    if not os.path.exists(os.path.join(out_dir, "kraken_out")):
        work_dir = os.path.dirname(out_dir)
        utils.safe_makedir(work_dir)
        num_cores = data["config"]["algorithm"].get("num_cores", 1)
        fn_file = data["files_orig"][0] if dd.get_save_diskspace(data) else data["files"][0]
        if fn_file.endswith("bam"):
            logger.info("kraken: need fastq files as input")
            return {"kraken_report": "null"}
        with tx_tmpdir(data) as tx_tmp_dir:
            with utils.chdir(tx_tmp_dir):
                out = os.path.join(tx_tmp_dir, "kraken_out")
                out_stats = os.path.join(tx_tmp_dir, "kraken_stats")
                cat = "zcat" if fn_file.endswith(".gz") else "cat"
                cl = ("{cat} {fn_file} | {kraken_cmd} --db {db} --quick "
                      "--preload --min-hits 2 "
                      "--threads {num_cores} "
                      "--output {out} --fastq-input /dev/stdin  2> {out_stats}").format(**locals())
                do.run(cl, "kraken: %s" % dd.get_sample_name(data))
                if os.path.exists(out_dir):
                    shutil.rmtree(out_dir)
                shutil.move(tx_tmp_dir, out_dir)
    metrics = _parse_kraken_output(out_dir, db, data)
    return metrics
1223411
def _break_reads(self, contig, position, fout, min_read_length=250):
        ''''''
        sam_reader = pysam.Samfile(self.bam, "rb")
        for read in sam_reader.fetch(contig):
            seqs = []
            if read.pos < position < read.reference_end - 1:
                split_point = position - read.pos
                if split_point - 1 >= min_read_length:
                    sequence = mapping.aligned_read_to_read(read, revcomp=False, ignore_quality=not self.fastq_out).subseq(0, split_point)
                    sequence.id += '.left'
                    seqs.append(sequence)
                if read.query_length - split_point >= min_read_length:
                    sequence = mapping.aligned_read_to_read(read, revcomp=False, ignore_quality=not self.fastq_out).subseq(split_point, read.query_length)
                    sequence.id += '.right'
                    seqs.append(sequence)
            else:
                seqs.append(mapping.aligned_read_to_read(read, revcomp=False, ignore_quality=not self.fastq_out))

            for seq in seqs:
                if read.is_reverse:
                    seq.revcomp()
                print(seq, file=fout)
679806
def subsample(data, sampling):
    
    if len(data) <= sampling:
        return data
    sampling_length = int(round(len(data) / float(sampling)))
    return [mean(data[s * sampling_length:(s + 1) * sampling_length]) for s in range(0, sampling)]
925001
def get_pretrain_dataset(data, batch_size, num_ctxes, shuffle, use_avg_len,
                         num_buckets, num_parts=1, part_idx=0, prefetch=True):
    """"""
    num_files = len(glob.glob(os.path.expanduser(data)))
    logging.debug('%d files found.', num_files)
    assert num_files >= num_parts, \
        'Number of training files must be greater than the number of partitions'
    split_sampler = nlp.data.SplitSampler(num_files, num_parts=num_parts, part_index=part_idx)
    stream = nlp.data.SimpleDatasetStream(nlp.data.NumpyDataset, data, split_sampler)
    if prefetch:
        stream = nlp.data.PrefetchingStream(stream)

    def get_dataloader(dataset):
        """create data loader based on the dataset chunk"""
        lengths = dataset.get_field('valid_lengths')
        # A batch includes: input_id, masked_id, masked_position, masked_weight,
        #                   next_sentence_label, segment_id, valid_length
        batchify_fn = Tuple(Pad(), Pad(), Pad(), Pad(), Stack(), Pad(), Stack())
        if use_avg_len:
            # sharded data loader
            sampler = nlp.data.FixedBucketSampler(lengths=lengths,
                                                  # batch_size per shard
                                                  batch_size=batch_size,
                                                  num_buckets=num_buckets,
                                                  shuffle=shuffle,
                                                  use_average_length=True,
                                                  num_shards=num_ctxes)
            dataloader = nlp.data.ShardedDataLoader(dataset,
                                                    batch_sampler=sampler,
                                                    batchify_fn=batchify_fn,
                                                    num_workers=num_ctxes)
        else:
            sampler = nlp.data.FixedBucketSampler(lengths,
                                                  batch_size=batch_size * num_ctxes,
                                                  num_buckets=num_buckets,
                                                  ratio=0,
                                                  shuffle=shuffle)
            dataloader = DataLoader(dataset=dataset,
                                    batch_sampler=sampler,
                                    batchify_fn=batchify_fn,
                                    num_workers=1)
        logging.debug('Sampler created for a new dataset:\n%s', sampler.stats())
        return dataloader

    stream = stream.transform(get_dataloader)
    return stream
613879
def _proc_gnusparse_00(self, next, pax_headers, buf):
        
        offsets = []
        for match in re.finditer(br"\d+ GNU.sparse.offset=(\d+)\n", buf):
            offsets.append(int(match.group(1)))
        numbytes = []
        for match in re.finditer(br"\d+ GNU.sparse.numbytes=(\d+)\n", buf):
            numbytes.append(int(match.group(1)))
        next.sparse = list(zip(offsets, numbytes))
647928
def items(self, match="*", count=1000):
        
        for field, value in self._client.hscan_iter(
          self.key_prefix, match=match, count=count):
            yield self._decode(field), self._loads(value)
20302
def acquire_fake_data(number_of_points=1000):
    
    # time array
    t = _n.linspace(0,10,number_of_points)

    return(t, [_n.cos(t)*(1.0+0.2*_n.random.random(number_of_points)),
               _n.sin(t      +0.5*_n.random.random(number_of_points))])
582582
def read_examples(input_files, batch_size, shuffle, num_epochs=None):
  """"""
  files = []
  for e in input_files:
    for path in e.split(','):
      files.extend(file_io.get_matching_files(path))
  thread_count = multiprocessing.cpu_count()

  # The minimum number of instances in a queue from which examples are drawn
  # randomly. The larger this number, the more randomness at the expense of
  # higher memory requirements.
  min_after_dequeue = 1000

  # When batching data, the queue's capacity will be larger than the batch_size
  # by some factor. The recommended formula is (num_threads + a small safety
  # margin). For now, we use a single thread for reading, so this can be small.
  queue_size_multiplier = thread_count + 3

  # Convert num_epochs == 0 -> num_epochs is None, if necessary
  num_epochs = num_epochs or None

  # Build a queue of the filenames to be read.
  filename_queue = tf.train.string_input_producer(files, num_epochs, shuffle)

  example_id, encoded_example = tf.TextLineReader().read_up_to(
      filename_queue, batch_size)

  if shuffle:
    capacity = min_after_dequeue + queue_size_multiplier * batch_size
    return tf.train.shuffle_batch(
        [example_id, encoded_example],
        batch_size,
        capacity,
        min_after_dequeue,
        enqueue_many=True,
        num_threads=thread_count)

  else:
    capacity = queue_size_multiplier * batch_size
    return tf.train.batch(
        [example_id, encoded_example],
        batch_size,
        capacity=capacity,
        enqueue_many=True,
        num_threads=thread_count)
1338052
def bait(self, maskmiddle='f', k='19'):
        
        logging.info('Performing kmer baiting of fastq files with {at} targets'.format(at=self.analysistype))
        # There seems to be some sort of issue with java incorrectly calculating the total system memory on certain
        # computers. For now, calculate the memory, and feed it into the bbduk call
        if self.kmer_size is None:
            kmer = k
        else:
            kmer = self.kmer_size
        with progressbar(self.runmetadata) as bar:
            for sample in bar:
                if sample.general.bestassemblyfile != 'NA' and sample[self.analysistype].runanalysis:
                    # Create the folder (if necessary)
                    make_path(sample[self.analysistype].outputdir)
                    # Make the system call
                    if len(sample.general.fastqfiles) == 2:
                        # Create the command to run the baiting - paired inputs and a single, zipped output
                        sample[self.analysistype].bbdukcmd = \
                            'bbduk.sh -Xmx{mem} ref={ref} in1={in1} in2={in2} k={kmer} maskmiddle={mm} ' \
                            'threads={c} outm={om}' \
                            .format(mem=self.mem,
                                    ref=sample[self.analysistype].baitfile,
                                    in1=sample.general.trimmedcorrectedfastqfiles[0],
                                    in2=sample.general.trimmedcorrectedfastqfiles[1],
                                    kmer=kmer,
                                    mm=maskmiddle,
                                    c=str(self.cpus),
                                    om=sample[self.analysistype].baitedfastq)
                    else:
                        sample[self.analysistype].bbdukcmd = \
                            'bbduk.sh -Xmx{mem} ref={ref} in={in1} k={kmer} maskmiddle={mm} ' \
                            'threads={cpus} outm={outm}' \
                            .format(mem=self.mem,
                                    ref=sample[self.analysistype].baitfile,
                                    in1=sample.general.trimmedcorrectedfastqfiles[0],
                                    kmer=kmer,
                                    mm=maskmiddle,
                                    cpus=str(self.cpus),
                                    outm=sample[self.analysistype].baitedfastq)
                    # Run the system call (if necessary)
                    if not os.path.isfile(sample[self.analysistype].baitedfastq):
                        out, err = run_subprocess(sample[self.analysistype].bbdukcmd)
                        write_to_logfile(sample[self.analysistype].bbdukcmd,
                                         sample[self.analysistype].bbdukcmd,
                                         self.logfile, sample.general.logout, sample.general.logerr,
                                         sample[self.analysistype].logout, sample[self.analysistype].logerr)
                        write_to_logfile(out,
                                         err,
                                         self.logfile, sample.general.logout, sample.general.logerr,
                                         sample[self.analysistype].logout, sample[self.analysistype].logerr)