Java Code Examples for htsjdk.samtools.util.StringUtil#toUpperCase()

public static int[] calculateRefWindowsByGc(final int windows, final File referenceSequence, final int windowSize) {
    final ReferenceSequenceFile refFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(referenceSequence);
    ReferenceSequence ref;

    final int [] windowsByGc = new int [windows];

    while ((ref = refFile.nextSequence()) != null) {
        final byte[] refBases = ref.getBases();
        StringUtil.toUpperCase(refBases);
        final int refLength = refBases.length;
        final int lastWindowStart = refLength - windowSize;

        final CalculateGcState state = new GcBiasUtils().new CalculateGcState();

        for (int i = 1; i < lastWindowStart; ++i) {
            final int windowEnd = i + windowSize;
            final int gcBin = calculateGc(refBases, i, windowEnd, state);
            if (gcBin != -1) windowsByGc[gcBin]++;
        }
    }

    return windowsByGc;
}
 

/**
 * Create one SAMSequenceRecord from a single fasta sequence
 */
private SAMSequenceRecord makeSequenceRecord(final ReferenceSequence refSeq) {
    final SAMSequenceRecord ret = new SAMSequenceRecord(refSeq.getName(), refSeq.length());

    // Compute MD5 of upcased bases
    final byte[] bases = refSeq.getBases();
    for (int i = 0; i < bases.length; ++i) {
        bases[i] = StringUtil.toUpperCase(bases[i]);
    }

    ret.setAttribute(SAMSequenceRecord.MD5_TAG, md5Hash(bases));
    if (genomeAssembly != null) {
        ret.setAttribute(SAMSequenceRecord.ASSEMBLY_TAG, genomeAssembly);
    }
    ret.setAttribute(SAMSequenceRecord.URI_TAG, uri);
    if (species != null) {
        ret.setAttribute(SAMSequenceRecord.SPECIES_TAG, species);
    }
    return ret;
}
 

/**
 * Create one SAMSequenceRecord from a single fasta sequence
 */
private SAMSequenceRecord makeSequenceRecord(final ReferenceSequence refSeq) {
  final SAMSequenceRecord ret = new SAMSequenceRecord(refSeq.getName(), refSeq.length());

  // Compute MD5 of upcased bases
  final byte[] bases = refSeq.getBases();
  for (int i = 0; i < bases.length; ++i) {
    bases[i] = StringUtil.toUpperCase(bases[i]);
  }

  ret.setAttribute(SAMSequenceRecord.MD5_TAG, md5Hash(bases));
  if (GENOME_ASSEMBLY != null) {
    ret.setAttribute(SAMSequenceRecord.ASSEMBLY_TAG, GENOME_ASSEMBLY);
  }
  ret.setAttribute(SAMSequenceRecord.URI_TAG, URI);
  if (SPECIES != null) {
    ret.setAttribute(SAMSequenceRecord.SPECIES_TAG, SPECIES);
  }
  return ret;
}
 

public Snp(final String name, final String chrom, final int pos, final byte allele1, final byte allele2,
           final double maf, final List<String> fingerprintPanels) {
    this.name = name;
    this.chrom = chrom;
    this.pos = pos;
    this.allele1 = StringUtil.toUpperCase(allele1);
    this.allele2 = StringUtil.toUpperCase(allele2);
    this.maf = maf;
    this.fingerprintPanels = fingerprintPanels == null ? new ArrayList<String>() : fingerprintPanels;

    // Construct the genotypes for ease of comparison
    this.genotypes[0] = DiploidGenotype.fromBases(allele1, allele1);
    this.genotypes[1] = DiploidGenotype.fromBases(allele1, allele2);
    this.genotypes[2] = DiploidGenotype.fromBases(allele2, allele2);
}
 

@Override
public void acceptRecord(final GcBiasCollectorArgs args) {
    final SAMRecord rec = args.getRec();
    if (logCounter < 100 && rec.getReadBases().length == 0) {
        log.warn("Omitting read " + rec.getReadName() + " with '*' in SEQ field.");
        if (++logCounter == 100) {
            log.warn("There are more than 100 reads with '*' in SEQ field in file.");
        }
        return;
    }
    if (!rec.getReadUnmappedFlag()) {
        if (referenceIndex != rec.getReferenceIndex() || gc == null) {
            final ReferenceSequence ref = args.getRef();
            refBases = ref.getBases();
            StringUtil.toUpperCase(refBases);
            final int refLength = refBases.length;
            final int lastWindowStart = refLength - scanWindowSize;
            gc = GcBiasUtils.calculateAllGcs(refBases, lastWindowStart, scanWindowSize);
            referenceIndex = rec.getReferenceIndex();
        }

        addReadToGcData(rec, this.gcData);
        if (ignoreDuplicates && !rec.getDuplicateReadFlag()) {
            addReadToGcData(rec, this.gcDataNonDups);
        }
    } else {
        updateTotalClusters(rec, this.gcData);
        if (ignoreDuplicates && !rec.getDuplicateReadFlag()) {
            updateTotalClusters(rec, this.gcDataNonDups);
        }
    }
}
 

