zzh
/
gatk-3.8
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Cleaned up documentation on SecondaryBaseTransitionTableWalker and added Read Group and Allele Balance to the info. 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@2827 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
mmelgar 2010-02-11 17:20:35 +00:00
parent 34e566c90d
commit ad608d0e9d
1 changed files with 29 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

41
java/src/org/broadinstitute/sting/playground/gatk/walkers/secondaryBases/SecondaryBaseTransitionTableWalker.java
View File

@ -3,23 +3,23 @@ package org.broadinstitute.sting.playground.gatk.walkers.secondaryBases;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.walkers.*;
import org.broadinstitute.sting.gatk.walkers.LocusWalker;
import org.broadinstitute.sting.gatk.walkers.Reference;
import org.broadinstitute.sting.gatk.walkers.Window;
import org.broadinstitute.sting.gatk.walkers.genotyper.*;
import org.broadinstitute.sting.playground.utils.NamedTable;
import org.broadinstitute.sting.utils.BaseUtils;
import org.broadinstitute.sting.utils.genotype.Genotype;
import org.broadinstitute.sting.utils.pileup.PileupElement;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.genotype.Genotype;
import java.util.HashMap;
/**
* Created By User: Michael Melgar
* Date: Nov 2, 2009
* Time: 8:45:32 PM
* Given a secondary base annotated .bam file and a reference, this walker generates a table of secondary base counts
* for all called loci in the .bam. Each base call made is an instance included in the table. Specifically, the walker
* maps the following vector to a count of secondary bases:
* <Called Genotype, Reference Base, Primary Base, Previous Genome Base, Read Group, Secondary Base>.
* <HomRef/Het/HomVar, Read Group, Called Genotype, AlleleBalance, Reference Base, Primary Base, Previous Read Base, Secondary Base>.
*/
@Reference(window=@Window(start=-1,stop=1))
@ -41,6 +41,9 @@ public class SecondaryBaseTransitionTableWalker extends LocusWalker<Integer, Int
public Integer map(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
char refBase = Character.toUpperCase(ref.getBase());
ReadBackedPileup pileup = context.getBasePileup();
int[] baseCounts = pileup.getBaseCounts();
int length = 0;
for (int i : baseCounts) {length += i;}
char[] contextBases = ref.getBases();
char prevBase = Character.toUpperCase(contextBases[0]);
char nextBase = Character.toUpperCase(contextBases[contextBases.length - 1]);
@ -51,11 +54,22 @@ public class SecondaryBaseTransitionTableWalker extends LocusWalker<Integer, Int
Genotype res = ugResult.genotypes.get(0);
String call = res.getBases();
String type;
if (!res.isVariant(refBase)) {type = "homref";}
String alleleBalance = "N/A";
if (!res.isVariant(refBase)) {
type = "homref";
}
else if (!res.isHet()) {type = "homvar";}
else if (call.contains(Character.toString(refBase))) {type = "het";}
else if (call.contains(Character.toString(refBase))) {
type = "het";
char alt;
if (call.charAt(0) == refBase) {alt = call.charAt(1);}
else {alt = call.charAt(0);}
double refCount = baseCounts[BaseUtils.simpleBaseToBaseIndex(refBase)];
double altCount = baseCounts[BaseUtils.simpleBaseToBaseIndex(alt)];
alleleBalance = Double.toString(Math.round(100.0*refCount/(refCount + altCount))/100.0);
}
else {type = "bad";}
if (type != "bad") {
if (!type.equals("bad")) {
for (PileupElement element : pileup) {
char primaryBase = Character.toUpperCase((char)element.getBase());
char secondaryBase = Character.toUpperCase((char)element.getSecondBase());
@ -65,20 +79,23 @@ public class SecondaryBaseTransitionTableWalker extends LocusWalker<Integer, Int
String strandPrimary;
String strandPrev;
String strandSecondary;
String strandCall;
if (!element.getRead().getReadNegativeStrandFlag()) {
strandRef = Character.toString(refBase);
strandPrimary = Character.toString(primaryBase);
strandPrev = Character.toString(prevBase);
strandSecondary = Character.toString(secondaryBase);
strandCall = call;
}
else {
strandRef = Character.toString(BaseUtils.simpleComplement(refBase));
strandPrimary = Character.toString(BaseUtils.simpleComplement(primaryBase));
strandPrev = Character.toString(BaseUtils.simpleComplement(nextBase));
strandSecondary = Character.toString(BaseUtils.simpleComplement(secondaryBase));
strandCall = BaseUtils.simpleReverseComplement(call);
}
if (strandPrev.charAt(0) != 'N') {
String key = RG+' '+type+' '+call+' '+strandRef+' '+strandPrimary+' '+strandPrev+' '+strandSecondary;
String key = type+' '+RG+' '+strandCall+' '+alleleBalance+' '+strandRef+' '+strandPrimary+' '+strandPrev+' '+strandSecondary;
if (counts.containsKey(key)) {
counts.put(key, counts.get(key) + Long.valueOf(1));
}
@ -100,10 +117,10 @@ public class SecondaryBaseTransitionTableWalker extends LocusWalker<Integer, Int
public void onTraversalDone(Integer result) {
out.println(">>>");
out.println("ReadGroup CallType CalledGenotype ReferenceBase PrimaryBase PreviousBase SecondaryBase Count");
out.println("Type ReadGroup CalledGenotype AlleleBalance ReferenceBase PrimaryBase PreviousBase SecondaryBase Count");
for (String key : counts.keySet()) {
out.println(key + ' ' + counts.get(key).toString());
}
out.println("Processed " + result.toString() + " loci.");
}
}
}