Retrieve gVCF blocks within our current evaluation region.

    gvcfgenotyper does not support calling larger regions with individual
    coverage blocks, so we split our big region into potentially multiple.
978610
def _find_gvcf_blocks(vcf_file, region, tmp_dir):
    
    region_file = os.path.join(tmp_dir, "cur_region.bed")
    with open(region_file, "w") as out_handle:
        chrom, coords = region.split(":")
        start, end = coords.split("-")
        out_handle.write("\t".join([chrom, start, end]) + "\n")
    final_file = os.path.join(tmp_dir, "split_regions.bed")
    cmd = "gvcf_regions.py {vcf_file} | bedtools intersect -a - -b {region_file} > {final_file}"
    do.run(cmd.format(**locals()))
    regions = []
    with open(final_file) as in_handle:
        for line in in_handle:
            chrom, start, end = line.strip().split("\t")
            regions.append("%s:%s-%s" % (chrom, start, end))
    return regions
820509
def _call_variants(example_dir, region_bed, data, out_file):
    
    tf_out_file = "%s-tfrecord.gz" % utils.splitext_plus(out_file)[0]
    if not utils.file_exists(tf_out_file):
        with file_transaction(data, tf_out_file) as tx_out_file:
            model = "wes" if strelka2.coverage_interval_from_bed(region_bed) == "targeted" else "wgs"
            cmd = ["dv_call_variants.py", "--cores", dd.get_num_cores(data),
                   "--outfile", tx_out_file, "--examples", example_dir,
                   "--sample", dd.get_sample_name(data), "--model", model]
            do.run(cmd, "DeepVariant call_variants %s" % dd.get_sample_name(data))
    return tf_out_file
379744
def _get_full_block(grouped_dicoms):
    
    # For each slice / mosaic create a data volume block
    data_blocks = []
    for index in range(0, len(grouped_dicoms)):
        logger.info('Creating block %s of %s' % (index + 1, len(grouped_dicoms)))
        data_blocks.append(_timepoint_to_block(grouped_dicoms[index]))

    # Add the data_blocks together to one 4d block
    size_x = numpy.shape(data_blocks[0])[0]
    size_y = numpy.shape(data_blocks[0])[1]
    size_z = numpy.shape(data_blocks[0])[2]
    size_t = len(data_blocks)
    full_block = numpy.zeros((size_x, size_y, size_z, size_t), dtype=data_blocks[0].dtype)
    for index in range(0, size_t):
        if full_block[:, :, :, index].shape != data_blocks[index].shape:
            logger.warning('Missing slices (slice count mismatch between timepoint %s and %s)' % (index - 1, index))
            logger.warning('---------------------------------------------------------')
            logger.warning(full_block[:, :, :, index].shape)
            logger.warning(data_blocks[index].shape)
            logger.warning('---------------------------------------------------------')
            raise ConversionError("MISSING_DICOM_FILES")
        full_block[:, :, :, index] = data_blocks[index]

    return full_block
533524
def isolate_region(sequences, start, end, gap_char='-'):
    
    # Check arguments
    if end <= start:
        raise ValueError("start of slice must precede end ({0} !> {1})".format(
            end, start))

    for sequence in sequences:
        seq = sequence.seq
        start_gap = gap_char * start
        end_gap = gap_char * (len(seq) - end)
        seq = Seq(start_gap + str(seq[start:end]) + end_gap,
                alphabet=seq.alphabet)
        sequence.seq = seq
        yield sequence
