import htsjdk.samtools.util.IntervalTree; //导入依赖的package包/类
private int processExon(IntervalTree<Block> intervalTree, int totalLength, GeneFeature feature,
                        int intronLength, int centerPosition, boolean forward) {
    int currentLength = 0;
    if (EXON_FEATURE_NAME.equalsIgnoreCase(feature.getFeature())) {
        Block exon = createExon(feature, intronLength);

        Iterator<IntervalTree.Node<Block>> nodeIterator = intervalTree.overlappers(exon.getStartIndex(),
                exon.getEndIndex());
        if (nodeIterator.hasNext()) {
            Block merged = exon;
            while (nodeIterator.hasNext()) {
                Block toMerge = nodeIterator.next().getValue();
                merged = mergeExons(merged, toMerge);
                currentLength = currentLength - calculateExonLength(toMerge, centerPosition, forward)
                              + calculateExonLength(merged, centerPosition, forward);
                nodeIterator.remove();
            }

            intervalTree.put(merged.getStartIndex(), merged.getEndIndex(), merged);
        } else {
            currentLength += calculateExonLength(exon, centerPosition, forward);
            intervalTree.put(exon.getStartIndex(), exon.getEndIndex(), exon);
        }
    }
    return totalLength + currentLength;
}
 

import htsjdk.samtools.util.IntervalTree; //导入依赖的package包/类
private void loadExonsForward(int centerPosition, int viewPortSize, Chromosome chromosome,
                              int intronLength, int endFeatures, IntervalTree<Block> intervalTree,
                              AbstractFeatureReader<GeneFeature, LineIterator> featureReader) throws IOException {
    CloseableIterator<GeneFeature> iterator = featureReader.query(chromosome.getName(), centerPosition,
            endFeatures);

    if (!iterator.hasNext()) {
        iterator = featureReader.query(Utils.changeChromosomeName(chromosome.getName()), centerPosition,
                endFeatures);
    }

    int totalLength = 0;
    while (iterator.hasNext() && totalLength < viewPortSize / 2) {
        GeneFeature feature = iterator.next();
        totalLength = processExon(intervalTree, totalLength, feature, intronLength, centerPosition, true);
    }
}
 

import htsjdk.samtools.util.IntervalTree; //导入依赖的package包/类
private IntervalTree<Annotation> retrieveTranscripts(Genomics genomics, StreamVariantsRequest request) {
  Stopwatch stopwatch = Stopwatch.createStarted();
  IntervalTree<Annotation> transcripts = new IntervalTree<>();
  Iterable<Annotation> transcriptIter =
      Paginator.Annotations.create(genomics, ShardBoundary.Requirement.OVERLAPS).search(
          new SearchAnnotationsRequest()
            .setAnnotationSetIds(transcriptSetIds)
            .setReferenceName(canonicalizeRefName(request.getReferenceName()))
            .setStart(request.getStart())
            .setEnd(request.getEnd()));
  for (Annotation annotation : transcriptIter) {
    transcripts.put(annotation.getStart().intValue(), annotation.getEnd().intValue(), annotation);
  }
  LOG.info(String.format("read %d transcripts in %s (%.2f / s)", transcripts.size(),
      stopwatch, (double)transcripts.size() / stopwatch.elapsed(TimeUnit.SECONDS)));
  return transcripts;
}
 

import htsjdk.samtools.util.IntervalTree; //导入依赖的package包/类
/**
 * Loads exon features in the requested viewport
 *
 * @param geneFileId an ID of GeneFile to load exons from
 * @param chromosomeId an ID of chromosome to load exons from
 * @param centerPosition a position of center of the screen
 * @param viewPortSize a size of the screen
 * @return a List of exon intervals. Overlapping exons are joined together
 * @throws IOException
 */
public List<Block> loadExonsInViewPort(long geneFileId, long chromosomeId, int centerPosition, int viewPortSize,
                                      int intronLength)
        throws IOException {
    double time1 = Utils.getSystemTimeMilliseconds();
    final GeneFile geneFile = geneFileManager.loadGeneFile(geneFileId);
    final Chromosome chromosome = referenceGenomeManager.loadChromosome(chromosomeId);
    Assert.notNull(chromosome, getMessage(MessagesConstants.ERROR_CHROMOSOME_ID_NOT_FOUND));

    final Map<String, Pair<Integer, Integer>> metaMap = fileManager.loadIndexMetadata(geneFile);
    Pair<Integer, Integer> bounds = metaMap.get(chromosome.getName());
    if (bounds == null) {
        bounds = metaMap.get(Utils.changeChromosomeName(chromosome.getName()));
    }

    Assert.notNull(bounds, getMessage(MessageCode.NO_SUCH_CHROMOSOME));
    int end = bounds.getRight();
    int start = bounds.getLeft();

    IntervalTree<Block> intervalTree = new IntervalTree<>();

    try (AbstractFeatureReader<GeneFeature, LineIterator> featureReader =
                 fileManager.makeGeneReader(geneFile, GeneFileType.ORIGINAL)) {
        loadExonsForward(centerPosition, viewPortSize, chromosome, intronLength, end, intervalTree,
                featureReader);
        loadExonsBackwards(centerPosition, viewPortSize, chromosome, intronLength, start, intervalTree,
                featureReader);
    }

    List<Block> exons = new ArrayList<>();
    for (IntervalTree.Node<Block> node : intervalTree) {
        exons.add(node.getValue());
    }
    double time2 = Utils.getSystemTimeMilliseconds();
    LOGGER.info(getMessage(MessagesConstants.DEBUG_GENE_EXONS_LOAD, time2 - time1));
    return exons;
}
 

