Use bbnorm from the bbmap suite of tools to perform read normalisation
186677
def normalise_reads(self):
        
        logging.info('Normalising reads to a kmer depth of 100')
        for sample in self.metadata:
            # Set the name of the normalised read files
            sample.general.normalisedreads = [fastq.split('.fastq.gz')[0] + '_normalised.fastq.gz'
                                              for fastq in sorted(sample.general.fastqfiles)]
            try:
                # Run the normalisation command
                out, err, cmd = bbtools.bbnorm(forward_in=sorted(sample.general.trimmedcorrectedfastqfiles)[0],
                                               forward_out=sample.general.normalisedreads[0],
                                               returncmd=True,
                                               threads=self.cpus)
                sample[self.analysistype].normalisecmd = cmd
                write_to_logfile(out, err, self.logfile, sample.general.logout, sample.general.logerr, None, None)
            except CalledProcessError:
                sample.general.normalisedreads = sample.general.trimmedfastqfiles
            except IndexError:
                sample.general.normalisedreads = list()
1162886
def select_unaligned_read_pairs(in_bam, extra, out_dir, config):
    
    runner = broad.runner_from_config(config)
    base, ext = os.path.splitext(os.path.basename(in_bam))
    nomap_bam = os.path.join(out_dir, "{}-{}{}".format(base, extra, ext))
    if not utils.file_exists(nomap_bam):
        with file_transaction(nomap_bam) as tx_out:
            runner.run("FilterSamReads", [("INPUT", in_bam),
                                          ("OUTPUT", tx_out),
                                          ("EXCLUDE_ALIGNED", "true"),
                                          ("WRITE_READS_FILES", "false"),
                                          ("SORT_ORDER", "queryname")])
    has_reads = False
    with pysam.Samfile(nomap_bam, "rb") as in_pysam:
        for read in in_pysam:
            if read.is_paired:
                has_reads = True
                break
    if has_reads:
        out_fq1, out_fq2 = ["{}-{}.fq".format(os.path.splitext(nomap_bam)[0], i) for i in [1, 2]]
        runner.run_fn("picard_bam_to_fastq", nomap_bam, out_fq1, out_fq2)
        return out_fq1, out_fq2
    else:
        return None, None
641892
def _get_region(self, contig, start, end, fout, min_length=250):
        ''''''
        sam_reader = pysam.Samfile(self.bam, "rb")
        trimming_end = (start == 0)
        for read in sam_reader.fetch(contig, start, end):
            read_interval = pyfastaq.intervals.Interval(read.pos, read.reference_end - 1)
            seq = mapping.aligned_read_to_read(read, ignore_quality=not self.fastq_out, revcomp=False)

            if trimming_end:
                bases_off_start = 0
                bases_off_end = max(0, read.reference_end - 1 - end)
                #seq.seq = seq.seq[:read.query_alignment_end - bases_off_end]
                seq = seq.subseq(0, read.query_alignment_end - bases_off_end)
            else:
                bases_off_start = max(0, start - read.pos + 1)
                #seq.seq = seq.seq[bases_off_start  + read.query_alignment_start:]
                seq = seq.subseq(bases_off_start  + read.query_alignment_start, len(seq))

            if read.is_reverse:
                seq.revcomp()

            if len(seq) >= min_length:
                print(seq, file=fout)
40600
def normalize_bbox(coords, ymax, scaler=2):
    
    return [
        coords[0] * scaler,
        ymax - (coords[3] * scaler),
        coords[2] * scaler,
        ymax - (coords[1] * scaler),
    ]
