use of org.broadinstitute.hellbender.utils.haplotype.Haplotype in project gatk by broadinstitute.
the class VectorLoglessPairHMM method computeLog10Likelihoods.
/**
* {@inheritDoc}
*/
@Override
public void computeLog10Likelihoods(final LikelihoodMatrix<Haplotype> logLikelihoods, final List<GATKRead> processedReads, final Map<GATKRead, byte[]> gcp) {
if (processedReads.isEmpty()) {
return;
}
if (doProfiling) {
startTime = System.nanoTime();
}
int readListSize = processedReads.size();
int numHaplotypes = logLikelihoods.numberOfAlleles();
ReadDataHolder[] readDataArray = new ReadDataHolder[readListSize];
int idx = 0;
for (GATKRead read : processedReads) {
readDataArray[idx] = new ReadDataHolder();
readDataArray[idx].readBases = read.getBases();
readDataArray[idx].readQuals = read.getBaseQualities();
readDataArray[idx].insertionGOP = ReadUtils.getBaseInsertionQualities(read);
readDataArray[idx].deletionGOP = ReadUtils.getBaseDeletionQualities(read);
readDataArray[idx].overallGCP = gcp.get(read);
++idx;
}
//to store results
mLogLikelihoodArray = new double[readListSize * numHaplotypes];
if (doProfiling) {
threadLocalSetupTimeDiff = (System.nanoTime() - startTime);
}
//for(reads)
// for(haplotypes)
// compute_full_prob()
pairHmm.computeLikelihoods(readDataArray, mHaplotypeDataArray, mLogLikelihoodArray);
int readIdx = 0;
for (int r = 0; r < readListSize; r++) {
int hapIdx = 0;
for (final Haplotype haplotype : logLikelihoods.alleles()) {
//Since the order of haplotypes in the List<Haplotype> and alleleHaplotypeMap is different,
//get idx of current haplotype in the list and use this idx to get the right likelihoodValue
final int idxInsideHaplotypeList = haplotypeToHaplotypeListIdxMap.get(haplotype);
logLikelihoods.set(hapIdx, r, mLogLikelihoodArray[readIdx + idxInsideHaplotypeList]);
++hapIdx;
}
readIdx += numHaplotypes;
}
if (doProfiling) {
threadLocalPairHMMComputeTimeDiff = (System.nanoTime() - startTime);
pairHMMComputeTime += threadLocalPairHMMComputeTimeDiff;
pairHMMSetupTime += threadLocalSetupTimeDiff;
}
}
use of org.broadinstitute.hellbender.utils.haplotype.Haplotype in project gatk by broadinstitute.
the class ReadThreadingAssemblerUnitTest method testAssemblyWithVariant.
private void testAssemblyWithVariant(final ReadThreadingAssembler assembler, final byte[] refBases, final SimpleInterval loc, final int nReadsToUse, final VariantContext site) {
final String preRef = new String(refBases).substring(0, site.getStart());
final String postRef = new String(refBases).substring(site.getEnd() + 1, refBases.length);
final byte[] altBases = (preRef + site.getAlternateAllele(0).getBaseString() + postRef).getBytes();
// logger.warn("ref " + new String(refBases));
// logger.warn("alt " + new String(altBases));
final List<GATKRead> reads = new LinkedList<>();
for (int i = 0; i < nReadsToUse; i++) {
final byte[] bases = altBases.clone();
final byte[] quals = Utils.dupBytes((byte) 30, altBases.length);
final String cigar = altBases.length + "M";
final GATKRead read = ArtificialReadUtils.createArtificialRead(header, loc.getContig(), loc.getContig(), loc.getStart(), bases, quals, cigar);
reads.add(read);
}
final Haplotype refHaplotype = new Haplotype(refBases, true);
final Haplotype altHaplotype = new Haplotype(altBases, false);
final List<Haplotype> haplotypes = assemble(assembler, refBases, loc, reads);
Assert.assertEquals(haplotypes, Arrays.asList(refHaplotype, altHaplotype));
}
use of org.broadinstitute.hellbender.utils.haplotype.Haplotype in project gatk by broadinstitute.
the class ReadThreadingAssemblerUnitTest method testAssembleRef.
@Test(dataProvider = "AssembleIntervalsData")
public void testAssembleRef(final ReadThreadingAssembler assembler, final SimpleInterval loc, final int nReadsToUse) {
final byte[] refBases = seq.getSubsequenceAt(loc.getContig(), loc.getStart(), loc.getEnd()).getBases();
final List<GATKRead> reads = new LinkedList<>();
for (int i = 0; i < nReadsToUse; i++) {
final byte[] bases = refBases.clone();
final byte[] quals = Utils.dupBytes((byte) 30, refBases.length);
final String cigar = refBases.length + "M";
final GATKRead read = ArtificialReadUtils.createArtificialRead(header, loc.getContig(), loc.getContig(), loc.getStart(), bases, quals, cigar);
reads.add(read);
}
// TODO -- generalize to all assemblers
final Haplotype refHaplotype = new Haplotype(refBases, true);
final List<Haplotype> haplotypes = assemble(assembler, refBases, loc, reads);
Assert.assertEquals(haplotypes, Collections.singletonList(refHaplotype));
}
use of org.broadinstitute.hellbender.utils.haplotype.Haplotype in project gatk by broadinstitute.
the class ReadThreadingAssemblerUnitTest method testSimpleAssembly.
@Test(dataProvider = "SimpleAssemblyTestData")
public void testSimpleAssembly(final String name, final ReadThreadingAssembler assembler, final SimpleInterval loc, final String ref, final String alt) {
final byte[] refBases = ref.getBytes();
final byte[] altBases = alt.getBytes();
final List<GATKRead> reads = new LinkedList<>();
for (int i = 0; i < 20; i++) {
final byte[] bases = altBases.clone();
final byte[] quals = Utils.dupBytes((byte) 30, altBases.length);
final String cigar = altBases.length + "M";
final GATKRead read = ArtificialReadUtils.createArtificialRead(header, loc.getContig(), loc.getContig(), loc.getStart(), bases, quals, cigar);
reads.add(read);
}
final Haplotype refHaplotype = new Haplotype(refBases, true);
final Haplotype altHaplotype = new Haplotype(altBases, false);
final List<Haplotype> haplotypes = assemble(assembler, refBases, loc, reads);
Assert.assertTrue(haplotypes.size() > 0, "Failed to find ref haplotype");
Assert.assertEquals(haplotypes.get(0), refHaplotype);
Assert.assertEquals(haplotypes.size(), 2, "Failed to find single alt haplotype");
Assert.assertEquals(haplotypes.get(1), altHaplotype);
}
use of org.broadinstitute.hellbender.utils.haplotype.Haplotype in project gatk by broadinstitute.
the class AssemblyRegionTestDataSet method cigarToHaplotype.
private Haplotype cigarToHaplotype(final String reference, final String cigar, final int offset, final boolean global) {
final String sequence = applyCigar(reference, cigar, offset, global);
final Haplotype haplotype = new Haplotype(sequence.getBytes(), reference.equals(sequence));
haplotype.setGenomeLocation(genomeLocParser.createGenomeLoc("chr1", 1, reference.length()));
haplotype.setCigar(Civar.fromCharSequence(cigar).toCigar(reference.length()));
return haplotype;
}
Aggregations