use of uk.ac.babraham.SeqMonk.DataTypes.Genome.Feature in project SeqMonk by s-andrews.
the class BisulphiteFeaturePipeline method startPipeline.
protected void startPipeline() {
// We first need to generate probes over all of the features listed in
// the feature types. The probes should cover the whole area of the
// feature regardless of where it splices.
Vector<Probe> probes = new Vector<Probe>();
int minCount = optionsPanel.getMinCount();
int minFeatureCountValue = optionsPanel.getMinFeatureCount();
String combineString = optionsPanel.getCombineType();
int combineType = 0;
if (combineString.equals("Mean")) {
combineType = MEAN;
} else if (combineString.equals("Median")) {
combineType = MEDIAN;
} else if (combineString.equals("Max")) {
combineType = MAX;
} else if (combineString.equals("Min")) {
combineType = MIN;
} else {
throw new IllegalStateException("Unknown combine type '" + combineString + "' found in methylation pipeline options");
}
Chromosome[] chrs = collection().genome().getAllChromosomes();
if (optionsPanel.getSelectedFeatureType() != "[Existing Probes]") {
for (int c = 0; c < chrs.length; c++) {
if (cancel) {
progressCancelled();
return;
}
progressUpdated("Making features for chr" + chrs[c].name(), c, chrs.length * 2);
Feature[] features = collection().genome().annotationCollection().getFeaturesForType(chrs[c], optionsPanel.getSelectedFeatureType());
Arrays.sort(features);
for (int f = 0; f < features.length; f++) {
if (cancel) {
progressCancelled();
return;
}
probes.add(new Probe(chrs[c], features[f].location().start(), features[f].location().end(), features[f].location().strand(), features[f].name()));
}
}
Probe[] allProbes = probes.toArray(new Probe[0]);
collection().setProbeSet(new ProbeSet("Features over " + optionsPanel.getSelectedFeatureType(), allProbes));
}
for (int c = 0; c < chrs.length; c++) {
if (cancel) {
progressCancelled();
return;
}
progressUpdated("Quantitating features on chr" + chrs[c].name(), chrs.length + c, chrs.length * 2);
Probe[] thisChrProbes = collection().probeSet().getProbesForChromosome(chrs[c]);
for (int p = 0; p < thisChrProbes.length; p++) {
// This data structure is going to hold the counts. The int array will hold
// pairs of values for the different data stores (meth and unmeth) so the length
// of this array will be data store length * 2.
Hashtable<Integer, int[]> counts = new Hashtable<Integer, int[]>();
for (int d = 0; d < data.length; d++) {
if (cancel) {
progressCancelled();
return;
}
long[] reads = data[d].getReadsForProbe(thisChrProbes[p]);
for (int r = 0; r < reads.length; r++) {
for (int pos = SequenceRead.start(reads[r]); pos <= SequenceRead.end(reads[r]); pos++) {
if (pos < thisChrProbes[p].start())
continue;
if (pos > thisChrProbes[p].end())
break;
if (!counts.containsKey(pos)) {
counts.put(pos, new int[data.length * 2]);
}
if (SequenceRead.strand(reads[r]) == Location.FORWARD) {
counts.get(pos)[2 * d]++;
} else if (SequenceRead.strand(reads[r]) == Location.REVERSE) {
counts.get(pos)[(2 * d) + 1]++;
}
}
}
}
if (optionsPanel.applyMinCountAcrossAllStores()) {
Enumeration<Integer> en = counts.keys();
OUTER: while (en.hasMoreElements()) {
Integer i = en.nextElement();
int[] values = counts.get(i);
for (int d = 0; d < data.length; d++) {
if (values[2 * d] + values[(2 * d) + 1] < minCount) {
counts.remove(i);
continue OUTER;
}
}
}
}
for (int d = 0; d < data.length; d++) {
Enumeration<Integer> en = counts.keys();
Vector<Float> percentages = new Vector<Float>();
while (en.hasMoreElements()) {
Integer pos = en.nextElement();
int total = counts.get(pos)[2 * d] + counts.get(pos)[(2 * d) + 1];
if (total < minCount)
continue;
float percent = 100f * (counts.get(pos)[2 * d] / (float) total);
percentages.add(percent);
}
// Now combine the methylation percentages we have in the way they
// wanted.
float[] validMeasures = new float[percentages.size()];
float finalMeasure = Float.NaN;
for (int i = 0; i < validMeasures.length; i++) {
validMeasures[i] = percentages.elementAt(i).floatValue();
}
if (validMeasures.length < minFeatureCountValue) {
finalMeasure = Float.NaN;
} else if (validMeasures.length == 0) {
finalMeasure = Float.NaN;
} else if (validMeasures.length == 1) {
finalMeasure = validMeasures[0];
} else {
Arrays.sort(validMeasures);
switch(combineType) {
case MEAN:
finalMeasure = SimpleStats.mean(validMeasures);
break;
case MEDIAN:
finalMeasure = SimpleStats.median(validMeasures);
break;
case MIN:
finalMeasure = validMeasures[0];
break;
case MAX:
finalMeasure = validMeasures[validMeasures.length - 1];
break;
}
}
data[d].setValueForProbe(thisChrProbes[p], finalMeasure);
}
}
}
if (optionsPanel.getSelectedFeatureType() != "[Existing Probes]") {
collection().probeSet().setDescription("Probes over " + optionsPanel.getSelectedFeatureType() + " features");
}
StringBuffer quantitationDescription = new StringBuffer();
quantitationDescription.append("Methylation feature pipeline quantitation using ");
quantitationDescription.append(optionsPanel.getSelectedFeatureType());
quantitationDescription.append(" features with min count ");
if (optionsPanel.applyMinCountAcrossAllStores()) {
quantitationDescription.append(" in all stores ");
}
quantitationDescription.append("of ");
quantitationDescription.append(optionsPanel.getMinCount());
quantitationDescription.append(" per position, and at least ");
quantitationDescription.append(optionsPanel.getMinFeatureCount());
quantitationDescription.append(" observations per feature, then combining using the ");
quantitationDescription.append(optionsPanel.getCombineType());
collection().probeSet().setCurrentQuantitation(quantitationDescription.toString());
quantitatonComplete();
}
use of uk.ac.babraham.SeqMonk.DataTypes.Genome.Feature in project SeqMonk by s-andrews.
the class GeneTrapQuantitationPipeline method startPipeline.
protected void startPipeline() {
// We first need to generate probes over all of the features listed in
// the feature types. The probes should cover the whole area of the
// feature regardless of where it splices.
Vector<Probe> probes = new Vector<Probe>();
boolean logTransform = optionsPanel.logTransform();
boolean rawCounts = optionsPanel.rawCountsBox.isSelected();
if (rawCounts) {
logTransform = false;
}
Chromosome[] chrs = collection().genome().getAllChromosomes();
for (int c = 0; c < chrs.length; c++) {
// System.err.println("Processing chr "+chrs[c].name());
if (cancel) {
progressCancelled();
return;
}
progressUpdated("Making features for chr" + chrs[c].name(), c, chrs.length * 2);
Feature[] features = collection().genome().annotationCollection().getFeaturesForType(chrs[c], optionsPanel.getSelectedFeatureType());
Arrays.sort(features);
FeatureGroup[] mergedTranscripts = mergeTranscripts(features);
for (int f = 0; f < mergedTranscripts.length; f++) {
if (cancel) {
progressCancelled();
return;
}
probes.add(new Probe(chrs[c], mergedTranscripts[f].start(), mergedTranscripts[f].end(), mergedTranscripts[f].strand(), mergedTranscripts[f].name));
}
}
Probe[] allProbes = probes.toArray(new Probe[0]);
collection().setProbeSet(new ProbeSet("Transcript features over " + optionsPanel.getSelectedFeatureType(), allProbes));
// Having made probes we now need to quantitate them. We'll fetch the
// probes overlapping each sub-feature and then aggregate these together
// to get the final quantitation.
int currentIndex = 0;
for (int c = 0; c < chrs.length; c++) {
if (cancel) {
progressCancelled();
return;
}
progressUpdated("Quantitating features on chr" + chrs[c].name(), chrs.length + c, chrs.length * 2);
Feature[] features = collection().genome().annotationCollection().getFeaturesForType(chrs[c], optionsPanel.getSelectedFeatureType());
Arrays.sort(features);
FeatureGroup[] mergedTranscripts = mergeTranscripts(features);
for (int f = 0; f < mergedTranscripts.length; f++) {
Location[] subLocations = mergedTranscripts[f].getSubLocations();
for (int d = 0; d < data.length; d++) {
if (cancel) {
progressCancelled();
return;
}
int totalCount = 0;
long[] reads = data[d].getReadsForProbe(new Probe(chrs[c], mergedTranscripts[f].start(), mergedTranscripts[f].end()));
READ: for (int r = 0; r < reads.length; r++) {
if (SequenceRead.strand(reads[r]) == mergedTranscripts[f].strand()) {
// TODO: Should we check if we're within the CDS?
++totalCount;
continue;
}
// We'll still count this if it overlaps with an exon
for (int s = 0; s < subLocations.length; s++) {
if (SequenceRead.overlaps(reads[r], subLocations[s].packedPosition())) {
++totalCount;
continue READ;
}
}
}
float value = totalCount;
// If we're log transforming then we need to set zero values to 0.9
if (logTransform && value == 0) {
value = 0.9f;
}
// We also correct by the total read count
if (!rawCounts) {
value /= (data[d].getTotalReadCount() / 1000000f);
}
// Finally we do the log transform if we've been asked to
if (logTransform) {
value = (float) Math.log(value) / log2;
}
data[d].setValueForProbe(allProbes[currentIndex], value);
}
currentIndex++;
}
}
StringBuffer quantitationDescription = new StringBuffer();
quantitationDescription.append("Gene trap quantitation on transcripts of type ");
quantitationDescription.append(optionsPanel.getSelectedFeatureType());
if (logTransform) {
quantitationDescription.append(". Log transformed");
}
collection().probeSet().setCurrentQuantitation(quantitationDescription.toString());
quantitatonComplete();
}
use of uk.ac.babraham.SeqMonk.DataTypes.Genome.Feature in project SeqMonk by s-andrews.
the class SplicingEfficiencyPipeline method startPipeline.
protected void startPipeline() {
// We first need to generate probes over all of the features listed in
// the feature types. The probes should cover the whole area of the
// feature regardless of where it splices.
Vector<Probe> probes = new Vector<Probe>();
boolean correctReadLength = optionsPanel.correctForReadLength();
boolean logTransform = optionsPanel.logTransform();
boolean applyTranscriptLengthCorrection = optionsPanel.applyTranscriptLengthCorrection();
// We need to make a lookup for chromosomes from their names for later on.
Chromosome[] chrs = collection().genome().getAllChromosomes();
Hashtable<String, Chromosome> chrLookup = new Hashtable<String, Chromosome>();
for (int c = 0; c < chrs.length; c++) {
chrLookup.put(chrs[c].name(), chrs[c]);
}
// We're going to cache the gene to transcript mappings so we don't have to re-do them
// later on.
Hashtable<String, Vector<Feature>> transcriptLookup = new Hashtable<String, Vector<Feature>>();
Vector<Feature> usedGeneFeatures = new Vector<Feature>();
// We'll record the number of failures so we can report it
int noTranscriptCount = 0;
int noIntronCount = 0;
for (int c = 0; c < chrs.length; c++) {
if (cancel) {
progressCancelled();
return;
}
progressUpdated("Making features for chr" + chrs[c].name(), c, chrs.length * 2);
Feature[] geneFeatures = collection().genome().annotationCollection().getFeaturesForType(chrs[c], optionsPanel.getSelectedGeneFeatureType());
Arrays.sort(geneFeatures);
Feature[] transcriptFeatures = collection().genome().annotationCollection().getFeaturesForType(chrs[c], optionsPanel.getSelectedTranscriptFeatureType());
Arrays.sort(transcriptFeatures);
// We need to figure out whether we can match based on name, or whether we need to do it based on location.
boolean matchOnNames = true;
for (int t = 0; t < transcriptFeatures.length; t++) {
if (!transcriptFeatures[t].name().matches("-\\d\\d\\d$")) {
matchOnNames = false;
break;
}
}
if (!matchOnNames) {
progressWarningReceived(new SeqMonkException("Non-Ensembl transcript names found - matching based on position"));
}
if (matchOnNames) {
for (int i = 0; i < transcriptFeatures.length; i++) {
String name = transcriptFeatures[i].name();
name = name.replaceAll("-\\d\\d\\d$", "");
if (!transcriptLookup.containsKey(name)) {
transcriptLookup.put(name, new Vector<Feature>());
}
transcriptLookup.get(name).add(transcriptFeatures[i]);
}
} else {
// We need to go through the genes and transcripts in parallel and match them up based on containment
// and direction.
int lastGeneIndex = 0;
for (int t = 0; t < transcriptFeatures.length; t++) {
for (int g = lastGeneIndex; g < geneFeatures.length; g++) {
if (transcriptFeatures[t].location().start() > geneFeatures[g].location().end()) {
lastGeneIndex = g;
continue;
}
// If the gene is beyond the end of the transcript we can stop looking
if (geneFeatures[g].location().start() > transcriptFeatures[t].location().end()) {
break;
}
// If we're on the same strand and contained within the gene then we have a match
if (geneFeatures[g].location().strand() == transcriptFeatures[t].location().strand() && transcriptFeatures[t].location().start() >= geneFeatures[g].location().start() && transcriptFeatures[t].location().end() <= geneFeatures[g].location().end()) {
String name = geneFeatures[g].name();
if (!transcriptLookup.containsKey(name)) {
transcriptLookup.put(name, new Vector<Feature>());
}
transcriptLookup.get(name).add(transcriptFeatures[t]);
}
}
}
}
for (int f = 0; f < geneFeatures.length; f++) {
if (cancel) {
progressCancelled();
return;
}
if (!transcriptLookup.containsKey(geneFeatures[f].name())) {
++noTranscriptCount;
// No good making features for genes with no transcript.
continue;
}
Feature[] localTranscripts = transcriptLookup.get(geneFeatures[f].name()).toArray(new Feature[0]);
boolean[] exonPositions = new boolean[geneFeatures[f].location().length()];
for (int t = 0; t < localTranscripts.length; t++) {
if (!SequenceRead.overlaps(geneFeatures[f].location().packedPosition(), localTranscripts[t].location().packedPosition())) {
// It's not a match for this gene
continue;
}
if (!(localTranscripts[t].location() instanceof SplitLocation)) {
for (int p = localTranscripts[t].location().start() - geneFeatures[f].location().start(); p <= localTranscripts[t].location().end() - geneFeatures[f].location().start(); p++) {
if (p < 0)
continue;
if (p >= exonPositions.length)
continue;
exonPositions[p] = true;
}
continue;
}
Location[] exons = ((SplitLocation) localTranscripts[t].location()).subLocations();
for (int e = 0; e < exons.length; e++) {
for (int p = exons[e].start() - geneFeatures[f].location().start(); p <= exons[e].end() - geneFeatures[f].location().start(); p++) {
if (p < 0)
continue;
if (p >= exonPositions.length)
continue;
exonPositions[p] = true;
}
}
}
int exonPositionCount = 0;
int intronPositionCount = 0;
for (int p = 0; p < exonPositions.length; p++) {
if (exonPositions[p])
exonPositionCount++;
else
intronPositionCount++;
}
if (exonPositionCount == 0) {
++noTranscriptCount;
continue;
}
if (intronPositionCount == 0) {
++noIntronCount;
continue;
}
probes.add(new Probe(chrs[c], geneFeatures[f].location().start(), geneFeatures[f].location().end(), geneFeatures[f].location().strand(), geneFeatures[f].name()));
usedGeneFeatures.add(geneFeatures[f]);
}
}
if (noTranscriptCount > 0) {
progressWarningReceived(new SeqMonkException("" + noTranscriptCount + " genes had no associated transcripts"));
}
if (noIntronCount > 0) {
progressWarningReceived(new SeqMonkException("" + noIntronCount + " genes had no intron only areas"));
}
Probe[] allProbes = probes.toArray(new Probe[0]);
collection().setProbeSet(new ProbeSet("Gene features over " + optionsPanel.getSelectedGeneFeatureType(), allProbes));
// Having made probes we now need to quantitate them. We'll get the full set of reads for the full gene
// and then work out a combined set of exons. We can then quantitate the bases which fall into introns
// and exons separately and then work out a ratio from there.
QuantitationStrandType readFilter = optionsPanel.readFilter();
Feature[] geneFeatures = usedGeneFeatures.toArray(new Feature[0]);
for (int g = 0; g < geneFeatures.length; g++) {
if (cancel) {
progressCancelled();
return;
}
if (g % 100 == 0) {
progressUpdated("Quantitating features", g, geneFeatures.length);
}
int[] readLengths = new int[data.length];
if (correctReadLength) {
for (int d = 0; d < data.length; d++) {
readLengths[d] = data[d].getMaxReadLength();
}
}
// Find the set of transcripts which relate to this gene.
Feature[] transcripts = transcriptLookup.get(geneFeatures[g].name()).toArray(new Feature[0]);
Vector<Feature> validTranscripts = new Vector<Feature>();
for (int t = 0; t < transcripts.length; t++) {
if (SequenceRead.overlaps(geneFeatures[g].location().packedPosition(), transcripts[t].location().packedPosition())) {
validTranscripts.add(transcripts[t]);
}
}
transcripts = validTranscripts.toArray(new Feature[0]);
// before so if we do then something has gone wrong.
if (transcripts.length == 0) {
throw new IllegalStateException("No transcripts for gene " + geneFeatures[g] + " on chr " + geneFeatures[g].chromosomeName());
}
Vector<Location> allExonsVector = new Vector<Location>();
for (int t = 0; t < transcripts.length; t++) {
if (transcripts[t].location() instanceof SplitLocation) {
Location[] sublocs = ((SplitLocation) transcripts[t].location()).subLocations();
for (int i = 0; i < sublocs.length; i++) allExonsVector.add(sublocs[i]);
} else {
allExonsVector.add(transcripts[t].location());
}
}
// if (geneFeatures[f].name().equals("Cpa6")) {
// System.err.println("Cpa6 had "+allExonsVector.size()+" total exons");
// }
Collections.sort(allExonsVector);
// Now go through and make a merged set of exons to remove any redundancy
Vector<Location> mergedExonsVector = new Vector<Location>();
Location lastLocation = null;
Enumeration<Location> en = allExonsVector.elements();
while (en.hasMoreElements()) {
Location l = en.nextElement();
if (lastLocation == null) {
// if (geneFeatures[f].name().equals("Cpa6")) {
// System.err.println("Setting as first location");
// }
lastLocation = l;
continue;
}
// Check if it's the same as the last, which is likely
if (l.start() == lastLocation.start() && l.end() == lastLocation.end()) {
continue;
}
// Check if they overlap and can be merged
if (l.start() <= lastLocation.end() && l.end() >= lastLocation.start()) {
// if (geneFeatures[f].name().equals("Cpa6")) {
// System.err.println("Overlaps with last location");
// }
// It overlaps with the last location so merge them
lastLocation = new Location(Math.min(l.start(), lastLocation.start()), Math.max(l.end(), lastLocation.end()), geneFeatures[g].location().strand());
// }
continue;
}
// Start a new location
// if (geneFeatures[f].name().equals("Cpa6")) {
// System.err.println("Doesn't overlap - adding last location and creating new one");
// }
mergedExonsVector.add(lastLocation);
lastLocation = l;
}
if (lastLocation != null) {
mergedExonsVector.add(lastLocation);
}
// if (geneFeatures[f].name().equals("Cpa6")) {
// System.err.println("Cpa6 had "+mergedExonsVector.size()+" merged exons");
// }
// Now we can start the quantitation.
int intronLength = geneFeatures[g].location().length();
int exonLength = 0;
Location[] subLocs = mergedExonsVector.toArray(new Location[0]);
for (int l = 0; l < subLocs.length; l++) {
exonLength += subLocs[l].length();
}
// if (geneFeatures[f].name().equals("Cpa6")) {
// System.err.println("Cpa6 total intron length="+intronLength+" exon length="+exonLength);
// }
intronLength -= exonLength;
if (intronLength <= 0) {
progressWarningReceived(new IllegalStateException("Intron length of " + intronLength + " for gene " + geneFeatures[g]));
for (int d = 0; d < data.length; d++) {
data[d].setValueForProbe(allProbes[g], Float.NaN);
}
continue;
}
if (exonLength <= 0) {
progressWarningReceived(new IllegalStateException("Exon length of " + exonLength + " for gene " + geneFeatures[g]));
for (int d = 0; d < data.length; d++) {
data[d].setValueForProbe(allProbes[g], Float.NaN);
}
continue;
}
for (int d = 0; d < data.length; d++) {
if (cancel) {
progressCancelled();
return;
}
int totalIntronCount = 0;
int totalExonCount = 0;
long[] reads = data[d].getReadsForProbe(new Probe(chrLookup.get(geneFeatures[g].chromosomeName()), geneFeatures[g].location().start(), geneFeatures[g].location().end()));
for (int r = 0; r < reads.length; r++) {
if (!readFilter.useRead(geneFeatures[g].location(), reads[r])) {
continue;
}
int overlap = (Math.min(geneFeatures[g].location().end(), SequenceRead.end(reads[r])) - Math.max(geneFeatures[g].location().start(), SequenceRead.start(reads[r]))) + 1;
totalIntronCount += overlap;
// Now we see if we overlap with any of the exons
for (int s = 0; s < subLocs.length; s++) {
if (subLocs[s].start() <= SequenceRead.end(reads[r]) && subLocs[s].end() >= SequenceRead.start(reads[r])) {
int exonOverlap = (Math.min(subLocs[s].end(), SequenceRead.end(reads[r])) - Math.max(subLocs[s].start(), SequenceRead.start(reads[r]))) + 1;
if (exonOverlap > 0) {
totalExonCount += exonOverlap;
totalIntronCount -= exonOverlap;
}
}
}
}
if (totalIntronCount < 0) {
progressWarningReceived(new SeqMonkException("Intron count of " + totalIntronCount + " for " + geneFeatures[g].name()));
continue;
}
// reads which this count could comprise, rounding down to a whole number.
if (correctReadLength) {
totalIntronCount /= readLengths[d];
totalExonCount /= readLengths[d];
}
float intronValue = totalIntronCount;
float exonValue = totalExonCount;
// If we're log transforming then we need to set zero values to 0.9
if (logTransform && intronValue == 0) {
intronValue = 0.9f;
}
if (logTransform && exonValue == 0) {
exonValue = 0.9f;
}
// been asked to.
if (applyTranscriptLengthCorrection) {
intronValue /= (intronLength / 1000f);
exonValue /= (exonLength / 1000f);
}
// We also correct by the total read count, or length
if (correctReadLength) {
intronValue /= (data[d].getTotalReadCount() / 1000000f);
exonValue /= (data[d].getTotalReadCount() / 1000000f);
} else {
intronValue /= (data[d].getTotalReadLength() / 1000000f);
exonValue /= (data[d].getTotalReadLength() / 1000000f);
}
// Finally we do the log transform if we've been asked to
if (logTransform) {
if (intronValue == 0) {
intronValue = 0.0001f;
}
if (exonValue == 0) {
exonValue = 0.0001f;
}
intronValue = (float) Math.log(intronValue) / log2;
exonValue = (float) Math.log(exonValue) / log2;
}
// Now we check what value they actually wanted to record
switch(optionsPanel.quantitationType()) {
case (EXONS):
{
data[d].setValueForProbe(allProbes[g], exonValue);
break;
}
case (INTRONS):
{
data[d].setValueForProbe(allProbes[g], intronValue);
break;
}
case (RATIO):
{
if (logTransform) {
data[d].setValueForProbe(allProbes[g], exonValue - intronValue);
} else {
if (intronValue == 0) {
intronValue = 0.0001f;
}
data[d].setValueForProbe(allProbes[g], exonValue / intronValue);
}
break;
}
default:
throw new IllegalStateException("Unknonwn quantitation type " + optionsPanel.quantitationType());
}
// if (geneFeatures[f].name().equals("Cpa6")) {
// try {
// System.err.println("Stored value was "+data[d].getValueForProbe(allProbes[currentIndex]));
// } catch (SeqMonkException e) {
// e.printStackTrace();
// }
// }
}
}
StringBuffer quantitationDescription = new StringBuffer();
quantitationDescription.append("Splicing efficiency quantitation");
if (correctReadLength) {
quantitationDescription.append(" counting reads");
} else {
quantitationDescription.append(" counting bases");
}
if (optionsPanel.logTransform()) {
quantitationDescription.append(" log transformed");
}
if (applyTranscriptLengthCorrection) {
quantitationDescription.append(" correcting for feature length");
}
// TODO: Add more description
collection().probeSet().setCurrentQuantitation(quantitationDescription.toString());
quantitatonComplete();
}
use of uk.ac.babraham.SeqMonk.DataTypes.Genome.Feature in project SeqMonk by s-andrews.
the class AntisenseTranscriptionPipeline method startPipeline.
protected void startPipeline() {
// We first need to generate probes over all of the features listed in
// the feature types. The probes should cover the whole area of the
// feature regardless of where it splices.
Vector<Probe> probes = new Vector<Probe>();
double pValue = optionsPanel.pValue();
QuantitationStrandType readFilter = optionsPanel.readFilter();
long[] senseCounts = new long[data.length];
long[] antisenseCounts = new long[data.length];
Chromosome[] chrs = collection().genome().getAllChromosomes();
for (int c = 0; c < chrs.length; c++) {
if (cancel) {
progressCancelled();
return;
}
progressUpdated("Getting total antisense rate for chr" + chrs[c].name(), c, chrs.length * 2);
Feature[] features = getValidFeatures(chrs[c]);
for (int f = 0; f < features.length; f++) {
if (cancel) {
progressCancelled();
return;
}
Probe p = new Probe(chrs[c], features[f].location().start(), features[f].location().end(), features[f].location().strand(), features[f].name());
probes.add(p);
for (int d = 0; d < data.length; d++) {
long[] reads = data[d].getReadsForProbe(p);
for (int r = 0; r < reads.length; r++) {
if (readFilter.useRead(p, reads[r])) {
senseCounts[d] += SequenceRead.length(reads[r]);
} else {
antisenseCounts[d] += SequenceRead.length(reads[r]);
}
}
}
}
}
Probe[] allProbes = probes.toArray(new Probe[0]);
collection().setProbeSet(new ProbeSet("Features over " + optionsPanel.getSelectedFeatureType(), allProbes));
// Now we can work out the overall antisense rate
double[] antisenseProbability = new double[data.length];
for (int d = 0; d < data.length; d++) {
System.err.println("Antisense counts are " + antisenseCounts[d] + " sense counts are " + senseCounts[d]);
antisenseProbability[d] = antisenseCounts[d] / (double) (antisenseCounts[d] + senseCounts[d]);
System.err.println("Antisense probability for " + data[d].name() + " is " + antisenseProbability[d]);
}
// Now we can quantitate each individual feature and test for whether it is significantly
// showing antisense expression
ArrayList<Vector<ProbeTTestValue>> significantProbes = new ArrayList<Vector<ProbeTTestValue>>();
for (int d = 0; d < data.length; d++) {
significantProbes.add(new Vector<ProbeTTestValue>());
}
int[] readLengths = new int[data.length];
for (int d = 0; d < readLengths.length; d++) {
readLengths[d] = data[d].getMaxReadLength();
System.err.println("For " + data[d].name() + " max read len is " + readLengths[d]);
}
for (int c = 0; c < chrs.length; c++) {
if (cancel) {
progressCancelled();
return;
}
progressUpdated("Quantitating features on chr" + chrs[c].name(), chrs.length + c, chrs.length * 2);
Probe[] thisChrProbes = collection().probeSet().getProbesForChromosome(chrs[c]);
for (int p = 0; p < thisChrProbes.length; p++) {
for (int d = 0; d < data.length; d++) {
if (cancel) {
progressCancelled();
return;
}
long senseCount = 0;
long antisenseCount = 0;
long[] reads = data[d].getReadsForProbe(thisChrProbes[p]);
for (int r = 0; r < reads.length; r++) {
if (readFilter.useRead(thisChrProbes[p], reads[r])) {
// TODO: Just count overlap?
senseCount += SequenceRead.length(reads[r]);
} else {
antisenseCount += SequenceRead.length(reads[r]);
}
}
// if (thisChrProbes[p].name().equals("RP4-798A10.2")) {
// System.err.println("Raw base counts are sense="+senseCount+" anti="+antisenseCount+" from "+reads.length+" reads");
// }
int senseReads = (int) (senseCount / readLengths[d]);
int antisenseReads = (int) (antisenseCount / readLengths[d]);
// if (thisChrProbes[p].name().equals("RP4-798A10.2")) {
// System.err.println("Raw read counts are sense="+senseReads+" anti="+antisenseReads+" from "+reads.length+" reads");
// }
BinomialDistribution bd = new BinomialDistribution(senseReads + antisenseReads, antisenseProbability[d]);
// Since the binomial distribution gives the probability of getting a value higher than
// this we need to subtract one so we get the probability of this or higher.
double thisPValue = 1 - bd.cumulativeProbability(antisenseReads - 1);
if (antisenseReads == 0)
thisPValue = 1;
// We have to add all results at this stage so we don't mess up the multiple
// testing correction later on.
significantProbes.get(d).add(new ProbeTTestValue(thisChrProbes[p], thisPValue));
double expected = ((senseReads + antisenseReads) * antisenseProbability[d]);
if (expected < 1)
expected = 1;
float obsExp = antisenseReads / (float) expected;
data[d].setValueForProbe(thisChrProbes[p], obsExp);
}
}
}
// filtering those which pass our p-value cutoff
for (int d = 0; d < data.length; d++) {
ProbeTTestValue[] ttestResults = significantProbes.get(d).toArray(new ProbeTTestValue[0]);
BenjHochFDR.calculateQValues(ttestResults);
ProbeList newList = new ProbeList(collection().probeSet(), "Antisense < " + pValue + " in " + data[d].name(), "Probes showing significant antisense transcription from a basal level of " + antisenseProbability[d] + " with a cutoff of " + pValue, "FDR");
for (int i = 0; i < ttestResults.length; i++) {
if (ttestResults[i].probe.name().equals("RP4-798A10.2")) {
System.err.println("Raw p=" + ttestResults[i].p + " q=" + ttestResults[i].q);
}
if (ttestResults[i].q < pValue) {
newList.addProbe(ttestResults[i].probe, (float) ttestResults[i].q);
}
}
}
StringBuffer quantitationDescription = new StringBuffer();
quantitationDescription.append("Antisense transcription pipeline quantitation ");
quantitationDescription.append(". Directionality was ");
quantitationDescription.append(optionsPanel.libraryTypeBox.getSelectedItem());
if (optionsPanel.ignoreOverlaps()) {
quantitationDescription.append(". Ignoring existing overlaps");
}
quantitationDescription.append(". P-value cutoff was ");
quantitationDescription.append(optionsPanel.pValue());
collection().probeSet().setCurrentQuantitation(quantitationDescription.toString());
quantitatonComplete();
}
use of uk.ac.babraham.SeqMonk.DataTypes.Genome.Feature in project SeqMonk by s-andrews.
the class IntronRegressionPipeline method startPipeline.
protected void startPipeline() {
// We first need to generate probes over all of the features listed in
// the feature types. The probes should cover the whole area of the
// feature regardless of where it splices.
Vector<Probe> probes = new Vector<Probe>();
int minDensity = optionsPanel.minDensity();
int minLength = optionsPanel.minLength();
double maxPValue = optionsPanel.maxPValue();
int binSize = optionsPanel.measurementBinSize();
QuantitationStrandType readFilter = optionsPanel.readFilter();
Chromosome[] chrs = collection().genome().getAllChromosomes();
for (int c = 0; c < chrs.length; c++) {
if (cancel) {
progressCancelled();
return;
}
Vector<Probe> probesForThisChromosome = new Vector<Probe>();
progressUpdated("Making probes", c, chrs.length);
Feature[] features = getValidFeatures(chrs[c]);
for (int f = 0; f < features.length; f++) {
if (cancel) {
progressCancelled();
return;
}
// There are no introns here
if (!(features[f].location() instanceof SplitLocation))
continue;
Location[] subLocations = ((SplitLocation) features[f].location()).subLocations();
// TODO: Reverse the subLocations if its a reverse feature
for (int intron = 1; intron < subLocations.length; intron++) {
int start = subLocations[intron - 1].end();
int end = subLocations[intron].start();
if ((end - start) + 1 < minLength) {
// This intron is too short.
continue;
}
// TODO: We could throw away any probes which didn't have enough reads in any feature
Probe p = new Probe(chrs[c], start, end, features[f].location().strand(), features[f].name() + "_" + intron);
probesForThisChromosome.add(p);
}
}
// Now we can deduplicate the probes for this chromosome and add them to the main collection
Probe[] dupProbes = probesForThisChromosome.toArray(new Probe[0]);
Arrays.sort(dupProbes);
for (int p = 0; p < dupProbes.length; p++) {
if (p > 0 && dupProbes[p].packedPosition() == dupProbes[p - 1].packedPosition())
continue;
probes.add(dupProbes[p]);
}
}
Probe[] allProbes = probes.toArray(new Probe[0]);
collection().setProbeSet(new ProbeSet("Features over " + optionsPanel.getSelectedFeatureType(), allProbes));
// Now we go back through the probes and quantitate them
for (int p = 0; p < allProbes.length; p++) {
if (cancel) {
progressCancelled();
return;
}
if (p % 1000 == 0) {
progressUpdated("Quantitated " + p + " out of " + allProbes.length + " probes", p, allProbes.length);
}
for (int d = 0; d < data.length; d++) {
long[] reads = data[d].getReadsForProbe(allProbes[p]);
int[] countsPerSite = new int[allProbes[p].length()];
int usableCounts = 0;
for (int r = 0; r < reads.length; r++) {
if (readFilter.useRead(allProbes[p], reads[r])) {
++usableCounts;
for (int pos = Math.max(0, SequenceRead.start(reads[r]) - allProbes[p].start()); pos <= Math.min(countsPerSite.length - 1, SequenceRead.end(reads[r]) - allProbes[p].start()); pos++) {
++countsPerSite[pos];
}
}
}
if (usableCounts / (allProbes[p].length() / 1000d) >= minDensity) {
// We're going to do a linear regression rather than a correlation
// We're analysing in bins so we'll work out the bin counts and
// add them dynamically to the regression.
SimpleRegression regression = new SimpleRegression();
int binCount = 0;
for (int i = 0; i < countsPerSite.length; i++) {
if (i > 0 && i % binSize == 0) {
regression.addData(i, binCount);
binCount = 0;
}
binCount += countsPerSite[i];
}
float slope = (float) (regression.getSlope() * 1000000);
double pValue = regression.getSignificance();
if (allProbes[p].strand() == Location.REVERSE) {
slope = 0 - slope;
}
if (pValue <= maxPValue) {
data[d].setValueForProbe(allProbes[p], slope);
} else {
data[d].setValueForProbe(allProbes[p], Float.NaN);
}
} else {
data[d].setValueForProbe(allProbes[p], Float.NaN);
}
}
}
StringBuffer quantitationDescription = new StringBuffer();
quantitationDescription.append("Intron regression pipeline quantitation ");
quantitationDescription.append(". Directionality was ");
quantitationDescription.append(optionsPanel.libraryTypeBox.getSelectedItem());
quantitationDescription.append(". Min intron length was ");
quantitationDescription.append(minLength);
quantitationDescription.append(". Min read density was ");
quantitationDescription.append(minDensity);
quantitationDescription.append(". Max slope p-value was ");
quantitationDescription.append(maxPValue);
collection().probeSet().setCurrentQuantitation(quantitationDescription.toString());
quantitatonComplete();
}
Aggregations