875956
def read_quality_trim(self):
        
        logging.info('Quality trim')
        for sample in self.metadata:
            sample.sampled_reads = GenObject()
            sample.sampled_reads.outputdir = os.path.join(sample.outputdir, 'sampled')
            sample.sampled_reads.trimmed_dir = os.path.join(sample.sampled_reads.outputdir, 'qualitytrimmed_reads')
            make_path(sample.sampled_reads.trimmed_dir)
            for depth in self.read_depths:
                # Create the depth GenObject
                setattr(sample.sampled_reads, depth, GenObject())
                # Set the depth and output directory attributes for the depth GenObject
                sample.sampled_reads[depth].depth = depth
                sample.sampled_reads[depth].depth_dir = os.path.join(sample.sampled_reads.outputdir, depth)
                # Create the output directory
                make_path(sample.sampled_reads[depth].depth_dir)
                for read_pair in self.read_lengths:
                    # Create the read_pair GenObject within the depth GenObject
                    setattr(sample.sampled_reads[depth], read_pair, GenObject())
                    # Set and create the output directory
                    sample.sampled_reads[depth][read_pair].outputdir = \
                        os.path.join(sample.sampled_reads[depth].depth_dir, read_pair)
                    make_path(sample.sampled_reads[depth][read_pair].outputdir)
                    # Create both forward_reads and reverse_reads sub-GenObjects
                    sample.sampled_reads[depth][read_pair].forward_reads = GenObject()
                    sample.sampled_reads[depth][read_pair].reverse_reads = GenObject()
                    sample.sampled_reads[depth][read_pair].trimmed_dir = \
                        os.path.join(sample.sampled_reads.trimmed_dir,
                                     read_pair)
                    make_path(sample.sampled_reads[depth][read_pair].trimmed_dir)
                    # Extract the forward and reverse reads lengths from the read_pair variable
                    sample.sampled_reads[depth][read_pair].forward_reads.length, \
                        sample.sampled_reads[depth][read_pair].reverse_reads.length = read_pair.split('_')
                    logging.info('Performing quality trimming on reads from sample {name} at depth {depth} '
                                 'for minimum read length {forward}'
                                 .format(name=sample.name,
                                         depth=depth,
                                         forward=sample.sampled_reads[depth][read_pair].forward_reads.length))
                    # Set the attributes for the trimmed forward and reverse reads to use for subsampling
                    sample.sampled_reads[depth][read_pair].trimmed_forwardfastq = \
                        os.path.join(sample.sampled_reads[depth][read_pair].trimmed_dir,
                                     '{name}_{length}_R1.fastq.gz'
                                     .format(name=sample.name,
                                             length=sample.sampled_reads[depth][read_pair].forward_reads.length))
                    sample.sampled_reads[depth][read_pair].trimmed_reversefastq = \
                        os.path.join(sample.sampled_reads[depth][read_pair].trimmed_dir,
                                     '{name}_{length}_R2.fastq.gz'
                                     .format(name=sample.name,
                                             length=sample.sampled_reads[depth][read_pair].forward_reads.length))
                    # Create the trimmed output directory attribute
                    sample.sampled_reads[depth][read_pair].sampled_trimmed_outputdir \
                        = os.path.join(sample.sampled_reads[depth][read_pair].outputdir,
                                       'sampled_trimmed')
                    # Set the name of the forward trimmed reads - include the depth and read length information
                    # This is set now, as the untrimmed files will be removed, and a check is necessary
                    sample.sampled_reads[depth][read_pair].forward_reads.trimmed_sampled_fastq = \
                        os.path.join(sample.sampled_reads[depth][read_pair].sampled_trimmed_outputdir,
                                     '{name}_sampled_{depth}_{read_pair}_R1.fastq.gz'
                                     .format(name=sample.name,
                                             depth=depth,
                                             read_pair=read_pair))
                    # Reverse reads
                    sample.sampled_reads[depth][read_pair].reverse_reads.trimmed_sampled_fastq = \
                        os.path.join(sample.sampled_reads[depth][read_pair].sampled_trimmed_outputdir,
                                     '{name}_sampled_{depth}_{read_pair}_R2.fastq.gz'
                                     .format(name=sample.name,
                                             depth=depth,
                                             read_pair=read_pair))
                    # Sample if the forward output file does not already exist
                    if not os.path.isfile(sample.sampled_reads[depth][read_pair].trimmed_forwardfastq) and \
                            not os.path.isfile(
                                sample.sampled_reads[depth][read_pair].forward_reads.trimmed_sampled_fastq):
                        out, \
                            err, \
                            sample.sampled_reads[depth][read_pair].sample_cmd = \
                            bbtools.bbduk_trim(forward_in=sample.forward_fastq,
                                               forward_out=sample.sampled_reads[depth][read_pair]
                                               .trimmed_forwardfastq,
                                               reverse_in=sample.reverse_fastq,
                                               reverse_out=sample.sampled_reads[depth][read_pair]
                                               .trimmed_reversefastq,
                                               minlength=sample.sampled_reads[depth][read_pair]
                                               .forward_reads.length,
                                               forcetrimleft=0,
                                               returncmd=True,
                                               **{'ziplevel': '9',
                                                  'Xmx': self.mem})
            # Update the JSON file
            self.write_json(sample)