import htsjdk.samtools.util.IntervalTree; //导入依赖的package包/类
private void testOverlapping(List<Block> exons) {
    IntervalTree<Block> intervalTree = new IntervalTree<>();
    for (Block exon : exons) {
        Iterator<IntervalTree.Node<Block>> nodeIterator = intervalTree.overlappers(exon.getStartIndex(),
                                                                                   exon.getEndIndex());
        Assert.assertFalse("Should be no overlapping exons", nodeIterator.hasNext());
        intervalTree.put(exon.getStartIndex(), exon.getEndIndex(), exon);
    }
}
 

import htsjdk.samtools.util.IntervalTree; //导入依赖的package包/类
public IntervalTree getSCNVintervalTreeWithExtendBase(ArrayList<VarscanCNVFormat> SCNVlist, int extendBaseNumber){
     IntervalTree CNVIntervalTree = new IntervalTree();
     VarscanCNVFormat vscnv=null;
     for (int i = 0; i < SCNVlist.size(); i++) {
      vscnv=SCNVlist.get(i);
      CNVIntervalTree.put(Integer.parseInt(vscnv.getChr_start())-extendBaseNumber, Integer.parseInt(vscnv.getChr_stop())+extendBaseNumber, vscnv);
    }
     return CNVIntervalTree;
}
 

import htsjdk.samtools.util.IntervalTree; //导入依赖的package包/类
public IntervalTree getSCNVintervalTree(ArrayList<VarscanCNVFormat> SCNVlist){
     IntervalTree CNVIntervalTree = new IntervalTree();
     VarscanCNVFormat vscnv=null;
     for (int i = 0; i < SCNVlist.size(); i++) {
      vscnv=SCNVlist.get(i);
      CNVIntervalTree.put(Integer.parseInt(vscnv.getChr_start()), Integer.parseInt(vscnv.getChr_stop()), vscnv);
    }
     return CNVIntervalTree;
}
 

import htsjdk.samtools.util.IntervalTree; //导入依赖的package包/类
public IntervalTree getInDelIntervalTreeWithExtendBase(ArrayList<VCFrecord> Indellist, int extendBaseNumber){
     IntervalTree CNVIntervalTree = new IntervalTree();
     VCFrecord vcf=null;
     int start;
     int end;
     for (int i = 0; i < Indellist.size(); i++) {
      vcf=Indellist.get(i);
      //start=Integer.parseInt(vcf.get)
     // CNVIntervalTree.put(-extendBaseNumber, Integer.parseInt(vscnv.getChr_stop())+extendBaseNumber, vcf);
    }
     return CNVIntervalTree;
}
 

import htsjdk.samtools.util.IntervalTree; //导入依赖的package包/类
public IntervalTree getBEDIntervalTreeWithExtendBase(ArrayList<BED> SCNVlist, int extendBaseNumber){
     IntervalTree CNVIntervalTree = new IntervalTree();
     BED bed=null;
     for (int i = 0; i < SCNVlist.size(); i++) {
      bed=SCNVlist.get(i);
      CNVIntervalTree.put(bed.getChromStart()-extendBaseNumber, bed.getChromEnd()+extendBaseNumber, bed);
    }
     return CNVIntervalTree;
}
 

