use of org.apache.commons.math3.geometry.partitioning.Region.Location in project SeqMonk by s-andrews.
the class CodonBiasPipeline 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();
String libraryType = optionsPanel.libraryType();
Chromosome[] chrs = collection().genome().getAllChromosomes();
for (int c = 0; c < chrs.length; c++) {
if (cancel) {
progressCancelled();
return;
}
progressUpdated("Making probes for chr" + chrs[c].name(), c, chrs.length * 2);
Feature[] features = collection().genome().annotationCollection().getFeaturesForType(chrs[c], optionsPanel.getSelectedFeatureType());
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);
}
}
allProbes = probes.toArray(new Probe[0]);
collection().setProbeSet(new ProbeSet("Features over " + optionsPanel.getSelectedFeatureType(), allProbes));
// Now we can quantitate each individual feature and test for whether it is significantly
// showing codon bias
ArrayList<Vector<ProbeTTestValue>> significantProbes = new ArrayList<Vector<ProbeTTestValue>>();
// data contains the data stores that this pipeline is going to use. We need to test each data store.
for (int d = 0; d < data.length; d++) {
significantProbes.add(new Vector<ProbeTTestValue>());
}
int probeCounter = 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());
for (int p = 0; p < features.length; p++) {
// Get the corresponding feature and work out the mapping between genomic position and codon sub position.
int[] mappingArray = createGenomeMappingArray(features[p]);
DATASTORE_LOOP: for (int d = 0; d < data.length; d++) {
if (cancel) {
progressCancelled();
return;
}
long[] reads = data[d].getReadsForProbe(allProbes[probeCounter]);
// TODO: make this configurable
if (reads.length < 5) {
data[d].setValueForProbe(allProbes[probeCounter], Float.NaN);
continue DATASTORE_LOOP;
}
int pos1Count = 0;
int pos2Count = 0;
int pos3Count = 0;
READ_LOOP: for (int r = 0; r < reads.length; r++) {
int genomicReadStart = SequenceRead.start(reads[r]);
int genomicReadEnd = SequenceRead.end(reads[r]);
int readStrand = SequenceRead.strand(reads[r]);
int relativeReadStart = -1;
// forward reads
if (readStrand == 1) {
if (libraryType == "Same strand specific") {
if (features[p].location().strand() == Location.FORWARD) {
// The start of the read needs to be within the feature
if (genomicReadStart - features[p].location().start() < 0) {
continue READ_LOOP;
} else {
// look up the read start pos in the mapping array
relativeReadStart = mappingArray[genomicReadStart - features[p].location().start()];
}
}
} else if (libraryType == "Opposing strand specific") {
if (features[p].location().strand() == Location.REVERSE) {
// The "start" of a reverse read/probe is actually the end
if (features[p].location().end() - genomicReadEnd < 0) {
continue READ_LOOP;
} else {
relativeReadStart = mappingArray[features[p].location().end() - genomicReadEnd];
}
}
}
}
// reverse reads
if (readStrand == -1) {
if (libraryType == "Same strand specific") {
if (features[p].location().strand() == Location.REVERSE) {
if (features[p].location().end() - genomicReadEnd < 0) {
continue READ_LOOP;
} else {
// System.out.println("features[p].location().end() is " + features[p].location().end() + ", genomicReadEnd is " + genomicReadEnd);
// System.out.println("mapping array[0] is " + mappingArray[0]);
relativeReadStart = mappingArray[features[p].location().end() - genomicReadEnd];
}
}
} else if (libraryType == "Opposing strand specific") {
if (features[p].location().strand() == Location.FORWARD) {
// The start of the read needs to be within the feature
if (genomicReadStart - features[p].location().start() < 0) {
continue READ_LOOP;
} else {
// look up the read start position in the mapping array
relativeReadStart = mappingArray[genomicReadStart - features[p].location().start()];
}
}
}
}
// find out which position the read is in
if (relativeReadStart == -1) {
continue READ_LOOP;
} else if (relativeReadStart % 3 == 0) {
pos3Count++;
continue READ_LOOP;
} else if ((relativeReadStart + 1) % 3 == 0) {
pos2Count++;
continue READ_LOOP;
} else if ((relativeReadStart + 2) % 3 == 0) {
pos1Count++;
}
}
// closing bracket for read loop
// System.out.println("pos1Count for "+ features[p].name() + " is " + pos1Count);
// System.out.println("pos2Count for "+ features[p].name() + " is " + pos2Count);
// System.out.println("pos3Count for "+ features[p].name() + " is " + pos3Count);
int interestingCodonCount = 0;
int otherCodonCount = 0;
if (optionsPanel.codonSubPosition() == 1) {
interestingCodonCount = pos1Count;
otherCodonCount = pos2Count + pos3Count;
} else if (optionsPanel.codonSubPosition() == 2) {
interestingCodonCount = pos2Count;
otherCodonCount = pos1Count + pos3Count;
} else if (optionsPanel.codonSubPosition() == 3) {
interestingCodonCount = pos3Count;
otherCodonCount = pos1Count + pos2Count;
}
int totalCount = interestingCodonCount + otherCodonCount;
// BinomialDistribution bd = new BinomialDistribution(interestingCodonCount+otherCodonCount, 1/3d);
BinomialDistribution bd = new BinomialDistribution(totalCount, 1 / 3d);
// 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(interestingCodonCount - 1);
if (interestingCodonCount == 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(allProbes[probeCounter], thisPValue));
float percentageCount;
if (totalCount == 0) {
percentageCount = 0;
} else {
percentageCount = ((float) interestingCodonCount / (float) totalCount) * 100;
}
data[d].setValueForProbe(allProbes[probeCounter], percentageCount);
// System.out.println("totalCount = " + totalCount);
// System.out.println("interestingCodonCount " + interestingCodonCount);
// System.out.println("pValue = " + thisPValue);
// System.out.println("percentageCount = " + percentageCount);
// System.out.println("");
}
probeCounter++;
}
}
// 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(), "Codon bias < " + pValue + " in " + data[d].name(), "Probes showing significant codon bias for position " + optionsPanel.codonSubPosition() + " with a cutoff of " + pValue, "FDR");
for (int i = 0; i < ttestResults.length; i++) {
if (ttestResults[i].q < pValue) {
newList.addProbe(ttestResults[i].probe, (float) ttestResults[i].q);
}
}
}
StringBuffer quantitationDescription = new StringBuffer();
quantitationDescription.append("Codon bias pipeline using codon position " + optionsPanel.codonSubPosition() + " for " + optionsPanel.libraryType() + " library.");
collection().probeSet().setCurrentQuantitation(quantitationDescription.toString());
quantitatonComplete();
}
use of org.apache.commons.math3.geometry.partitioning.Region.Location in project j6dof-flight-sim by chris-ali.
the class Engine method calculateEngMoments.
/**
* Calculates the moment generated by the engine as a function of its thrust and location
* relative to the aircraft's center of gravity. Used in {@link Engine#updateEngineState(EnumMap, EnumMap, double[])}
*/
protected void calculateEngMoments() {
Vector3D forceVector = new Vector3D(engineThrust);
Vector3D armVector = new Vector3D(enginePosition);
this.engineMoment = Vector3D.crossProduct(forceVector, armVector).toArray();
}
use of org.apache.commons.math3.geometry.partitioning.Region.Location in project MindsEye by SimiaCryptus.
the class ObjectLocation method renderAlpha.
/**
* Render alpha tensor.
*
* @param alphaPower the alpha power
* @param img the img
* @param locationResult the location result
* @param classification the classification
* @param category the category
* @return the tensor
*/
public Tensor renderAlpha(final double alphaPower, final Tensor img, final Result locationResult, final Tensor classification, final int category) {
TensorArray tensorArray = TensorArray.wrap(new Tensor(classification.getDimensions()).set(category, 1));
DeltaSet<Layer> deltaSet = new DeltaSet<>();
locationResult.accumulate(deltaSet, tensorArray);
double[] rawDelta = deltaSet.getMap().entrySet().stream().filter(x -> x.getValue().target == img.getData()).findAny().get().getValue().getDelta();
Tensor deltaColor = new Tensor(rawDelta, img.getDimensions()).mapAndFree(x -> Math.abs(x));
Tensor delta1d = blur(reduce(deltaColor), 3);
return TestUtil.normalizeBands(TestUtil.normalizeBands(delta1d, 1).mapAndFree(x -> Math.pow(x, alphaPower)));
}
use of org.apache.commons.math3.geometry.partitioning.Region.Location in project chordatlas by twak.
the class Concarnie method removeOverlaps.
private void removeOverlaps(List<Problem> current) {
Closer<Problem> closer = new Closer();
for (Problem a : current) {
try {
Region<Euclidean2D> ar = a.chull.createRegion();
b: for (Problem b : current) if (a != b)
for (Vector2D v : b.chull.getVertices()) {
Location vInA = ar.checkPoint(v);
if (vInA == Location.BOUNDARY || vInA == Location.INSIDE) {
closer.add(a, b);
continue b;
}
}
} catch (InsufficientDataException th) {
} catch (MathIllegalArgumentException th) {
}
}
for (Set<Problem> close : closer.close()) {
List<Problem> intersecting = new ArrayList<Problem>(close);
Collections.sort(intersecting, (a, b) -> -Double.compare(a.area(), b.area()));
for (int i = 1; i < intersecting.size(); i++) {
Problem togo = intersecting.get(i);
togo.addPortal();
current.remove(togo);
}
}
}
use of org.apache.commons.math3.geometry.partitioning.Region.Location 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();
}
Aggregations