private static VariantContext getVariantContext(final ReferenceSequenceFile reference,
                                                final String sample,
                                                final HaplotypeProbabilities haplotypeProbabilities) {
    final Snp snp = haplotypeProbabilities.getRepresentativeSnp();
    final byte refAllele = StringUtil.toUpperCase(reference.getSubsequenceAt(
            snp.getChrom(),
            snp.getPos(),
            snp.getPos()).getBases()[0]);

    if (snp.getAllele1() != refAllele && snp.getAllele2() != refAllele) {
        throw new PicardException("Don't know how to deal with missing reference allele in fingerprinting map");
    }

    final Allele alleleRef;
    final Allele alleleAlt;
    final int obsRef, obsAlt;
    final boolean swap12 = snp.getAllele2() == refAllele;

    if (swap12) {
        alleleRef = Allele.create(snp.getAllele2(), true);
        alleleAlt = Allele.create(snp.getAllele1(), false);
        obsRef = haplotypeProbabilities.getObsAllele2();
        obsAlt = haplotypeProbabilities.getObsAllele1();
    } else {
        alleleRef = Allele.create(snp.getAllele1(), true);
        alleleAlt = Allele.create(snp.getAllele2(), false);
        obsRef = haplotypeProbabilities.getObsAllele1();
        obsAlt = haplotypeProbabilities.getObsAllele2();
    }

    final double[] origPLs = haplotypeProbabilities.getLogLikelihoods();
    final double[] PLs =  Arrays.copyOf(origPLs,origPLs.length);
    if (swap12) {
        ArrayUtils.reverse(PLs);
    }
    final List<Allele> alleles = Arrays.asList(alleleRef, alleleAlt);

    final Genotype gt = new GenotypeBuilder()
            .DP(haplotypeProbabilities.getTotalObs())
            .noAttributes()
            .PL(PLs)
            .AD(new int[]{obsAlt, obsRef})
            .name(sample)
            .make();
    try {
        return new VariantContextBuilder(
                snp.getName(),
                snp.getChrom(),
                snp.getPos(),
                snp.getPos(),
                alleles)
                .log10PError(VariantContext.NO_LOG10_PERROR)
                .genotypes(gt)
                .unfiltered().make();
    } catch (IllegalArgumentException e) {
        throw new IllegalArgumentException(String.format("Trouble creating variant at %s-%d", snp.getChrom(), snp.getPos()), e);
    }
}
 

@Override
protected int doWork() {
    IOUtil.assertFileIsReadable(INPUT);
    IOUtil.assertFileIsWritable(OUTPUT);

    // set up the reference and a mask so that we only count the positions requested by the user
    final ReferenceSequenceFile ref = ReferenceSequenceFileFactory.getReferenceSequenceFile(INPUT);
    final ReferenceSequenceMask referenceSequenceMask;
    if (INTERVALS != null) {
        IOUtil.assertFileIsReadable(INTERVALS);
        final IntervalList intervalList = IntervalList.fromFile(INTERVALS);
        referenceSequenceMask = new IntervalListReferenceSequenceMask(intervalList);
    } else {
        final SAMFileHeader header = new SAMFileHeader();
        header.setSequenceDictionary(ref.getSequenceDictionary());
        referenceSequenceMask = new WholeGenomeReferenceSequenceMask(header);
    }

    long nonNbases = 0L;

    for (final SAMSequenceRecord rec : ref.getSequenceDictionary().getSequences()) {
        // pull out the contig and set up the bases
        final ReferenceSequence sequence = ref.getSequence(rec.getSequenceName());
        final byte[] bases = sequence.getBases();
        StringUtil.toUpperCase(bases);

        for (int i = 0; i < bases.length; i++) {
            // only investigate this position if it's within our mask
            if (referenceSequenceMask.get(sequence.getContigIndex(), i+1)) {
                nonNbases += bases[i] == SequenceUtil.N ? 0 : 1;
            }
        }
    }

    try {
        final BufferedWriter out = IOUtil.openFileForBufferedWriting(OUTPUT);
        out.write(nonNbases + "\n");
        out.close();
    }
    catch (IOException ioe) {
        throw new PicardException("Error writing to file " + OUTPUT.getAbsolutePath(), ioe);
    }

    return 0;
}