902027
def variantcall_batch_region(items):
    
    items = [utils.to_single_data(x) for x in items]
    align_bams = [dd.get_align_bam(x) for x in items]
    variantcaller = _get_batch_variantcaller(items)
    region_blocks = list(set([tuple(x.get("region_block")) for x in items if "region_block" in x]))
    assert len(region_blocks) == 1, region_blocks
    region_block = region_blocks[0]
    # Pre-called input variant files
    if not variantcaller and all(d.get("vrn_file") for d in items):
        return {"vrn_file_region": None, "region_block": region_block}
    caller_fn = get_variantcallers()[variantcaller]
    assoc_files = tz.get_in(("genome_resources", "variation"), items[0], {})
    region = _region_to_coords(region_block[0])
    chrom, start, end = region
    region_str = "_".join(str(x) for x in region)
    batch_name = _get_batch_name(items)
    out_file = os.path.join(dd.get_work_dir(items[0]), variantcaller, chrom,
                            "%s-%s-block.vcf.gz" % (batch_name, region_str))
    utils.safe_makedir(os.path.dirname(out_file))
    with pshared.bedtools_tmpdir(items[0]):
        if variantcaller in SUPPORT_MULTICORE:
            call_file = caller_fn(align_bams, items, dd.get_ref_file(items[0]), assoc_files,
                                  [_region_to_coords(r) for r in region_block], out_file)
        else:
            call_file = _run_variantcall_batch_multicore(items, region_block, out_file)
    return {"vrn_file_region": call_file, "region_block": region_block}
267826
def blocks(self):
        

        for block_addr, block in self._local_blocks.items():
            try:
                yield self._get_block(block_addr, size=block.size,
                                      byte_string=block.bytestr if isinstance(block, BlockNode) else None)
            except (SimEngineError, SimMemoryError):
                pass
221338
def _get_project_vcf(x, suffix=""):
    
    vcfs = _get_variant_file(x, ("population", "vcf"), suffix=suffix)
    if not vcfs:
        vcfs = _get_variant_file(x, ("vrn_file_batch", ), suffix=suffix, ignore_do_upload=True)
        if not vcfs and x.get("variantcaller") == "ensemble":
            vcfs = _get_variant_file(x, ("vrn_file", ), suffix=suffix)
    return vcfs
349849
private void processParityHarPartBlock(FileSystem dfs, Path partFile,
      long blockOffset,
      FileStatus partFileStat,
      HarIndex harIndex,
      File localBlockFile,
      Progressable progress)
  throws IOException {
    String partName = partFile.toUri().getPath(); // Temporarily.
    partName = partName.substring(1 + partName.lastIndexOf(Path.SEPARATOR));

    OutputStream out = new FileOutputStream(localBlockFile);

    try {
      // A HAR part file block could map to several parity files. We need to
      // use all of them to recover this block.
      final long blockEnd = Math.min(blockOffset + 
          partFileStat.getBlockSize(),
          partFileStat.getLen());
      for (long offset = blockOffset; offset < blockEnd; ) {
        HarIndex.IndexEntry entry = harIndex.findEntry(partName, offset);
        if (entry == null) {
          String msg = "Lost index file has no matching index entry for " +
          partName + ":" + offset;
          LOG.warn(msg);
          throw new IOException(msg);
        }
        Path parityFile = new Path(entry.fileName);
        Encoder encoder = null;
        for (Codec codec : Codec.getCodecs()) {
          if (isParityFile(parityFile, codec)) {
            encoder = new Encoder(getConf(), codec);
          }
        }
        if (encoder == null) {
          String msg = "Could not figure out codec correctly for " + parityFile;
          LOG.warn(msg);
          throw new IOException(msg);
        }
        Path srcFile = RaidUtils.sourcePathFromParityPath(parityFile, dfs);
        if (null == srcFile) {
          String msg = "Can not find the source path for parity file: " +
                          parityFile;
          LOG.warn(msg);
          throw new IOException(msg);
        }
        FileStatus srcStat = dfs.getFileStatus(srcFile);
        if (srcStat.getModificationTime() != entry.mtime) {
          String msg = "Modification times of " + parityFile + " and " +
          srcFile + " do not match.";
          LOG.warn(msg);
          throw new IOException(msg);
        }
        long lostOffsetInParity = offset - entry.startOffset;
        LOG.info(partFile + ":" + offset + " maps to " +
            parityFile + ":" + lostOffsetInParity +
            " and will be recovered from " + srcFile);
        encoder.recoverParityBlockToStream(dfs, srcStat,
            srcStat.getBlockSize(), parityFile,
            lostOffsetInParity, out, progress);
        // Finished recovery of one parity block. Since a parity block has the
        // same size as a source block, we can move offset by source block 
        // size.
        offset += srcStat.getBlockSize();
        LOG.info("Recovered " + srcStat.getBlockSize() + " part file bytes ");
        if (offset > blockEnd) {
          String msg =
            "Recovered block spills across part file blocks. Cannot continue";
          throw new IOException(msg);
        }
        progress.progress();
      }
    } finally {
      out.close();
    }
  }