1071523
def _get_all_unmapped_reads(self, fout):
        ''''''
        sam_reader = pysam.Samfile(self.bam, "rb")
        for read in sam_reader.fetch(until_eof=True):
            if read.is_unmapped:
                print(mapping.aligned_read_to_read(read, ignore_quality=not self.fastq_out), file=fout)
486473
def tiered_alignment(in_bam, tier_num, multi_mappers, extra_args,
                     genome_build, pair_stats,
                     work_dir, dirs, config):
    
    nomap_fq1, nomap_fq2 = select_unaligned_read_pairs(in_bam, "tier{}".format(tier_num),
                                                       work_dir, config)
    if nomap_fq1 is not None:
        base_name = "{}-tier{}out".format(os.path.splitext(os.path.basename(in_bam))[0],
                                          tier_num)
        config = copy.deepcopy(config)
        dirs = copy.deepcopy(dirs)
        config["algorithm"]["bam_sort"] = "queryname"
        config["algorithm"]["multiple_mappers"] = multi_mappers
        config["algorithm"]["extra_align_args"] = ["-i", int(pair_stats["mean"]),
                                               int(pair_stats["std"])] + extra_args
        out_bam, ref_file = align_to_sort_bam(nomap_fq1, nomap_fq2,
                                              lane.rg_names(base_name, base_name, config),
                                              genome_build, "novoalign",
                                              dirs, config,
                                              dir_ext=os.path.join("hydra", os.path.split(nomap_fq1)[0]))
        return out_bam
    else:
        return None
348522
def postprocess_alignment(data):
    
    data = cwlutils.normalize_missing(utils.to_single_data(data))
    data = cwlutils.unpack_tarballs(data, data)
    bam_file = data.get("align_bam") or data.get("work_bam")
    ref_file = dd.get_ref_file(data)
    if vmulti.bam_needs_processing(data) and bam_file and bam_file.endswith(".bam"):
        out_dir = utils.safe_makedir(os.path.join(dd.get_work_dir(data), "align",
                                                  dd.get_sample_name(data)))
        bam_file_ready = os.path.join(out_dir, os.path.basename(bam_file))
        if not utils.file_exists(bam_file_ready):
            utils.symlink_plus(bam_file, bam_file_ready)
        bam.index(bam_file_ready, data["config"])
        covinfo = callable.sample_callable_bed(bam_file_ready, ref_file, data)
        callable_region_bed, nblock_bed = \
            callable.block_regions(covinfo.raw_callable, bam_file_ready, ref_file, data)
        data["regions"] = {"nblock": nblock_bed,
                           "callable": covinfo.raw_callable,
                           "sample_callable": covinfo.callable,
                           "mapped_stats": readstats.get_cache_file(data)}
        data["depth"] = covinfo.depth_files
        data = coverage.assign_interval(data)
        data = samtools.run_and_save(data)
        data = recalibrate.prep_recal(data)
        data = recalibrate.apply_recal(data)
    elif dd.get_variant_regions(data):
        callable_region_bed, nblock_bed = \
            callable.block_regions(dd.get_variant_regions(data), bam_file, ref_file, data)
        data["regions"] = {"nblock": nblock_bed, "callable": dd.get_variant_regions(data),
                           "sample_callable": dd.get_variant_regions(data)}
    return [[data]]
966079
def _replace_bm(self):
        """"""
        self._block_matcher = cv2.StereoBM(preset=self._bm_preset,
                                          ndisparities=self._search_range,
                                          SADWindowSize=self._window_size)
665016
def normalize(s, replace_spaces=True):
    
    whitelist = (' -' + string.ascii_letters + string.digits)

    if type(s) == six.binary_type:
        s = six.text_type(s, 'utf-8', 'ignore')

    table = {}
    for ch in [ch for ch in s if ch not in whitelist]:
        if ch not in table:
            try:
                replacement = unicodedata.normalize('NFKD', ch)[0]
                if replacement in whitelist:
                    table[ord(ch)] = replacement
                else:
                    table[ord(ch)] = u'_'
            except:
                table[ord(ch)] = u'_'

    if replace_spaces:
        return s.translate(table).replace(u'_', u'').replace(' ', '_')
    else:
        return s.translate(table).replace(u'_', u'')
