BioJava:CookbookFrench:Distribution:Gibbs
Comment construire un échantillonneur de Gibbs à l’aide de Distributions?
L’échantillonnage de Gibbs est une technique statistique apparentée à l’échantillonnage Monte Carlo et chaînes de Markov. On l’utilise afin de trouver un ensemble de solutions à une optimisation ou au moins une solution qui soit optimale dans un espace local. C’est essentiellement une technique itérative: un seul paramètre est sélectionné aléatoirement et sa valeur fixée également de manière aléatoire (ou à partir d’une distribution de valeurs possibles) alors que tous les autres paramètres demeurent inchangés. si la nouvelle solution est meilleure que l’ancienne, celle-ci devient le nouveau modèle; sinon, l’ancien est conservé. Le processus de sélection des paramètres et de leurs valeurs continue jusqu’à ce qu’une certaine valeur-seuil soit atteinte comme par exemple, la convergence du modèle vers une solution optimale localement ou lorsqu’un certina nombre d’itérations ont été effectués. En biologie, l’échantillonnage de Gibbs a été appliqué avec succès pour des tâches tel que la découverte de motifs conservés dans de grandes séquences. On appelle également cette technique l’alignement de Gibbs.
Il est très facile de construire un simple automate d’alignement de Gibbs en utilisant le package org.biojava.bio.dist de BioJava. C’est également une excellente opportunité d’explorer certaines des classes de la famille Distribution. Dans le code de démonstration ci-dessous, des Distributions sont utilisés afin de randomiser les écarts (offsets) d’alignement et pour calculer le contenu en information d’un motif. Le premier exemple peut paraitre surprenant parce que lDistribution se fait sur un alphabet d’entiers; le deuxième emploit un alphabet d’ADN ou de protéines. Ceci démontre qu’il est assez simple d’utiliser et d’échantillonner une Distribution sur n’importe quel Alphabet pouvant être construit avec BioJava. Dans une tel cas, BioJava n’est pas simplement ‘bio’ mais peut être utilisé afi nde représenter et de manipuler n’importe quelle donnée symbolique.
La première classe se nomme SimpleGibbsAligner. C’est le moteur de base, faisant tout le travail d’échantillonnage et d’évaluation des motifs. Elle utilise une interface, GibbsStoppingCriteria, qui qui l’assiste en figurant quand arrêter l’itération. L’interface présenté fournit égalemnt quelques implémentations simples. Finalement, une application de démonstration avec la méthode main() assemble le tout pour effectuer le travail à la console.
SimpleGibbsAligner
```java package gibbs;
import java.util.HashMap; import java.util.Map; import java.util.Random; import java.util.Vector; import org.biojava.bio.BioException; import org.biojava.bio.dist.Distribution; import org.biojava.bio.dist.DistributionFactory; import org.biojava.bio.dist.DistributionTools; import org.biojava.bio.dist.DistributionTrainerContext; import org.biojava.bio.dist.SimpleDistributionTrainerContext; import org.biojava.bio.seq.Sequence; import org.biojava.bio.seq.SequenceIterator; import org.biojava.bio.symbol.Alignment; import org.biojava.bio.symbol.Alphabet; import org.biojava.bio.symbol.IllegalAlphabetException; import org.biojava.bio.symbol.IllegalSymbolException; import org.biojava.bio.symbol.IntegerAlphabet; import org.biojava.bio.symbol.SimpleAlignment; import org.biojava.bio.symbol.Symbol; import org.biojava.bio.symbol.SymbolList;
/**
* A class that uses Gibbs Sampling to generate a local alignment of an over
* represented motif.
*/
public class SimpleGibbsAligner {
private Sequence[] s; // sequence array.
private int w; //window size.
private int[] a; //starting indices.
private int iterations = 0;
private Distribution[] pattern; //the probabilistic pattern description.
private Distribution background; //the probabilistic background description.
private Random rand; //random number generator
private Alphabet alphabet; //the alphabet in which the sampler operates.
private GibbsStoppingCriteria criteria; //determines when to stop sampling.
/**
* Constructs the gibbs aligner to find a common motif in a collection
* of sequences. It is assumed that all the sequences are constructed
* from the same Alphabet. If this is not the case then calls
* to iterate will throw exceptions. This class is designed to be single use
* and is not thread safe. To use in a threaded environment each thread
* should be given its own SimpleGibbsAligner.
*
* @param windowSize the expected size of the motif
* @param it a collection of sequences in which to search for a motif.
* @param criteria an object which specifies when sampling should stop.
*/
public SimpleGibbsAligner(int windowSize,
SequenceIterator it,
GibbsStoppingCriteria criteria){
w = windowSize;
this.criteria = criteria;
rand = new Random();
//get the sequences
Vector v = new Vector();
while(it.hasNext()){
try{
v.add(it.nextSequence());
}catch(BioException e){
//cannot retreive the sequence from the iterator, not likely to happen.
e.printStackTrace();
}
}
v.trimToSize();
s = new Sequence[v.size()];
v.copyInto(s);
//intitialize the offsets
a = new int[s.length];
a = initIndices();
//set the alphabet
alphabet = s[0].getAlphabet();
}
/**
* Initialize an array of random offsets.
* @return the array of offsets
*/
private int[] initIndices(){
int[] indices = new int[s.length];
for (int i = 0; i < indices.length; i++) {
int index = rand.nextInt(s[i].length() - w-1);
// as we are making offset indices to symbollists
// they must be from 1 not 0
index++;
indices[i] = index;
}
return indices;
}
/**
* Iterates through a procedure of predictive updates and sampling until
* the stopping criteria defined in the stop() method are met.
* Once the method returns the getXXX methods can be used to
* determine the results.
*/
public void iterate(){
try {
//choose a sequence at random
int index = rand.nextInt(s.length);
do{
//calculate pattern in all but the chosen sequence
pattern = updatePattern(index, a);
//occasionaly try a phase shift
if(rand.nextDouble() < 0.1){
tryPhaseShift(index);
}
//calculate the background
background = updateBackground(index);
//sample the randomly chosen sequence to find the best start index a.
a[index] = sampleSequence(index);
//reportMatch(a[index], s[index]);
iterations++;
index = (++index)%s.length;
}while(stop() == false);
}
catch (Exception ex) {
ex.printStackTrace();
}
}
/**
* Determines when to stop iterating.
* @return true if the StoppingCriteria says to stop and false otherwise.
*/
protected boolean stop(){
return criteria.stop(this);
}
/**
* Produces a pattern to describe the motif covered by the window
* @param excludeIndex the index of the sequence to be excluded from sampling.
* @param offsets the matrix of offset positions
* @return the updated motif pattern
*/
private Distribution[] updatePattern(int excludeIndex, int[] offsets){
Distribution[] d = null;
Map label2Res = new HashMap(s.length);
for (int i = 0; i < s.length; i++) {//for each sequence
if(i == excludeIndex) continue; //except this sequence
SymbolList subSeq = s[i].subList(offsets[i],
offsets[i] +w -1);//take the subsequence
label2Res.put(new Integer(i),subSeq); //put it in the hashmap
}
Alignment al = new SimpleAlignment(label2Res);//make an alignment of subseqs
try {
d = DistributionTools.distOverAlignment(al, false,1.0);//make the pattern
}
catch (IllegalAlphabetException ex) {
ex.printStackTrace();
}
return d;
}
/**
* produces a distribution to describe the background distribution
* @param excludeIndex the index of the sequence to exclude
* @return the updated background distribution.
*/
private Distribution updateBackground(int excludeIndex){
Distribution d = null;
try {
DistributionTrainerContext dtc = new SimpleDistributionTrainerContext();
d = DistributionFactory.DEFAULT.createDistribution(alphabet);
dtc.setNullModelWeight(1.0);
dtc.registerDistribution(d);
for (int i = 0; i < s.length; i++) {//for each sequence
if(i == excludeIndex) continue; //except this sequence
for(int j = 1; j <= s[i].length(); j++){//count each base
if(j >= a[i] && j < a[i] + w-1) continue; //except these ones
dtc.addCount(d, s[i].symbolAt(j), 1.0);
}
}
dtc.train();
}
catch (Exception ex) {
ex.printStackTrace();
}
return d;
}
/**
* Attempts to prevent the pattern getting locked in a local optimum by
* shifting the pattern one step to the left or right and seeing if it is
* better than the current pattern. If the phase shift improves the model
* the pattern and offsets will be updated.
* @param excludeIndex the index of the sequence to be excluded.
*/
private void tryPhaseShift(int excludeIndex){
int[] newOffSets = new int[a.length];
System.arraycopy(a,0,newOffSets,0,a.length); // copy offsets
Distribution[] newPattern;
if (rand.nextBoolean()) {//shift left
for (int i = 0; i < newOffSets.length; i++) {
if(i == excludeIndex) continue; //skip this sequence
if(newOffSets[i] > 1) newOffSets[i]--;
}
}
else {// shift right
for (int i = 0; i < newOffSets.length; i++) {
if(i == excludeIndex) continue; //skip this sequence
if(newOffSets[i] < s[i].length() - w-2) newOffSets[i]++;
}
}
newPattern = updatePattern(excludeIndex, newOffSets);
if(getInfoContent(newPattern) > getInfoContent(pattern)){
a = newOffSets;
pattern = newPattern;
}
}
/**
* Determines a weighted distribution of offsets in the sequence to be
* sampled and randomly selects an offset from that distribution to be used
* in the next pattern update.
* @param sequenceIndex the sequence to be sampled.
* @return the selected offset
*/
private int sampleSequence(int sequenceIndex){
Distribution d = null;
try {
SymbolList seq = s[sequenceIndex];
//make an alphabet of the possible offsets
IntegerAlphabet.SubIntegerAlphabet alpha =
IntegerAlphabet.getSubAlphabet(1, seq.length()-w-1);
//make a distribution to hold the weighted probabilities of each offset.
d = DistributionFactory.DEFAULT.createDistribution(alpha);
DistributionTrainerContext dtc = new SimpleDistributionTrainerContext();
dtc.setNullModelWeight(1.0);
dtc.registerDistribution(d);
//score each subsequence
for(int i = 1; i <= seq.length()-w-1; i++){
double score = scoreSequence(seq.subList(i, i+w-1));
//add the weight to the distribution of offsets
dtc.addCount(d,alpha.getSymbol(i),score);
}
dtc.train();
}
catch (Exception ex) {
ex.printStackTrace();
}
//sample the distribution of offsets
int offset = ((IntegerAlphabet.IntegerSymbol)d.sampleSymbol()).intValue();
return offset;
}
/**
* Scores a potential motif against the pattern description and background
* distribution.
* @param sl the potential motif to score
* @return the score
*/
private double scoreSequence(SymbolList sl){
double pMotif = 1.0;
double pBackGround = 1.0;
for(int i = 0; i < sl.length(); i++){
Symbol s = sl.symbolAt(i+1); //+1 as we are indexing from zero this time
try {
pMotif *= pattern[i].getWeight(s); //probability of s at position i
pBackGround *= background.getWeight(s); //probability of s in background
}
catch (IllegalSymbolException ex) {
ex.printStackTrace();
}
}
return pMotif/pBackGround;
}
/**
* Determines the information content (in bits) of the motif inclding pseudo
* counts.
* @return the Information content.
*/
public double getInfoContent(){
return getInfoContent(pattern);
}
/**
* determines the information content (in bits) of the specified pattern
* including pseudo counts.
* @param d the pattern of the motif
* @return the information content
*/
private double getInfoContent(Distribution[] d){
double info = 0.0;
for (int i = 0; i < d.length; i++) {
info += DistributionTools.bitsOfInformation(d[i]);
}
return info;
}
/**
* Returns the current Alphabet being used.
* @return an Alphabet
*/
public Alphabet getAlphabet(){
return alphabet;
}
/**
* Get the background distribution.
* @return a Distribution of background frequencies.
*/
public Distribution getBackground() {
return background;
}
/**
* The current iteration of the sampler
* @return an int >= 0
*/
public int getIterations() {
return iterations;
}
/**
* The current pattern at this iteration of the sampler
* @return the pattern as a Distribution[].
* Effectively a weight matrix.
*/
public Distribution[] getPattern() {
return pattern;
}
/**
* Tje set of sequence offsets being used for this iteration of
* sampling
* @return an array of ints ≥ 1
*/
public int[] getOffSets(){
return a;
}
/**
* The set of Sequences being sampled
* @return a Sequence[]
*/
public Sequence[] getSequences(){
return s;
}
/**
* The size of the pattern being sampled for.
* @return an int > 0
*/
public int getWindowSize(){
return w;
}
} ```
GibbsStoppingCriteria
```java package gibbs;
import org.biojava.bio.BioException; import org.biojava.bio.dist.Distribution; import org.biojava.bio.dist.DistributionTools;
/**
* Defines the criteria under which Gibbs Sampling should stop
*/
public interface GibbsStoppingCriteria {
/**
* Uses a heuristic proceedure to determine when to stop. If the information
* content of the motif has failed to increase above its previous maximum for
* 100 iterations then the method will return true. NOTE: it is expected that
* the same SimpleGibbsSampler will be passed to the stop() method at each
* call.
*/
public static GibbsStoppingCriteria HEURISTIC = new Heuristic();
/**
* Returns true when the emission spectra of the last iteration equals that
* of this iteration. Note that this may never return if convergence is not
* reached. Thus the method has a built in stopping point of 10,000
* iterations. NOTE: it is expected that the same SimpleGibbsSampler will be
* passed to the stop() method at each call.
*/
public static GibbsStoppingCriteria CONVERGE = new Converge();
/**
* This method should return true when stopping criteria have been reached.
* @param sga the GibbsAligner that is being tested for stopping conditions
* @return true if it should stop, false otherwise.
*/
public boolean stop(SimpleGibbsAligner sga);
/**
* Implementation of GibbsStoppingCriteria
*/
class Heuristic implements GibbsStoppingCriteria{
double bestInfo = 0.0; //the level of conservation
int bestIteration = 0; //the most conserved pattern
public boolean stop(SimpleGibbsAligner sga){
double info = sga.getInfoContent();
if(info > bestInfo){
bestInfo = info;
bestIteration = sga.getIterations();
return false; //don"t stop
}else if(sga.getIterations() >= bestIteration+99){
return true;
}
return false; //don"t stop
}
}// end of Heuristic
/**
* Implementation of GibbsStoppingCriteria
*/
class Converge implements GibbsStoppingCriteria{
Distribution[] previous = null; //the last pattern
public boolean stop(SimpleGibbsAligner sga){
if(previous == null) return false; //there is no previous yet.
if(sga.getIterations() == 10000) return true; //max iterations.
try{
if (DistributionTools.areEmissionSpectraEqual(previous,sga.getPattern())){
return true; // patterns have converged.
}
else {
previous = sga.getPattern();
return false; //don"t stop
}
}catch(BioException e){
//this can"t really happen but...
e.printStackTrace();
return false;
}
}
}// end of converge
}// end of GibbsStoppingCriteria ```
SimpleGibbsAlignerDemo
```java package gibbs;
import java.io.BufferedReader; import java.io.File; import java.io.FileReader; import org.biojava.bio.seq.Sequence; import org.biojava.bio.seq.SequenceIterator; import org.biojava.bio.seq.io.SeqIOTools;
public class SimpleGibbsAlignerDemo {
/**
* Usage information
*/
public static void help(){
System.out.println(
"Usage: java SimpleGibbsAlignerDemo
public static void main(String[] args) throws Exception{
if(args.length != 4) help();
//a file of sequences sequences
File f = new File(args[0]);
//am I dealing with protein?
boolean protein = new Boolean(args[1]).booleanValue();
//the size of the motif I am looking for.
int window = Integer.parseInt(args[2]);
//the number of times to attempt a motif identification.
int trials = Integer.parseInt(args[3]);
SequenceIterator it;
for(int i = 0; i < trials; i++){
BufferedReader br = new BufferedReader(new FileReader(f));
if(protein){
it =(SequenceIterator)SeqIOTools.fileToBiojava("fasta", "protein", br);
}else{
it =(SequenceIterator)SeqIOTools.fileToBiojava("fasta", "DNA", br);
}
//make an aligner wih Heuristic stopping criteria
SimpleGibbsAligner gibbs = new SimpleGibbsAligner(window,
it, GibbsStoppingCriteria.HEURISTIC);
//start the aligner running
gibbs.iterate();
//how many iterations till convergence?
System.out.println("Converged after "+gibbs.getIterations()+" iterations");
//What is the information content of the motif?
System.out.println("Information (bits): "+gibbs.getInfoContent());
//get the sequences, offsets and window size to print out the motif
Sequence[] seqs = gibbs.getSequences();
int[] offSets = gibbs.getOffSets();
int wind = gibbs.getWindowSize();
//print out the motif
for (int j = 0; j < offSets.length; j++) {
System.out.println(seqs[j].subStr(offSets[j],offSets[j]+wind -1));
}
System.out.println();
}
}
} ```