1244422
def get_region(self, x,z):
        """"""
        if (x,z) not in self.regions:
            if (x,z) in self.regionfiles:
                self.regions[(x,z)] = region.RegionFile(self.regionfiles[(x,z)])
            else:
                # Return an empty RegionFile object
                # TODO: this does not yet allow for saving of the region file
                # TODO: this currently fails with a ValueError!
                # TODO: generate the correct name, and create the file
                # and add the fie to self.regionfiles
                self.regions[(x,z)] = region.RegionFile()
            self.regions[(x,z)].loc = Location(x=x,z=z)
        return self.regions[(x,z)]
598593
def _get_interval(variant_regions, region, out_file, items):
    
    target = shared.subset_variant_regions(variant_regions, region, out_file, items)
    if target:
        if isinstance(target, six.string_types) and os.path.isfile(target):
            return "--interval %s" % target
        else:
            return "--interval %s" % bamprep.region_to_gatk(target)
    else:
        return ""
478206
def _variant_vc(checkpoints):
    
    if not checkpoints.get("vc"):
        return [], []
    vc_wf = [s("get_parallel_regions", "batch-split",
               [["batch_rec"]],
               [cwlout(["region_block"], {"type": "array", "items": "string"})],
               "bcbio-vc",
               disk={"files": 2.0}, cores=1),
             s("variantcall_batch_region", "batch-parallel",
               [["batch_rec"], ["region_block"]],
               [cwlout(["vrn_file_region"], ["File", "null"], [".tbi"]),
                cwlout(["region_block"], {"type": "array", "items": "string"})],
               "bcbio-vc", ["bcftools", "bedtools", "freebayes=1.1.0.46",
                            "gatk4", "vqsr_cnn", "deepvariant;env=dv", "sentieon;env=python2",
                            "htslib", "octopus", "picard", "platypus-variant;env=python2", "pythonpy",
                            "samtools", "pysam>=0.13.0", "strelka;env=python2", "vardict", "vardict-java",
                            "varscan", "moreutils", "vcfanno", "vcflib", "vt", "r=3.5.1", "r-base",
                            "perl"],
               disk={"files": 2.0}),
             s("concat_batch_variantcalls", "batch-merge",
               [["batch_rec"], ["region_block"], ["vrn_file_region"]],
               [cwlout(["vrn_file"], "File", [".tbi"])],
               "bcbio-vc", ["bcftools", "htslib", "gatk4"],
               disk={"files": 1.5}, cores=1)]
    if not checkpoints.get("jointvc"):
        vc_wf += [s("postprocess_variants", "batch-single",
                    [["batch_rec"], ["vrn_file"]],
                    [cwlout(["vrn_file"], "File", [".tbi"])],
                    "bcbio-vc", ["snpeff=4.3.1t"], disk={"files": 0.5})]
    vc_rec_exclude = [["align_bam"]]
    if not checkpoints.get("jointvc"):
        vc_rec_exclude.append(["genome_resources", "variation"])
    vc_wf += [s("compare_to_rm", "batch-single",
                [["batch_rec"], ["vrn_file"]],
                [cwlout("vc_rec", "record",
                        fields=[cwlout(["batch_samples"], ["null", {"type": "array", "items": "string"}]),
                                cwlout(["validate", "summary"], ["File", "null"]),
                                cwlout(["validate", "tp"], ["File", "null"], [".tbi"]),
                                cwlout(["validate", "fp"], ["File", "null"], [".tbi"]),
                                cwlout(["validate", "fn"], ["File", "null"], [".tbi"]),
                                cwlout("inherit", exclude=vc_rec_exclude)])],
                "bcbio-vc", ["bcftools", "bedtools", "pythonpy", "gvcf-regions;env=python2",
                             "htslib", "rtg-tools", "vcfanno"],
                disk={"files": 1.5})]
    batch_in = [["analysis"], ["genome_build"], ["align_bam"], ["vrn_file"],
                ["metadata", "batch"], ["metadata", "phenotype"],
                ["config", "algorithm", "callable_regions"], ["regions", "sample_callable"],
                ["config", "algorithm", "variantcaller"],
                ["config", "algorithm", "ensemble"],
                ["config", "algorithm", "vcfanno"],
                ["config", "algorithm", "coverage_interval"],
                ["config", "algorithm", "effects"],
                ["config", "algorithm", "min_allele_fraction"],
                ["config", "algorithm", "exclude_regions"],
                ["config", "algorithm", "variant_regions"],
                ["config", "algorithm", "variant_regions_merged"],
                ["config", "algorithm", "validate"], ["config", "algorithm", "validate_regions"],
                ["config", "algorithm", "tools_on"],
                ["config", "algorithm", "tools_off"],
                ["reference", "fasta", "base"],
                ["reference", "rtg"], ["reference", "genome_context"],
                ["genome_resources", "variation", "clinvar"],
                ["genome_resources", "variation", "cosmic"], ["genome_resources", "variation", "dbsnp"],
                ["genome_resources", "variation", "esp"], ["genome_resources", "variation", "exac"],
                ["genome_resources", "variation", "gnomad_exome"],
                ["genome_resources", "variation", "1000g"],
                ["genome_resources", "variation", "lcr"], ["genome_resources", "variation", "polyx"],
                ["genome_resources", "variation", "encode_blacklist"],
                ["genome_resources", "aliases", "ensembl"], ["genome_resources", "aliases", "human"],
                ["genome_resources", "aliases", "snpeff"], ["reference", "snpeff", "genome_build"]]
    if checkpoints.get("umi"):
        batch_in.append(["config", "algorithm", "umi_type"])
    if checkpoints.get("rnaseq"):
        batch_in += [["genome_resources", "variation", "editing"]]
    else:
        batch_in += [["genome_resources", "variation", "train_hapmap"],
                     ["genome_resources", "variation", "train_indels"]]
    vc = [s("batch_for_variantcall", "multi-batch", batch_in,
            [cwlout("batch_rec", "record",
                    fields=[cwlout(["config", "algorithm", "variantcaller_order"], "int"),
                            cwlout("inherit")])],
            "bcbio-vc",
            disk={"files": 2.0}, cores=1,
            unlist=[["config", "algorithm", "variantcaller"]], no_files=True),
          w("variantcall", "multi-parallel", vc_wf,
            [["region"], ["region_block"], ["vrn_file_region"], ["vrn_file"], ["validate", "summary"]])]
    if checkpoints.get("jointvc"):
        vc += _variant_jointvc()
    if checkpoints.get("ensemble"):
        vc += _variant_ensemble(checkpoints)
    summarize_in = [["jointvc_rec" if checkpoints.get("jointvc") else "vc_rec"]]
    if checkpoints.get("ensemble"):
        summarize_in += [["ensemble_rec"]]
    vc += [s("summarize_vc", "multi-combined", summarize_in,
             [cwlout(["variants", "calls"], {"type": "array", "items": ["File", "null"]}),
              cwlout(["variants", "gvcf"], ["null", {"type": "array", "items": ["File", "null"]}]),
              cwlout(["variants", "samples"], {"type": "array", "items": {"type": "array",
                                                                          "items": ["File", "null"]}}),
              cwlout(["validate", "grading_summary"], ["File", "null"]),
              cwlout(["validate", "grading_plots"], {"type": "array", "items": ["File", "null"]})],
             "bcbio-vc",
             disk={"files": 2.0}, cores=1)]
    return vc, [["validate", "grading_summary"], ["variants", "calls"], ["variants", "gvcf"]]
