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

Added some methods to retreive the probability distributions of individual bases. 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@484 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
kiran 2009-04-21 22:26:25 +00:00
parent 5b2a7c9c23
commit dac76f041b
2 changed files with 55 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

50
java/src/org/broadinstitute/sting/utils/BasicPileup.java
View File

@ -140,6 +140,56 @@ abstract public class BasicPileup implements Pileup {
return quals2.toString();
}
public static double[][] probDistPileup( List<SAMRecord> reads, List<Integer> offsets ) {
double[][] dist = new double[reads.size()][4];
for (int readIndex = 0; readIndex < dist.length; readIndex++) {
SAMRecord read = reads.get(readIndex);
String bases = read.getReadString();
int offset = offsets.get(readIndex);
int bestBaseIndex = BaseUtils.simpleBaseToBaseIndex(bases.charAt(offset));
if (bestBaseIndex >= 0 && bestBaseIndex < 4) {
dist[readIndex][bestBaseIndex] = QualityUtils.qualToProb(read.getBaseQualities()[offset]);
byte[] sqs = (byte[]) read.getAttribute("SQ");
if (sqs != null) {
int secondBestBaseIndex = QualityUtils.compressedQualityToBaseIndex(sqs[offset]);
dist[readIndex][secondBestBaseIndex] = (1.0 - dist[readIndex][bestBaseIndex]);
} else {
for (int baseIndex = 0; baseIndex < 4; baseIndex++) {
if (baseIndex != bestBaseIndex) {
dist[readIndex][baseIndex] = (dist[readIndex][bestBaseIndex]/3.0);
}
}
}
} else {
for (int baseIndex = 0; baseIndex < 4; baseIndex++) {
dist[readIndex][baseIndex] = 0.25;
}
}
}
return dist;
}
public static String probDistPileupAsString( List<SAMRecord> reads, List<Integer> offsets ) {
double[][] dist = probDistPileup(reads, offsets);
String distString = "";
for (int readIndex = 0; readIndex < dist.length; readIndex++) {
distString += "[ ";
for (int baseIndex = 0; baseIndex < 4; baseIndex++) {
distString += String.format("%3.3f ", dist[readIndex][baseIndex]);
}
distString += "]\n";
}
return distString;
}
public static String pileupDiff(final Pileup a, final Pileup b)
{
return pileupDiff(a,b,true);

6
java/src/org/broadinstitute/sting/utils/ReadBackedPileup.java
View File

@ -69,7 +69,11 @@ public class ReadBackedPileup extends BasicPileup {
public String getSecondaryQualPileup() {
return secondaryQualPileupAsString(reads, offsets);
}
public String getProbDistPileup() {
return probDistPileupAsString(reads, offsets);
}
public String getPileupString()
{
return String.format("%s: %s %s %s", getLocation(), getRef(), getBases(), getQuals());