Java Code Examples for htsjdk.variant.variantcontext.VariantContextBuilder#noID()

private static VariantContext getContextFromVariant(VCFHeader header, Variant in) {
  // Create allele list
  String refBases = in.getReferenceBases();
  List<String> altBases = in.getAlternateBasesList();
  List<Allele> alleleList = new ArrayList<>();
  alleleList.add(Allele.create(refBases.getBytes(StandardCharsets.UTF_8), true));
  for (String base : altBases) {
    alleleList.add(Allele.create(base.getBytes(StandardCharsets.UTF_8), false));
  }

  // Note htsjdk start/end are both 1-based closed
  VariantContextBuilder builder = new VariantContextBuilder(SOURCE, in.getReferenceName(),
      in.getStart() + 1, in.getEnd(), alleleList);

  // Set Quality
  if (in.getQuality() != 0.0) {
    builder.log10PError(-in.getQuality() / 10);
  }

  // Set names
  int numNames = in.getNamesCount();
  int i = 0;
  String namesString = "";
  for (String name : in.getNamesList()) {
    if (i == numNames - 1) {
      break;
    }
    namesString += name + ",";
  }
  if (numNames == 0) {
    builder.noID();
  } else {
    // Also handle fence post
    builder.id(namesString + in.getNames(numNames - 1));
  }

  // Set filters
  boolean hadFilter = false;
  for (String filter : in.getFilterList()) {
    hadFilter = true;
    if (filter.equals(VCFConstants.PASSES_FILTERS_v4)) {
      builder.passFilters();
      break;
    }
    if (filter.equals(VCFConstants.MISSING_VALUE_v4)) {
      builder.unfiltered();
      break;
    }

    builder.filter(filter);
  }
  if (!hadFilter) {
    builder.unfiltered();
  }

  // Set info
  Map<String, Object> toSet = new HashMap<>();
  for (Map.Entry<String, ListValue> entry : in.getInfo().entrySet()) {
    String currKey = entry.getKey();
    List<Value> currValue = entry.getValue().getValuesList();
    int currSize = currValue.size();

    if (currSize == 0) {
      // If the key is present in the info map, there must
      // be non-zero attributes associated with it
      throw new IllegalStateException(String.format("info field %s has length 0", currKey));
    }

    VCFHeaderLineType type = header.getInfoHeaderLine(currKey).getType();
    if (currSize == 1) {
      toSet.put(currKey, parseObject(currKey, currValue.get(0), type));
      continue;
    }

    List<Object> objectList = new ArrayList<>();
    for (Value val : currValue) {
      objectList.add(parseObject(currKey, val, type));
    }
    toSet.put(currKey, objectList);
  }
  builder.attributes(toSet);

  // Set calls
  List<Genotype> genotypes = new ArrayList<>();
  for (VariantCall vc : in.getCallsList()) {
    genotypes.add(getGenotypeFromVariantCall(header, vc, alleleList));
  }
  if (genotypes.size() == 0) {
    builder.noGenotypes();
  } else {
    builder.genotypes(genotypes);
  }

  return builder.make();
}
 

public static VariantContext read(final DataInput in) throws IOException {
	final VariantContextBuilder builder = new VariantContextBuilder();

	int count, len;
	byte[] b;

	len = in.readInt();
	b = new byte[len];
	in.readFully(b);
	final String chrom = new String(b, "UTF-8");
	builder.chr(chrom);

	final int start = in.readInt();
	builder.start(start);
	builder.stop (in.readInt());

	len = in.readInt();
	if (len == 0)
		builder.noID();
	else {
		if (len > b.length) b = new byte[len];
		in.readFully(b, 0, len);
		builder.id(new String(b, 0, len, "UTF-8"));
	}

	count = in.readInt();
	final List<Allele> alleles = new ArrayList<Allele>(count);
	for (int i = 0; i < count; ++i) {
		len = in.readInt();
		if (len > b.length) b = new byte[len];
		in.readFully(b, 0, len);
		alleles.add(Allele.create(Arrays.copyOf(b, len), i == 0));
	}
	builder.alleles(alleles);

	final int qualInt = in.readInt();
	builder.log10PError(
		qualInt == 0x7f800001
			? VariantContext.NO_LOG10_PERROR
			: Float.intBitsToFloat(qualInt));

	count = in.readInt();
	switch (count) {
	case -2: builder.unfiltered(); break;
	case -1: builder.passFilters(); break;
	default:
		while (count-- > 0) {
			len = in.readInt();
			if (len > b.length) b = new byte[len];
			in.readFully(b, 0, len);
			builder.filter(new String(b, 0, len, "UTF-8"));
		}
		break;
	}

	count = in.readInt();
	final Map<String,Object> attrs = new HashMap<String,Object>(count, 1);
	while (count-- > 0) {
		len = in.readInt();
		if (len > b.length) b = new byte[len];
		in.readFully(b, 0, len);
		attrs.put(new String(b, 0, len, "UTF-8"), decodeAttrVal(in));
	}
	builder.attributes(attrs);

	count = in.readInt();
	final byte genoType = in.readByte();
	len = in.readInt();

	// Resize b even if it's already big enough, minimizing the amount of
	// memory LazyGenotypesContext hangs on to.
	b = new byte[len];
	in.readFully(b);

	switch (genoType) {
	case 0:
		builder.genotypesNoValidation(
			new LazyVCFGenotypesContext(alleles, chrom, start, b, count));
		break;

	case 1:
		builder.genotypesNoValidation(
			new LazyBCFGenotypesContext(alleles, in.readInt(), b, count));
		break;

	default:
		throw new IOException(
			"Invalid genotypes type identifier: cannot decode");
	}

	return builder.make();
}
 

