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

better printing in VE2. Added support for TiTv analysis 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@2766 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
depristo 2010-02-02 21:20:29 +00:00
parent cbbc0e98d2
commit fa2cd432fd
2 changed files with 116 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

65
java/src/org/broadinstitute/sting/oneoffprojects/variantcontext/varianteval2/TiTvVariantEvaluator.java 100755
View File

@ -0,0 +1,65 @@
package org.broadinstitute.sting.oneoffprojects.variantcontext.varianteval2;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.oneoffprojects.variantcontext.VariantContext;
import org.broadinstitute.sting.oneoffprojects.variantcontext.Genotype;
import org.broadinstitute.sting.utils.StingException;
import java.util.List;
import java.util.Arrays;
public class TiTvVariantEvaluator extends VariantEvaluator {
long nTi = 0, nTv = 0;
long nTiInStd = 0, nTvInStd = 0;
private double ratio(long num, long denom) {
return ((double)num) / (Math.max(denom, 1));
}
public String getName() {
return "titv";
}
public int getComparisonOrder() {
return 2; // we only need to see each eval track
}
public void updateTiTv(VariantContext vc, boolean updateStandard) {
if ( vc != null && vc.isSNP() && vc.isBiallelic() ) {
if ( vc.isTransition() ) {
if ( updateStandard ) nTiInStd++; else nTi++;
} else {
if ( updateStandard ) nTvInStd++; else nTv++;
}
}
}
public void update2(VariantContext vc1, VariantContext vc2, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
if ( vc1 != null ) updateTiTv(vc1, false);
if ( vc2 != null ) updateTiTv(vc2, true);
}
public String toString() {
return getName() + ": " + summaryLine();
}
private String summaryLine() {
return String.format("%d %d %.2f %d %d %.2f",
nTi, nTv, ratio(nTi, nTv),
nTiInStd, nTvInStd, ratio(nTiInStd, nTvInStd));
}
private static List<String> HEADER =
Arrays.asList("nTi", "nTv", "TiTvRatio", "nTiStandard", "nTvStandard", "TiTvRatioStandard");
// making it a table
public List<String> getTableHeader() {
return HEADER;
}
public List<List<String>> getTableRows() {
return Arrays.asList(Arrays.asList(summaryLine().split(" ")));
}
}

56
java/src/org/broadinstitute/sting/oneoffprojects/variantcontext/varianteval2/VariantEval2Walker.java
View File

@ -69,13 +69,15 @@ public class VariantEval2Walker extends RodWalker<Integer, Integer> {
for ( VariantContextUtils.MatchExp e : selectExps ) {
addNewContext(d.getName(), d.getName() + "." + e.name, e);
}
addNewContext(d.getName(), d.getName(), null);
addNewContext(d.getName(), d.getName() + ".all", null);
} else if ( d.getName().startsWith("dbsnp") || d.getName().startsWith("hapmap") || d.getName().startsWith("comp") ) {
compNames.add(d.getName());
} else {
logger.info("Not evaluating ROD binding " + d.getName());
}
}
determineContextNamePartSizes();
}
private void determineAllEvalations() {
@ -144,7 +146,7 @@ public class VariantEval2Walker extends RodWalker<Integer, Integer> {
// now call the single or paired update function
switch ( evaluation.getComparisonOrder() ) {
case 1:
if ( vc != null ) evaluation.update1(vc, tracker, ref, context);
if ( evalWantsVC && vc != null ) evaluation.update1(vc, tracker, ref, context);
break;
case 2:
for ( VariantContext comp : comps.values() ) {
@ -249,16 +251,59 @@ public class VariantEval2Walker extends RodWalker<Integer, Integer> {
return null;
}
public <T extends Comparable<T>> List<T> sorted(Collection<T> c ) {
private <T extends Comparable<T>> List<T> sorted(Collection<T> c ) {
List<T> l = new ArrayList<T>(c);
Collections.sort(l);
return l;
}
private final static String CONTEXT_HEADER = "track.subset.novelty.filter";
private final static int N_CONTEXT_NAME_PARTS = CONTEXT_HEADER.split("\\.").length;
private static int[] nameSizes = new int[N_CONTEXT_NAME_PARTS];
static {
int i = 0;
for ( String elt : CONTEXT_HEADER.split("\\.") )
nameSizes[i++] = elt.length();
}
private void determineContextNamePartSizes() {
for ( String contextName : sorted(contexts.keySet()) ) {
EvaluationContext group = contexts.get(contextName);
for ( String evalSubgroupName : sorted(group.keySet()) ) {
String keyWord = contextName + "." + evalSubgroupName;
String[] parts = keyWord.split("\\.");
if ( parts.length != N_CONTEXT_NAME_PARTS ) {
throw new StingException("Unexpected number of eval name parts " + keyWord + " length = " + parts.length + ", expected " + N_CONTEXT_NAME_PARTS);
} else {
for ( int i = 0; i < parts.length; i++ )
nameSizes[i] = Math.max(nameSizes[i], parts[i].length());
}
}
}
}
private String formatKeyword(String keyWord) {
//System.out.printf("keyword %s%n", keyWord);
StringBuilder s = new StringBuilder();
int i = 0;
for ( String part : keyWord.split("\\.") ) {
//System.out.printf("part %s %d%n", part, nameSizes[i]);
s.append(String.format("%"+nameSizes[i]+"s ", part));
i++;
}
return s.toString();
}
public void onTraversalDone(Integer result) {
// todo -- this really needs to be pretty printed; use some kind of table organization
// todo -- so that we can load up all of the data in one place, analyze the widths of the columns
// todo -- and pretty print it
for ( String evalName : variantEvaluationNames ) {
boolean first = true;
out.printf("%n%n");
// todo -- show that comp is dbsnp, etc. is columns
for ( String contextName : sorted(contexts.keySet()) ) {
EvaluationContext group = contexts.get(contextName);
@ -267,13 +312,14 @@ public class VariantEval2Walker extends RodWalker<Integer, Integer> {
Set<VariantEvaluator> evalSet = group.get(evalSubgroupName);
VariantEvaluator eval = getEvalByName(evalName, evalSet);
String keyWord = contextName + "." + evalSubgroupName;
if ( first ) {
out.printf("%20s\t%40s\t%s%n", evalName, "context", Utils.join("\t\t", eval.getTableHeader()));
out.printf("%20s %s %s%n", evalName, formatKeyword(CONTEXT_HEADER), Utils.join("\t\t", eval.getTableHeader()));
first = false;
}
for ( List<String> row : eval.getTableRows() )
out.printf("%20s\t%40s\t%s%n", evalName, keyWord, Utils.join("\t\t", row));
out.printf("%20s %s %s%n", evalName, formatKeyword(keyWord), Utils.join("\t\t", row));
}
}
}