50401
function parse_RgceArea(blob, length, opts) {
	var w = 2;
	if(opts) {
		if(opts.biff >= 2 && opts.biff <= 5) return parse_RgceArea_BIFF2(blob, length, opts);
		else if(opts.biff == 12) w = 4;
	}
	var r=blob.read_shift(w), R=blob.read_shift(w);
	var c=parse_ColRelU(blob, 2);
	var C=parse_ColRelU(blob, 2);
	return { s:{r:r, c:c[0], cRel:c[1], rRel:c[2]}, e:{r:R, c:C[0], cRel:C[1], rRel:C[2]} };
}
586078
def main(input_bed, output_file, output_features=False, genome=None,
         only_canonical=False, short=False, extended=False, high_confidence=False,
         ambiguities_method=False, coding_only=False, collapse_exons=False, work_dir=False, is_debug=False):
    
    logger.init(is_debug_=is_debug)

    if not genome:
        raise click.BadParameter('Error: please, specify genome build name with -g (e.g. `-g hg19`)', param='genome')

    if short:
        if extended:        raise click.BadParameter('--short and --extended can\'t be set both', param='extended')
        if output_features: raise click.BadParameter('--short and --output-features can\'t be set both', param='output_features')
    elif output_features or extended:
        extended = True
        short    = False

    if not verify_file(input_bed):
        click.BadParameter(f'Usage: {__file__} Input_BED_file -g hg19 -o Annotated_BED_file [options]', param='input_bed')
    input_bed = verify_file(input_bed, is_critical=True, description=f'Input BED file for {__file__}')

    if work_dir:
        work_dir = join(adjust_path(work_dir), os.path.splitext(basename(input_bed))[0])
        safe_mkdir(work_dir)
        info(f'Created work directory {work_dir}')
    else:
        work_dir = mkdtemp('bed_annotate')
        debug('Created temporary work directory {work_dir}')

    input_bed = clean_bed(input_bed, work_dir)
    input_bed = verify_bed(input_bed, is_critical=True, description=f'Input BED file for {__file__} after cleaning')

    output_file = adjust_path(output_file)

    output_file = annotate(
        input_bed, output_file, work_dir, genome=genome,
        only_canonical=only_canonical, short=short, extended=extended,
        high_confidence=high_confidence, collapse_exons=collapse_exons,
        output_features=output_features,
        ambiguities_method=ambiguities_method, coding_only=coding_only,
        is_debug=is_debug)

    if not work_dir:
        debug(f'Removing work directory {work_dir}')
        shutil.rmtree(work_dir)

    info(f'Done, saved to {output_file}')
913369
def run_sparser(pmid_list, tmp_dir, num_cores, start_index, end_index,
                force_read, force_fulltext, cleanup=True, verbose=True):
    ''
    reader_version = sparser.get_version()
    _, _, _, pmids_read, pmids_unread, _ =\
        get_content_to_read(
            pmid_list, start_index, end_index, tmp_dir, num_cores,
            force_fulltext, force_read, 'sparser', reader_version
            )

    logger.info('Adjusting num cores to length of pmid_list.')
    num_cores = min(len(pmid_list), num_cores)
    logger.info('Adjusted...')
    if num_cores is 1:
        stmts = get_stmts(pmids_unread, cleanup=cleanup)
        stmts.update({pmid: get_stmts_from_cache(pmid)[pmid]
                      for pmid in pmids_read.keys()})
    elif num_cores > 1:
        logger.info("Starting a pool with %d cores." % num_cores)
        pool = mp.Pool(num_cores)
        pmids_to_read = list(pmids_unread.keys())
        N = len(pmids_unread)
        dn = int(N/num_cores)
        logger.info("Breaking pmids into batches.")
        batches = []
        for i in range(num_cores):
            batches.append({
                k: pmids_unread[k]
                for k in pmids_to_read[i*dn:min((i+1)*dn, N)]
                })
        get_stmts_func = functools.partial(
            get_stmts,
            cleanup=cleanup,
            sparser_version=reader_version
            )
        logger.info("Mapping get_stmts onto pool.")
        unread_res = pool.map(get_stmts_func, batches)
        logger.info('len(unread_res)=%d' % len(unread_res))
        read_res = pool.map(get_stmts_from_cache, pmids_read.keys())
        logger.info('len(read_res)=%d' % len(read_res))
        pool.close()
        logger.info('Multiprocessing pool closed.')
        pool.join()
        logger.info('Multiprocessing pool joined.')
        stmts = {
            pmid: stmt_list for res_dict in unread_res + read_res
            for pmid, stmt_list in res_dict.items()
            }
        logger.info('len(stmts)=%d' % len(stmts))

    return (stmts, pmids_unread)