import htsjdk.samtools.util.IntervalTree; //导入依赖的package包/类
private void loadExonsBackwards(int centerPosition, int viewPortSize, Chromosome chromosome,
                                int intronLength, int featuresStart,
                                final IntervalTree<Block> intervalTree,
                                final AbstractFeatureReader<GeneFeature, LineIterator> featureReader)
        throws IOException {
    int totalLength = 0;
    // check if some of exons, got by forward lookup are good for backwards
    Iterator<IntervalTree.Node<Block>> nodeIterator = intervalTree.overlappers(featuresStart, centerPosition);
    while (nodeIterator.hasNext()) {
        Block exon = nodeIterator.next().getValue();
        totalLength += calculateExonLength(exon, centerPosition, false);
    }

    int i = 0;
    boolean lastChunk = false;

    while (totalLength < viewPortSize / 2) {
        if (lastChunk) {
            break;
        }

        int firstIndex = centerPosition - EXON_SEARCH_CHUNK_SIZE * (i + 1);
        final int lastIndex = centerPosition - 1 - EXON_SEARCH_CHUNK_SIZE * i;
        if (firstIndex < featuresStart) {
            firstIndex = featuresStart;
            lastChunk = true; // this is the last chunk to be traversed
        }

        CloseableIterator<GeneFeature> iterator = Utils.query(featureReader, chromosome, firstIndex, lastIndex);
            // instead traversing the whole file, read it by small chunks, 100000 bps
            // long. Hopefully, the desired window will be covered by first/second chunk

        if (iterator.hasNext()) {
            List<GeneFeature> featuresChunk = iterator.toList();
            ListIterator<GeneFeature> listIterator = featuresChunk.listIterator(featuresChunk.size() - 1);

            while (listIterator.hasPrevious() && totalLength < viewPortSize / 2) {
                GeneFeature feature = listIterator.previous();
                totalLength =
                    processExon(intervalTree, totalLength, feature, intronLength, centerPosition, false);
            }
        }

        i++;
    }
}
 

import htsjdk.samtools.util.IntervalTree; //导入依赖的package包/类
@Override
public IntervalTree<T> debugGetTree(String sequence) {
    return treeMap.debugGetTree(sequence);
}
 

import htsjdk.samtools.util.IntervalTree; //导入依赖的package包/类
@SuppressWarnings({ "rawtypes", "unchecked" })
@Override
public void update() throws InvalidGenomicCoordsException, IOException{
	
	if(this.getyMaxLines() == 0){
		return;
	}
	if( ! this.useSamtools()){
		String chrom= this.getGc().getChrom();
		
		if(! this.loci.containsKey(chrom)){
			this.loci.put(chrom, new HashMap<Integer, Locus>());
		}
		if(! this.zeroDepthIntervals.containsKey(chrom)){
			this.zeroDepthIntervals.put(chrom, new IntervalTree());
		}

		// Check cache is not growing too much
		if(this.loci.get(chrom).keySet().size() > 500000){
			this.loci.get(chrom).clear();
			this.zeroDepthIntervals.get(chrom).clear();
		}
		
		// Find the positions that we haven't visited before:
		List<Integer> missingPos= new ArrayList<Integer>();
		for(int pos= this.getGc().getFrom(); pos <= this.getGc().getTo(); pos++){
			if( ! this.loci.get(chrom).containsKey(pos) &&  
				! this.zeroDepthIntervals.get(chrom).overlappers(pos, pos).hasNext()){
				missingPos.add(pos);
			}
		}
		for(List<Integer> gap : mergePositionsInIntervals(missingPos)){

			int qryFrom= gap.get(0);
			int qryTo= gap.get(1);
			
			SamReader samReader= Utils.getSamReader(this.getWorkFilename());
			List<Boolean> passFilter= this.filterReads(samReader, chrom, qryFrom, qryTo);
			samReader= Utils.getSamReader(this.getWorkFilename());
			
			Iterator<SAMRecord> sam= samReader.query(chrom, qryFrom, qryTo, false);
	
			ListIterator<Boolean> pass = passFilter.listIterator();
			while(sam.hasNext()){
				SAMRecord rec= sam.next();
				if(pass.next()){
					this.add(rec, qryFrom, qryTo, this.loci.get(chrom));
				}
			}
			
			// Now add the loci that have been collected in this last update
			List<Integer> zeroDepthPos= new ArrayList<Integer>();
			for(int pos= qryFrom; pos <= qryTo; pos++){
				if( ! this.loci.get(chrom).containsKey(pos)){
					zeroDepthPos.add(pos);
				}
			}
			List<List<Integer>> zeroDepthGaps = mergePositionsInIntervals(zeroDepthPos);
			for(List<Integer> z : zeroDepthGaps){
				this.zeroDepthIntervals.get(chrom).put(z.get(0), z.get(1), null);
			}
		}
	}
	List<Float> screenScores= this.prepareScreenScores();
	this.setScreenScores(screenScores);
}
 