553758
private ByteBuffer fetchBlockByteRangeScatterGather(LocatedBlock block, 
                      long start, long len) throws IOException {
    //
    // Connect to best DataNode for desired Block, with potential offset
    //
    Socket dn = null;

    // cached block locations may have been updated by chooseDatNode()
    // or fetchBlockAt(). Always get the latest list of locations before the
    // start of the loop.
    BlockSeekContext seekContext = new BlockSeekContext(getBlockAt(
        block.getStartOffset(), false, true));
    while (true) {
      DNAddrPair retval = chooseDataNode(seekContext);
      DatanodeInfo chosenNode = retval.info;
      InetSocketAddress targetAddr = retval.addr;
      ByteBuffer result = null;
      BlockReaderLocalBase localReader = null;
      BlockReaderAccelerator remoteReader = null;

       try {
         if (DFSClient.LOG.isDebugEnabled()) {
           DFSClient.LOG.debug("fetchBlockByteRangeScatterGather " +
                    " localhst " + dfsClient.localHost +
                    " targetAddr " + targetAddr);
         }
         
         // first try reading the block locally.
         if (dfsClient.shortCircuitLocalReads && 
             NetUtils.isLocalAddressWithCaching(targetAddr.getAddress())) {
           localReader = BlockReaderLocalBase.newBlockReader(dfsClient.conf, src,
                                                dfsClient.namespaceId, block.getBlock(),
                                                chosenNode,
                                                start,
                                                len,
                                                dfsClient.metrics,
                                                verifyChecksum,
                                                this.clearOsBuffer,
                                                false);
           if (localReader != null) {
             localReader.setReadLocal(true);
             localReader.setFsStats(dfsClient.stats);
             result = localReader.readAll();
           } else if (!dfsClient.shortcircuitDisableWhenFail) {
            throw new IOException(
                "Short circuit local read not supported for this scase");
           }
         }
         if (localReader == null) {

           // go to the datanode
           dn = dfsClient.socketFactory.createSocket();
           NetUtils.connect(dn, targetAddr, dfsClient.socketTimeout,
               dfsClient.ipTosValue);
           dn.setSoTimeout(dfsClient.socketTimeout);
           remoteReader = new BlockReaderAccelerator(dfsClient.conf,
                                          targetAddr,
                                          chosenNode,
                                          dfsClient.getDataTransferProtocolVersion(),
                                          dfsClient.namespaceId, dfsClient.clientName,
                                          dn, src,
                                          block,
                                          start, len,
                                          verifyChecksum, dfsClient.metrics);
           result = remoteReader.readAll();
          }
          if (result.remaining() != len) {
            throw new IOException("truncated return from reader.read(): " +
                                "expected " + len + ", got " + 
                                  result.remaining());
          }
          if (NetUtils.isLocalAddressWithCaching(targetAddr.getAddress())) {
            dfsClient.stats.incrementLocalBytesRead(len);
            dfsClient.stats.incrementRackLocalBytesRead(len);
          } else if (dfsClient.isInLocalRack(targetAddr)) {
            dfsClient.stats.incrementRackLocalBytesRead(len);
          }

          return result;
      } catch (ChecksumException e) {
        DFSClient.LOG.warn("fetchBlockByteRangeScatterGather(). Got a checksum exception for " +
                 src + " at " + block.getBlock() + ":" +
                 e.getPos() + " from " + chosenNode.getName());
        dfsClient.reportChecksumFailure(src, block.getBlock(), chosenNode);
      } catch (IOException e) {
        DFSClient.LOG.warn("Failed to connect to " + targetAddr +
                 " for file " + src +
                 " for block " + block.getBlock().getBlockId() + ":"  +
                 StringUtils.stringifyException(e));
      } finally {
        IOUtils.closeStream(localReader);
        IOUtils.closeStream(remoteReader);
        IOUtils.closeSocket(dn);
      }
       dfsClient.incReadExpCntToStats();
      // Put chosen node into dead list, continue
      addToDeadNodes(chosenNode);
    }
  }