728693
def nl2br(self, text):
        
        if isinstance(text, bytes):
            return text.replace(b'\n', b'<br/>\n')
        else:
            return text.replace('\n', '<br/>\n')
103937
def _getReadAlignments(
            self, reference, start, end, readGroupSet, readGroup):
        
        # TODO If reference is None, return against all references,
        # including unmapped reads.
        samFile = self.getFileHandle(self._dataUrl)
        referenceName = reference.getLocalId().encode()
        # TODO deal with errors from htslib
        start, end = self.sanitizeAlignmentFileFetch(start, end)
        readAlignments = samFile.fetch(referenceName, start, end)
        for readAlignment in readAlignments:
            tags = dict(readAlignment.tags)
            if readGroup is None:
                if 'RG' in tags:
                    alignmentReadGroupLocalId = tags['RG']
                    readGroupCompoundId = datamodel.ReadGroupCompoundId(
                        readGroupSet.getCompoundId(),
                        str(alignmentReadGroupLocalId))
                yield self.convertReadAlignment(
                    readAlignment, readGroupSet, str(readGroupCompoundId))
            else:
                if self._filterReads:
                    if 'RG' in tags and tags['RG'] == self._localId:
                        yield self.convertReadAlignment(
                            readAlignment, readGroupSet,
                            str(readGroup.getCompoundId()))
                else:
                    yield self.convertReadAlignment(
                        readAlignment, readGroupSet,
                        str(readGroup.getCompoundId()))
452393
def novo2(args):
    
    p = OptionParser(novo2.__doc__)
    p.set_fastq_names()
    p.set_align(pctid=95)
    opts, args = p.parse_args(args)

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

    trimmed, pf = args
    pctid = opts.pctid
    reads, samples = scan_read_files(trimmed, opts.names)

    # Set up directory structure
    clustdir = "uclust"
    acdir = "allele_counts"
    for d in (clustdir, acdir):
        mkdir(d)

    mm = MakeManager()
    clustfiles = []
    # Step 0 - clustering within sample
    for s in samples:
        flist = [x for x in reads if op.basename(x).split(".")[0] == s]
        outfile = s + ".P{0}.clustS".format(pctid)
        outfile = op.join(clustdir, outfile)
        cmd = "python -m jcvi.apps.uclust cluster --cpus=8"
        cmd += " {0} {1}".format(s, " ".join(flist))
        cmd += " --outdir={0}".format(clustdir)
        cmd += " --pctid={0}".format(pctid)
        mm.add(flist, outfile, cmd)
        clustfiles.append(outfile)

    # Step 1 - make consensus within sample
    allcons = []
    for s, clustfile in zip(samples, clustfiles):
        outfile = s + ".P{0}.consensus".format(pctid)
        outfile = op.join(clustdir, outfile)
        cmd = "python -m jcvi.apps.uclust consensus"
        cmd += " {0}".format(clustfile)
        mm.add(clustfile, outfile, cmd)
        allcons.append(outfile)

    # Step 2 - clustering across samples
    clustSfile = pf + ".P{0}.clustS".format(pctid)
    cmd = "python -m jcvi.apps.uclust mcluster {0}".format(" ".join(allcons))
    cmd += " --prefix={0}".format(pf)
    mm.add(allcons, clustSfile, cmd)

    # Step 3 - make consensus across samples
    locifile = pf + ".P{0}.loci".format(pctid)
    cmd = "python -m jcvi.apps.uclust mconsensus {0}".format(" ".join(allcons))
    cmd += " --prefix={0}".format(pf)
    mm.add(allcons + [clustSfile], locifile, cmd)

    mm.write()
1058206
function parse_RgceArea_BIFF2(blob/*::, length, opts*/) {
	var r=parse_ColRelU(blob, 2), R=parse_ColRelU(blob, 2);
	var c=blob.read_shift(1);
	var C=blob.read_shift(1);
	return { s:{r:r[0], c:c, cRel:r[1], rRel:r[2]}, e:{r:R[0], c:C, cRel:R[1], rRel:R[2]} };
}