import htsjdk.samtools.util.IntervalTree; //导入依赖的package包/类
@ProcessElement
public void processElement(
    DoFn<StreamVariantsRequest, KV<String, VariantAnnotation>>.ProcessContext c) throws Exception {
  Genomics genomics = GenomicsFactory.builder().build().fromOfflineAuth(auth);

  StreamVariantsRequest request = StreamVariantsRequest.newBuilder(c.element())
      .addAllCallSetIds(callSetIds)
      .build();
  LOG.info("processing contig " + request);

  Iterator<StreamVariantsResponse> iter = VariantStreamIterator.enforceShardBoundary(auth,
      request, ShardBoundary.Requirement.STRICT, VARIANT_FIELDS);
  if (!iter.hasNext()) {
    LOG.info("region has no variants, skipping");
    return;
  }

  IntervalTree<Annotation> transcripts = retrieveTranscripts(genomics, request);
  ListMultimap<Range<Long>, Annotation> variantAnnotations =
      retrieveVariantAnnotations(genomics, request);

  Stopwatch stopwatch = Stopwatch.createStarted();
  int varCount = 0;
  while (iter.hasNext()) {
    Iterable<Variant> varIter = FluentIterable
        .from(iter.next().getVariantsList())
        .filter(VariantUtils.IS_SNP);
    for (Variant variant : varIter) {
      List<String> alleles = ImmutableList.<String>builder()
          .addAll(variant.getAlternateBasesList())
          .add(variant.getReferenceBases())
          .build();
      Range<Long> pos = Range.closedOpen(variant.getStart(), variant.getEnd());
      for (String allele : alleles) {
        String outKey = Joiner.on(":").join(
            variant.getReferenceName(), variant.getStart(), allele, variant.getId());
        for (Annotation match : variantAnnotations.get(pos)) {
          if (allele.equals(match.getVariant().getAlternateBases())) {
            // Exact match to a known variant annotation; straightforward join.
            c.output(KV.of(outKey, match.getVariant()));
          }
        }

        Iterator<Node<Annotation>> transcriptIter = transcripts.overlappers(
            pos.lowerEndpoint().intValue(), pos.upperEndpoint().intValue() - 1); // Inclusive.
        while (transcriptIter.hasNext()) {
          // Calculate an effect of this allele on the coding region of the given transcript.
          Annotation transcript = transcriptIter.next().getValue();
          VariantEffect effect = AnnotationUtils.determineVariantTranscriptEffect(
              variant.getStart(), allele, transcript,
              getCachedTranscriptBases(genomics, transcript));
          if (effect != null && !VariantEffect.SYNONYMOUS_SNP.equals(effect)) {
            c.output(KV.of(outKey, new VariantAnnotation()
            .setAlternateBases(allele)
            .setType("SNP")
            .setEffect(effect.toString())
            .setGeneId(transcript.getTranscript().getGeneId())
            .setTranscriptIds(ImmutableList.of(transcript.getId()))));
          }
        }
      }
      varCount++;
      if (varCount%1e3 == 0) {
        LOG.info(String.format("read %d variants (%.2f / s)",
            varCount, (double)varCount / stopwatch.elapsed(TimeUnit.SECONDS)));
      }
    }
  }
  LOG.info("finished reading " + varCount + " variants in " + stopwatch);
}
 

import htsjdk.samtools.util.IntervalTree; //导入依赖的package包/类
void removeBlackListed()
{
if(blackListedBedFile==null) return;
final BedLineCodec codec= new BedLineCodec();
BufferedReader r=null;
int countBefore = this.items.size();
try
	{
	r = IOUtils.openFileForBufferedReading(blackListedBedFile);
	final IntervalTree<Boolean> badRgns = new IntervalTree<>();
	r.lines().map(L->codec.decode(L)).
		filter(B->B!=null).
		filter(B->B.getContig().equals(getContig())).
		forEach(B->{
			badRgns.put(B.getStart(),B.getEnd(), true);
			});
	r.close();r=null;
	
	int i=0;
	while(i<this.items.size())
		{
		GCAndDepth item = this.items.get(i);
		boolean ok=true;
		Iterator<IntervalTree.Node<Boolean>> it = badRgns.overlappers(item.start,item.end);
		while(it.hasNext())
			{
			final IntervalTree.Node<Boolean> node=it.next();
			if(node.getStart() < item.start && item.start< node.getEnd() && node.getEnd()<=item.end)
				{
				item.start=node.getEnd()+1;
				}
			else if(item.start <=node.getStart()  && node.getStart() <=item.end && item.end<node.getEnd())
				{
				item.end=node.getStart()-1;
				}
			else
				{
				ok=false;
				}
			if(item.end-item.start< 0.8*windowSize)
				{
				ok=false;
				}
			}
		if(!ok)
			{
			this.items.remove(i);
			}
		else
			{
			i++;
			}
		}
	if(countBefore!=this.items.size())
		{
		LOG.info("removed "+(countBefore-this.items.size())+" blacklisted regions.");
		}
	}
catch(IOException err)
	{
	throw new RuntimeIOException(err);
	}
finally
	{
	CloserUtil.close(r);
	}	
}