92098
def _get_vcf_breakends(hydra_file, genome_2bit, options=None):
    
    if options is None: options = {}
    for features in group_hydra_breakends(hydra_parser(hydra_file, options)):
        if len(features) == 1 and is_deletion(features[0], options):
            yield build_vcf_deletion(features[0], genome_2bit)
        elif len(features) == 1 and is_tandem_dup(features[0], options):
            yield build_tandem_deletion(features[0], genome_2bit)
        elif len(features) == 2 and is_inversion(*features):
            yield build_vcf_inversion(features[0], features[1], genome_2bit)
        elif len(features) == 2 and is_translocation(*features):
            info = get_translocation_info(features[0], features[1])
            for feature in features:
                for brend in build_vcf_parts(feature, genome_2bit, info):
                    yield brend
        else:
            for feature in features:
                for brend in build_vcf_parts(feature, genome_2bit):
                    yield brend
1098679
def _get_variant_regions(items):
    
    return list(filter(lambda x: x is not None,
                       [tz.get_in(("config", "algorithm", "variant_regions"), data)
                        for data in items
                        if tz.get_in(["config", "algorithm", "coverage_interval"], data) != "genome"]))
1280576
def subset_variant_regions(variant_regions, region, out_file, items=None, do_merge=True, data=None):
    
    if region is None:
        return variant_regions
    elif variant_regions is None:
        return region
    elif not isinstance(region, (list, tuple)) and region.find(":") > 0:
        raise ValueError("Partial chromosome regions not supported")
    else:
        merge_text = "-unmerged" if not do_merge else ""
        subset_file = "{0}".format(utils.splitext_plus(out_file)[0])
        subset_file += "%s-regions.bed" % (merge_text)
        if not os.path.exists(subset_file):
            data = items[0] if items else data
            with file_transaction(data, subset_file) as tx_subset_file:
                if isinstance(region, (list, tuple)):
                    _subset_bed_by_region(variant_regions, tx_subset_file, to_multiregion(region),
                                          dd.get_ref_file(data), do_merge=do_merge)
                else:
                    _rewrite_bed_with_chrom(variant_regions, tx_subset_file, region)
        if os.path.getsize(subset_file) == 0:
            return region
        else:
            return subset_file
633382
def _subset_regions(region, base_file, items):
    
    variant_regions = bedutils.population_variant_regions(items, merged=True)
    target = pshared.subset_variant_regions(variant_regions, region, base_file, items)
    if isinstance(target, six.string_types) and os.path.isfile(target):
        return target
    else:
        return bamprep.region_to_gatk(target)
698239
def _run_gvcfgenotyper(data, region, vrn_files, out_file):
    
    if not utils.file_exists(out_file):
        with file_transaction(data, out_file) as tx_out_file:
            input_file = "%s-inputs.txt" % utils.splitext_plus(tx_out_file)[0]
            with open(input_file, "w") as out_handle:
                out_handle.write("%s\n" % "\n".join(vrn_files))
            cmd = ["gvcfgenotyper", "-f", dd.get_ref_file(data), "-l", input_file,
                   "-r", region, "-O", "z", "-o", tx_out_file]
            do.run(cmd, "gvcfgenotyper: %s %s" % (dd.get_sample_name(data), region))
    return out_file