@Override
public void apply(AlignmentContext alignmentContext, ReferenceContext referenceContext, FeatureContext featureContext) {
    final byte referenceBase = referenceContext.getBases()[0];
    final int unconvertedBases;
    final int convertedBases;
    final byte alt;
    byte[] context = null;

    // check the forward strand for methylated coverage
    if (referenceBase == (byte)'C') {
        alt = (byte)'T';
        final ReadPileup forwardBasePileup = alignmentContext.stratify(AlignmentContext.ReadOrientation.FORWARD).getBasePileup();
        // unconverted: C, index=1; converted: T, index=3
        final int[] forwardBaseCounts = forwardBasePileup.getBaseCounts();
        unconvertedBases = forwardBaseCounts[BaseUtils.simpleBaseToBaseIndex((byte)'C')];
        convertedBases = forwardBaseCounts[BaseUtils.simpleBaseToBaseIndex((byte)'T')];

        // if there is methylated coverage
        if (unconvertedBases + convertedBases > 0) {
            context = referenceContext.getBases(0, REFERENCE_CONTEXT_LENGTH);
        }
    }
    // check the reverse strand for methylated coverage
    else if (referenceBase == (byte)'G') {
        alt = (byte)'A';
        final ReadPileup reverseBasePileup = alignmentContext.stratify(AlignmentContext.ReadOrientation.REVERSE).getBasePileup();
        // unconverted: G, index=2; converted: A, index=0
        final int[] reverseBaseCounts = reverseBasePileup.getBaseCounts();
        unconvertedBases = reverseBaseCounts[BaseUtils.simpleBaseToBaseIndex((byte)'G')];
        convertedBases = reverseBaseCounts[BaseUtils.simpleBaseToBaseIndex((byte)'A')];

        // if there is methylated coverage
        if (unconvertedBases + convertedBases > 0) {
            // get the reverse complement for context b/c we are on the reverse strand
            context = BaseUtils.simpleReverseComplement(referenceContext.getBases(REFERENCE_CONTEXT_LENGTH,0));
        }
    }
    // if reference strand does not support methylation
    else {
        return;
    }

    // if there are reads that have methylated coverage
    if (context != null) {
        final LinkedHashSet<Allele> alleles = new LinkedHashSet<>();
        alleles.add(Allele.create(referenceBase, true));
        alleles.add(Allele.create(alt, false));

        final VariantContextBuilder vcb = new VariantContextBuilder();
        vcb.chr(alignmentContext.getContig());
        vcb.start(alignmentContext.getPosition());
        vcb.stop(alignmentContext.getPosition());
        vcb.noID();
        vcb.alleles(alleles);
        vcb.noGenotypes();
        vcb.unfiltered();
        vcb.attribute(GATKVCFConstants.UNCONVERTED_BASE_COVERAGE_KEY, unconvertedBases);
        vcb.attribute(GATKVCFConstants.CONVERTED_BASE_COVERAGE_KEY, convertedBases);
        vcb.attribute(GATKVCFConstants.METHYLATION_REFERENCE_CONTEXT_KEY, new String(context));
        vcb.attribute(VCFConstants.DEPTH_KEY, alignmentContext.size());

        // write to VCF
        final VariantContext vc = vcb.make();
        vcfWriter.add(vc);
    }
}