use of dr.evolution.alignment.ExtractPairs in project beast-mcmc by beast-dev.
the class AlignmentScore method main.
public static void main(String[] args) throws java.io.IOException, Importer.ImportException {
NexusImporter importer = new NexusImporter(new FileReader(args[0]));
Alignment alignment = importer.importAlignment();
ExtractPairs pairs = new ExtractPairs(alignment);
Parameter muParam = new Parameter.Default(1.0);
Parameter kappaParam = new Parameter.Default(1.0);
kappaParam.addBounds(new Parameter.DefaultBounds(100.0, 0.0, 1));
muParam.addBounds(new Parameter.DefaultBounds(1.0, 1.0, 1));
Parameter freqParam = new Parameter.Default(alignment.getStateFrequencies());
FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, freqParam);
SubstitutionModel substModel = new HKY(kappaParam, freqModel);
SiteModel siteModel = new GammaSiteModel(substModel, muParam, null, 1, null);
ScoreMatrix scoreMatrix = new ScoreMatrix(siteModel, 0.1);
double threshold = 0.1;
List<PairDistance> pairDistances = new ArrayList<PairDistance>();
Set<Integer> sequencesUsed = new HashSet<Integer>();
List<Integer> allGaps = new ArrayList<Integer>();
for (int i = 0; i < alignment.getSequenceCount(); i++) {
for (int j = i + 1; j < alignment.getSequenceCount(); j++) {
Alignment pairAlignment = pairs.getPairAlignment(i, j);
if (pairAlignment != null) {
SitePatterns patterns = new SitePatterns(pairAlignment);
double distance = getGeneticDistance(scoreMatrix, patterns);
if (distance < threshold) {
List gaps = new ArrayList();
GapUtils.getGapSizes(pairAlignment, gaps);
pairDistances.add(new PairDistance(i, j, distance, gaps, pairAlignment.getSiteCount()));
System.out.print(".");
} else {
System.out.print("*");
}
} else {
System.out.print("x");
}
}
System.out.println();
}
Collections.sort(pairDistances);
int totalLength = 0;
for (PairDistance pairDistance : pairDistances) {
Integer x = pairDistance.x;
Integer y = pairDistance.y;
if (!sequencesUsed.contains(x) && !sequencesUsed.contains(y)) {
allGaps.addAll(pairDistance.gaps);
sequencesUsed.add(x);
sequencesUsed.add(y);
System.out.println("Added pair (" + x + "," + y + ") d=" + pairDistance.distance + " L=" + pairDistance.alignmentLength);
totalLength += pairDistance.alignmentLength;
}
}
printFrequencyTable(allGaps);
System.out.println("total length=" + totalLength);
}
Aggregations