| $g->{'name'} (edit: regulators | gene ) | EOT } $content .= "
| position | name | regulator |
| $beg-$end | $name | $reg_links |
".nazina::sequence_annotate(\$seq, \@elements)."\n"; return $content; } # -- # -- # -- motif_show # -- return motifs for a given gene sub motif_show { my ($gene_id) = @_; return "
".$dao->motifs($gene_id).""; # my $sth = $dbh->prepare("select motifs from gene where gene_id = ?"); # $sth->execute($gene_id); # my ($motifs) = $sth->fetchrow_array(); # return "
$motifs"; } # -- # -- # -- regulator_links # -- given a list of regulatory elements, return links to display # -- their motifs in a single string. sub regulator_links { my ($regulators) = @_; my @links = (); for ( sort { $a->{'abbv'} cmp $b->{'abbv'} } @{ $regulators }) { push @links, "$_->{'abbv'}"; } return (join ", \n", @links); } # -- # -- regulator_show # -- show genes grouped by regulator elemetns # -- sub regulator_show { my $content = ""; my $genes = $dao->genes(); $content .= "
| $r->{name} ($r->{abbv}) | |
| ";
# fetch elements regulated by this
$esth->execute($r_id);
while (my ($gene_id, $e_name, $g_name, $abbv) = $esth->fetchrow_array()) {
$content .= "$g_name ($abbv); $e_name \n"; } $content .= " | |
| ${ $genes }{$_}->{'name'} ( reg | gene ) | EOT } $content .= "
| position | name | regulator |
| $beg-$end | $name | @regulators |
".nazina::sequence_annotate(\$seq, \@elements)."\n"; return $content; } Appendix 2 LWF Application and Analysis Code Java Files LwfException.java package edu.dartmouth.bglab.crm.lwf; /** * Simple exception class * * @author: $Author: zburke $ * @version: $Revision: 1.2 $ * @update $Date: 2006/02/03 06:56:29 $ * */ public class LwfException extends Exception { public LwfException() { super(); } public LwfException(String s) { super(s); } public LwfException(Exception e) { super(e); } } Model.java package edu.dartmouth.bglab.crm.lwf; import java.io.BufferedReader; import java.io.FileReader; import java.io.FileNotFoundException; import java.io.IOException; /** * LWF helper class; basically a bit-bucket for model data. * * * This code is based on Perl written by Dmitri Papatsenko. The source * paper is Nazina, A. and Papatsenko, D. (2003). Statistical extraction * of Drosophila cis-regulatory modules using exhaustive assessment of * local word frequency, BMC Bioinformatics 4:65. * * @author: $Author: zburke $ * @version: $Revision: 1.7 $ * @update: $Date: 2005/05/16 02:33:15 $ * */ public class Model { /** word length */ public long length = 0; /** how many mismatches are tolerated per word */ public long mismatchCount = 0; /** size of detection window */ public long window = 0; public long windowSmall = 0; public long windowFull = 0; /** maximum number of channels to create */ public long channelCount = 0; // TODO: what does expr represent? it's calculated as the binomial // distribution of mismatches in a window over the length of a string public double expr = 0; /** string containing list of files composing this model */ public String fileList = null; /** list of files composing this model */ public String trainsets[] = null; /** */ // TODO: why is probability 0.75? public double probability = 0.75; /** maximum frequency counts for each channel */ public long[] channels = null; /** L-scores for each window in each channel */ public double[][] desmatrix = null; // nickname for this model private String name = null; // filename containing model's params and data private String filename = null; private long countmdllines = 0; private long countrows = 0; /** * default constructor does nothing interesting * */ public Model() { } /** * Creates a Model named 'name' prepared to read model parameters * from the file 'filename'. * * @param name nickname for this model * @param filename name of a file containing this models params and data * */ public Model(String name, String filename) { this.name = name; this.filename = filename; this.countmdllines = 0; this.countrows = 0; } /** * Parse model's datafile into this object. * * @throws FileNotFoundException if the model's datafile cannot be found * @throws IOException if the model's datafile cannot be read */ public void fileParse() throws FileNotFoundException, IOException { BufferedReader br = null; String line = null; int j = 0; br = new BufferedReader(new FileReader(this.filename)); while ((line = br.readLine()) != null) { if (line.trim().equals("")) continue; String[] words = line.split(" "); // word-length if (line.matches("Word:.*")) this.length = (new Long(words[1])).longValue(); // mismatch count if (line.matches("Mism:.*")) this.mismatchCount = (new Long(words[1])).longValue(); // sliding window size if (line.matches("Win:.*")) { this.window = new Long(words[1]).longValue() / 2; this.windowSmall = this.window; this.windowFull = this.windowSmall * 2; } // frequency limits for each channel if (line.matches("Limits:.*")) { this.channels = new long[words.length - 1]; for (int i = 1; i < words.length; i++) this.channels[i-1] = (new Long(words[i])).longValue(); this.channelCount = this.channels.length - 1; } // L-scores for a channel define strength of a channel as an // predictor of membership in a training set if ( line.matches("DMatrix:.*")) { if (this.desmatrix == null) this.desmatrix = new double[(int)this.channelCount][]; this.desmatrix[j] = new double[words.length - 1]; for (int i = 1; i < words.length; i++) this.desmatrix[j][i-1] = (new Double(words[i])).doubleValue(); j++; } this.countmdllines++; } this.countrows = j; } } Scan.java package edu.dartmouth.bglab.crm.lwf; import edu.dartmouth.bglab.crm.lwf.scan.ProgressBar; import edu.dartmouth.bglab.crm.lwf.scan.ScanSequenceParams; import edu.dartmouth.bglab.crm.lwf.train.TrainingSet; import edu.dartmouth.bglab.crm.lwf.train.TrainSequenceParams; import java.io.BufferedReader; import java.io.BufferedWriter; import java.io.FileNotFoundException; import java.io.FileReader; import java.io.FileWriter; import java.io.IOException; import java.io.PrintWriter; import java.io.File; import java.io.FileInputStream; import java.io.BufferedInputStream; import java.text.DecimalFormat; import java.text.NumberFormat; import java.util.Arrays; import java.util.List; import java.util.ArrayList; import java.util.Iterator; import java.util.Hashtable; import java.util.Date; import java.util.Properties; /** * Scans test sequence using statistical models generated by Lwf-train and * generates a recognition profile * * Input: * - file containing input sequences in fasta format, *.fa (test set) * - at least two *.sts files generated by LWF-train and containing statistics * for the positive and the negative training sets. Both *.sts files must * contain the same input parameters for the scan (see *.sts file headers). * * Options: * - equalize channels: low values will supress difference between channels * - profile cutoff * - extraction option (fraction of sequence length) * - peak width cutoff * - smoothing window * - GNUPlot option - can display profile for sequences < 1MB * * Output: * - a directory containing files: * - *.mdl - model for each pair of the input *.sts files * - *.dat - recognition profiles for each sequence from the input *.fa file. * - *.xtr.dat - coordinates of the extracted sequences * - *.xtr.fa - extracted sequences * * This code is based on Perl written by Dmitri Papatsenko. The source * paper is Nazina, A. and Papatsenko, D. (2003). Statistical extraction * of Drosophila cis-regulatory modules using exhaustive assessment of * local word frequency, BMC Bioinformatics 4:65. * * @author: $Author: zburke $ * @version: $Revision: 1.14 $ * @update: $Date: 2006/02/03 06:56:29 $ * */ public class Scan { public static void main(String[] args) { try { Scan s = null; // use properties file for runtime params if (args != null && args.length > 0) s = new Scan(args[0]); else s = new Scan(); // get runtime parameters including the sequence to parse // and the models to scan against, and read the raw sequence // data. s.train(); // create an LWF profile of the sequence and calculate // the feature score for the windows starting at each position // in the sequence. results are saved in a temp-scan file. s.scan(); // parse temp-scan files to find windows with feature scores // beyond the cutoff representing stronger resemblance to the // positive, rather than the negative, training set, i.e. // sequence substrings classified as CRMs. s.classify(); System.out.print("\n\n"); } catch (FileNotFoundException e) { e.printStackTrace(); } catch (IOException e) { e.printStackTrace(); } catch (LwfException e) { System.err.println("ERROR: "+e.getMessage()); e.printStackTrace(System.err); } } // positive/negative training sequences used to parse this sequence private TrainingSet ts = null; // run-time params private ScanSequenceParams params = null; // lists putative CRM sequence subscripts private PrintWriter xtr = null; // lists putative CRM sequences in FASTA format private PrintWriter fa = null; // model properties file private Properties props = null; // for pretty printing private NumberFormat nf = null; // max length of sequence to use for GNUPlot scanning results private int GNUPLOT_MAX = 1000000; /** * Constructs a new CRM scanner and initializes the output files * @throws IOException */ public Scan() throws IOException { init(); } /** * Constructs a new CRM scanner and reads sequence and scanning * data from a properties file. * * @param pfile properties file to read */ public Scan(String pfile) { try { init(); Properties runtimeProps = new Properties(); BufferedInputStream b = null; b = new BufferedInputStream(new FileInputStream(pfile)); if (b != null) runtimeProps.load(b); this.props = runtimeProps; } // error reading properties file; set it to null and continue // as if we never had a file to read. catch (IOException e) { this.props = null; } } /** * Initialize instance variables. * */ private void init() { this.ts = new TrainingSet(); this.params = new ScanSequenceParams(); this.nf = NumberFormat.getInstance(); if (this.nf instanceof DecimalFormat) ((DecimalFormat) this.nf).applyPattern("#.##############################"); } /** * Get training params and sequence filenames, then create models * from the training data. * * @throws FileNotFoundException if any input file cannot be found * @throws IOException if any input file cannot be read * @throws LwfException if there are problems parsing a training set * (e.g. parse params for positive and negative sets don't match) */ private void train() throws FileNotFoundException, IOException, LwfException { // if we have a properties file, get sequence and model filenames // from it; otherwise, collect params from user if (this.props != null) { this.params.trainSet( this.props.getProperty("sequence"), this.props.getProperty("train"), (new Integer(this.props.getProperty("equalize"))).intValue()); } else { this.params.trainGet(); } System.out.print("\n\nCREATING MODELS...\n\n"); trainParse(); // get scanning params String propsScan = this.props.getProperty("scan"); if (this.props != null && propsScan != null && propsScan.equals("defaults")) { this.params.scanSet( this.params.profileCutoffSuggest(), this.params.peakWidthSuggest(), this.params.smoothSuggest(), this.params.plotSuggest()); } else if (this.props != null) { this.params.scanSet( (new Double(this.props.getProperty("profile-cutoff"))).doubleValue(), (new Integer(this.props.getProperty("peak-width"))).intValue(), (new Integer(this.props.getProperty("smooth"))).intValue(), (new Boolean(this.props.getProperty("plot"))).booleanValue()); } else { this.params.scanGet(); } // read sequence file this.ts.sequenceFileParse(this.params.seqfname, this.params); } /** * Iterate through pairs of positive and negative TrainingSet files, * parsing each into a model which is written to a model profile * file (*.mdl). * * @throws FileNotFoundException if a training file cannot be found * @throws IOException if a training file cannot be read * @throws LwfException if a training set parse params don't match */ private void trainParse() throws FileNotFoundException, IOException, LwfException { Iterator i = this.params.stsList.iterator(); while (i.hasNext()) { // get training set files Hashtable h = (Hashtable) i.next(); String stsPos = (String) h.get("pos"); String stsNeg = (String) h.get("neg"); // create model file String[] nameList = stsPos.split("\\."); String mdlFile = this.params.scanoutdir + "/" + nameList[0] + ".mdl"; PrintWriter mpw = new PrintWriter(new BufferedWriter(new FileWriter(mdlFile))); // parse training set into a model modelCreate(mpw, stsPos, stsNeg); this.params.modelNames.add(nameList[0]); mpw.flush(); mpw.close(); mpw = null; } } /** * Parse a positive and negative training-set pair into a * sequence model. Store the model data in a model profile * file (*.mdl). * * @param mpw where to write model data * @param stsPos filename of positive training set * @param stsNeg filename of negative training set * * @throws FileNotFoundException if a training set file cannot be found * @throws IOException if a training set file cannot be read * @throws LwfException if a training set parse params don't match */ private void modelCreate(PrintWriter mpw, String stsPos, String stsNeg) throws FileNotFoundException, IOException, LwfException { // stores all statistics from both training sets ArrayList statistics = new ArrayList(); // read training sets and make sure parse params match TrainSequenceParams tPos, tNeg = null; tPos = trainParsePos(statistics, stsPos, mpw); tNeg = trainParseNeg(statistics, stsNeg, mpw); if (tPos.channelCount != tNeg.channelCount || tPos.length != tNeg.length) throw new LwfException(stsPos+" and "+stsNeg+" matrix dimensions don't match; model excluded\n"); // weight cummulative E-values for feature scores within each channel double[][] matrix = new double[(int)tPos.channelCount][]; for (int j = 0; j < tPos.channelCount; j++) { double[] cPos = (double[]) statistics.get(j); // TODO: why grab a different channel for the negative set???? double[] cNeg = (double[]) statistics.get(j + (int)tPos.channelCount ); matrix[j] = trainWeightEValues(cPos, cNeg); } // write training model to an output file trainModelDistributionWrite(mpw, matrix); } /** * Read data from the positive training set and write it to the * model file. * * //TODO: rewrite with fieldnames instead of line numbers * * @param list collection of E-value distributions for each channel * in the training set. * @param filename trainingset filename to read * @param mpw model file to write to * * @return a TrainSequenceParams object with its parse-parameters * filled in * * @throws FileNotFoundException if the trainingset file cannot be found * @throws IOException if the trainingset file cannot be read * */ private TrainSequenceParams trainParsePos(ArrayList list, String filename, PrintWriter mpw) throws FileNotFoundException, IOException { TrainSequenceParams tsp = new TrainSequenceParams(); int i = 0; BufferedReader br = null; String line = null; br = new BufferedReader(new FileReader(filename)); while ((line = br.readLine()) != null) { line = line.trim(); // collect Word, Mism, Win, Exp, Chann, Limits params if (i > 0 && i < 7) mpw.println(line); // add suppression (how much to supress extreme L-scores) if (i == 7) { mpw.println("Sigma: "+this.params.supression); mpw.println(""); } String[] words = line.split(" "); // window size; "Win: " if (i == 3) { tsp.window = (new Long(words[1])).longValue(); this.params.mdlwins.add(new Long(tsp.window)); } // channel count; "Chann: " if (i == 5) tsp.channelCount = (new Long(words[1])).longValue(); // derived detection window size if (i == 10) tsp.length = words.length; if (i > 9) list.add(trainModelTrainWrite(mpw, words, "TrainPos")); i++; } br.close(); mpw.println(""); mpw.flush(); return tsp; } /** * Read data from the negative training set and write it to * the model file. * * // TODO: rewrite with fieldnames instead of line numbers * * @param list collection of E-value distributions for each channel * in the training set. * @param filename trainingset filename to read * @param mpw model file PrintWriter * * @return a TrainSequenceParams object with its parse-parameters * filled in * * @throws FileNotFoundException if the trainingset file cannot be found * @throws IOException if the trainingset file cannot be read */ private TrainSequenceParams trainParseNeg(ArrayList list, String filename, PrintWriter mpw) throws FileNotFoundException, IOException { int i = 0; TrainSequenceParams tsp = new TrainSequenceParams(); BufferedReader br = null; String line = null; br = new BufferedReader(new FileReader(filename)); while ((line = br.readLine()) != null) { line = line.trim(); String[] words = line.split(" "); // channel count; Chann: if (line.matches("Chann:.*")) tsp.channelCount = (new Long(words[1])).longValue(); // derived detection window size if (i == 10) tsp.length = words.length; // L-scores for each word in a window if (i > 9) list.add(trainModelTrainWrite(mpw, words, "TrainNeg")); i++; } br.close(); mpw.println(""); mpw.flush(); return tsp; } /** * Write feature scores for a channel in a training set to the * model file. * * @param mpw model file PrintWriter * @param evalues String representations of feature scores for a * channel, as read from a training file. * @param header field header indicating data-type on this line * * @return coverted feature scores (String -> double) for this channel */ private double[] trainModelTrainWrite(PrintWriter mpw, String[] evalues, String header) { double[] numbers = new double[evalues.length]; mpw.print(header +": "); for (int j = 0; j < evalues.length; j++) { //System.err.println("value is: "+evalues[j]); mpw.print(evalues[j] + " "); numbers[j] = ((Double) new Double(evalues[j])).doubleValue(); } mpw.println(""); return numbers; } /** * Weight the distribution (E) of feature scores (S) within a channel * (j) away from the value S-0 at which membership in either the * positive or negative training set is equally likely. * * From Nazina and Papatsenko regardging equations 2 and 3 and error * correction: * * Given functional classes omega-1 and omega-2 representing positive * and negative training sets, and given a feature score S for a test * sequence, the classification problem is solved for each frequency * channel (j) independently using a log-likelihood ratio: * * [Equation 2] * * L-j == log (E [omega-1 | S-j] / E [omega-2 | S-j] ) * * * where E is the distribution of S in each channel j. L-j is zero * where E [omega-1 | S-j] == E [omega-2 | S-j], i.e. where membership * in either class is equally likely. As noted in the Methods section, * this value is much closer to the expectation for the negative * training set than for the positive one. To account for this error, * the cummulative E values are weighted away from S-0 as follows: * * [Equation 3] * * L == log (E [omega-1 | S] / E [omega-2 | S] ) * + * log ((sum(E [omega-1 | S]) * (1 - sum(E [omega-1 | S]))) * / * (sum(E [omega-2 | S]) * (1 - sum(E [omega-2 | S]))) * ) * * * @param cPos E-values for one channel the positive training set * @param cNeg E-values for one channel the negative training set * * @return weighted E-values for this channel * */ private double[] trainWeightEValues(double[] cPos, double[] cNeg) { // weighted E-values for this channel double[] row = new double[(int)cPos.length]; // TODO: figure out why the weighted E-values have this skew attached double skew = 0.001; for (int i = 0; i < cPos.length; i++) { double sumPos = Util.sumArray(cPos, 0, i); double sumNeg = Util.sumArray(cNeg, 0, i); // log(a/b) + log ((sumA (1 - sumA + a)) / sumB (1 - sumB + b)) double value = Math.log( (cPos[i] + skew) / (cNeg[i] + skew) ) + Math.log( ((sumPos * (1 - sumPos + cPos[i] )) + skew) / ((sumNeg * (1 - sumNeg + cNeg[i] )) + skew) ); row[i] = value * Util.normalDistribution(value, 0, this.params.supression); } return row; } /** * Write weighted E-values for a channel in combined training set to * the model file. * * @param mpw model file PrintWriter * @param evalues weighted E-values for a channel as determined by * processing the positive and negative training sets. * * @see trainWeightEValues */ private void trainModelDistributionWrite(PrintWriter mpw, double[][] evalues) { for (int i = 0; i < evalues.length; i++) { mpw.print("DMatrix: "); for (int j = 0; j < evalues[i].length; j++) mpw.print(evalues[i][j] + " "); mpw.println(""); } } /** * Scan the sequences using the models created from the training sets. * Recognition data temporarily is stored in a file for later parsing. * Each sequence is parsed with all available models. * * @throws FileNotFoundException * @throws IOException */ private void scan() throws FileNotFoundException, IOException { System.out.print("\n\nPARSING MODELS...\n\n"); Iterator i = this.params.modelNames.iterator(); while (i.hasNext()) modelParse((String) i.next()); } /** * Parse each sequence with the named model. * * @param modelName name of model to use for parsing * * @throws FileNotFoundException if there is trouble finding a model file * @throws IOException if there is trouble reading a model file */ private void modelParse(String modelName) throws FileNotFoundException, IOException { // read model data from its source file String modelFile = this.params.scanoutdir + "/" + modelName + ".mdl"; Model model = new Model(modelName, modelFile); model.fileParse(); // process each sequence with this model Iterator iter = ts.sequences.iterator(); while (iter.hasNext()) { Sequence s = (Sequence) iter.next(); // find beginning of last window in the sequence s.maxPosition = s.length - model.windowSmall - model.length + 1; String outfile = this.params.scanoutdir + "/" + s.filename + "~" + modelName +".tmp"; PrintWriter out = new PrintWriter(new BufferedWriter(new FileWriter(outfile))); Date timeStart = new Date(); ProgressBar progress = new ProgressBar(s.data, s.length, model.window, modelName); // build channel distribution matrix for a sequence. // the index of the start of the final sliding window is returned // so word distributions for all positions beyond that index can // be zeroed out. long pos = sequenceParse(model, s, out, progress); long print0 = pos - model.window + 1; while (print0++ < pos + model.window + 3) out.println("0"); out.flush(); out.close(); out = null; Date timeEnd = new Date(); System.out.println(""); System.out.println(runtimeDisplay(timeStart, timeEnd, s)); } } /** * Create a channel distribution matrix for the given sequence * with the given model. The results are written to a file. * * @param model model to use for parsing * @param s sequence to parse * @param out PrintWriter for temp sequence recognition file * @param progress progress meter * * @return index of the start of the final sliding window within the * sequence, beyond which no frequency-distribution calculations * were performed. */ private long sequenceParse(Model model, Sequence s, PrintWriter out, ProgressBar progress) { List channels, channelWordCount; channels = new ArrayList(); channelWordCount = null; long pos = model.window; int i = 0; // TODO: why start at model.window instead of at 0? while (pos < s.maxPosition) { // detect channels for words in this window // i.e. list of word-counts, per channel channelWordCount = Util.channelDetect(s, pos, channels, model); // feature score of a channel representing the number of // words in a window with a frequency in a given range double featureScore = 0; if (channelWordCount != null && channelWordCount.size() > 1) { // channel number int j = 0; // loop over channels Iterator iter = channelWordCount.iterator(); while (iter.hasNext()) { long wordsInChannel = ((Long)iter.next()).longValue(); featureScore += model.desmatrix[j][(int)wordsInChannel]; j++; } } out.println(featureScore); // if (i == progress.progress) // progress.update(i); i++; pos++; } return pos; } /** * Classify windows as resembling the positive or negative training * set. Positive windows, i.e. putative CRMs, are extracted and * written to output files. * * @throws FileNotFoundException if the scan file cannot be found * @throws IOException if the scan file cannot be read */ private void classify() throws FileNotFoundException, IOException { System.out.print("\n\nWRITING SEQUENCE PROFILE DATA...\n\n"); // gathers sequence data for a GNUPlot file ArrayList plot = new ArrayList(); // prepare output files: // xtr contains substring indexes of putative CRMs in the input file // fa contains FASTA sequences of the putative CRMs in the input file String xtrFile = this.params.scanoutdir+"/" + this.params.seqfname + ".xtr.dat"; this.xtr = new PrintWriter(new BufferedWriter(new FileWriter(xtrFile))); String faFile = this.params.scanoutdir+"/" + this.params.seqfname + ".xtr.fa"; this.fa = new PrintWriter(new BufferedWriter(new FileWriter(faFile))); // parse sequence-profile temp files created during scanning to find // sequences with sufficient similarity to the positive training set. // write data to the xtr and fa files Iterator i = ts.sequences.iterator(); while (i.hasNext()) sequenceClassifyConfig((Sequence) i.next(), plot); xtr.flush(); xtr.close(); fa.flush(); fa.close(); // write GNUPlot data file for short-ish sequences if (ts.length < this.GNUPLOT_MAX) plot(plot); System.out.print("\n\n"); } /** * Determine positive-classification score threshold, then classify * the windows. Data for positive windows is saved in a file. * * A window's resemblance to the positive or negative training set * is calculated as a simple log-likelihood ratio of predicted * membership in the positive set over predicted membership in the * negative set. Windows with scores beyond a certain threshold are * classified as positive. * * Before this score is calculated, the S-scores (used to calculate * the L-score) are smoothed to remove extreme values that would * otherwise dominate a channel's recognition profile. * * @param s the Sequence containing windows to classify * @param plotData bit bucket for plotting sequence profile with GNUPlot * * @see classifyThreshold */ private void sequenceClassifyConfig(Sequence s, List plotData) throws FileNotFoundException, IOException { String filename = this.params.scanoutdir + "/" + s.filename + ".dat"; PrintWriter dat = new PrintWriter(new BufferedWriter(new FileWriter(filename))); double[] profileTmp = spTmpDataParse(s); double[] profile = Util.smooth(this.params.inputSmooth, profileTmp); double threshold = classifyThreshold(profile); sequenceClassify(threshold, profile, s, dat); if (s.length < GNUPLOT_MAX) plotData.add(profile); dat.flush(); dat.close(); dat = null; } /** * Retrieve Sequence's scan data (feature scores for each window in * the Sequence) from the tmp file. * * //TODO: rewrite this with file-filters * * @param s Sequence for which data will be retrieved * * @return scan profile data for the input sequences * * @throws FileNotFoundException if the Sequence's scan file cannot be found * @throws IOException if the Sequence's scan file cannot be read */ private double[] spTmpDataParse(Sequence s) throws FileNotFoundException, IOException { double[] profile = new double[(int)s.length]; File dir = new File(this.params.scanoutdir); String[] files = dir.list(); for (int i = 0; i < files.length; i++) { if (files[i].matches(s.data+"~.*")) { BufferedReader br = null; String line = null; br = new BufferedReader(new FileReader(this.params.scanoutdir + "/"+files[i])); int j = 0; while ((line = br.readLine()) != null) { line = line.trim(); if (! line.equals("")) profile[j] = (new Double(line)).doubleValue() + profile[j]; j++; } br.close(); File tmp = new File(this.params.scanoutdir + "/"+files[i]); tmp.delete(); tmp = null; } } return profile; } /** * Calculate the L-score threshold that serves as the border between * positive or negative classification as a CRM. * * A window's resemblance to the positive or negative training set * is calculated as a simple log-likelihood ratio of predicted * membership in the positive set over predicted membership in the * negative set. Windows with scores beyond a certain threshold are * classified as positive. * * @param bins * @param profiles feature scores for each window * * @return threshold distinguishing positive and negative CRM status */ private double classifyThreshold(double[] profile) { // number of bins to distribute profile[] items (feature scores) into int bins = 1000; Hashtable hash = distributeFeatureScores(bins, profile); double[] distribution = (double[]) hash.get("distribution"); double max = ((Double) hash.get("max")).doubleValue(); double min = ((Double) hash.get("min")).doubleValue(); int count = distribution.length - 1; for (double i = 1; i > 1 - this.params.profileCutoff; count--) i = distribution[count]; double cutoff = min + ( ( count * ( max - min ) ) / bins ); return cutoff; } /** * Distribute feature scores for each channel across a number * of bins. * * @param bins number of bins into which scores are sorted * @param scores feature scores for each window * * @return a hash containing the max and min feature scores and * an array of distributed feature scores */ private Hashtable distributeFeatureScores(int bins, double[] scores) { // copy profile so original isn't modified; sort copy descending double[] profileSort = new double[scores.length]; System.arraycopy(scores, 0, profileSort, 0, scores.length); Util.reverseSort(profileSort); double min = profileSort[profileSort.length - 1]; double max = profileSort[0]; // feature score histogram double[] distribution = distribute(profileSort, min, max, bins); double sumArray = Util.sumArray(distribution, 0, distribution.length); double sum = 0.0; for (int i = 0; i < distribution.length; i++) { distribution[i] /= sumArray; sum += distribution[i]; distribution[i] = sum; } Hashtable hash = new Hashtable(); hash.put("max", new Double(max)); hash.put("min", new Double(min)); hash.put("distribution", distribution); return hash; } /** * Distribute feature scores for each channel across a number * of bins, i.e. create a histogram of feature scores. * * @param scores all feature scores * @param min smallest feature scores * @param max largest feature score * @param bins number of bins into which scores are sorted * * @return a feature score histogram */ private double[] distribute(double[] scores, double min, double max, int bins) { double[] distribution = new double[(int)bins + 1]; Arrays.fill(distribution, 0); double binSize = (max - min) / bins; if (binSize <= 0) return distribution; for (int i = 0; i < scores.length; i++) { // TODO: what is the 0.000001 adjustment for? int current = (int) ((scores[i] - min - 0.000001) / binSize); if (current <= bins) distribution[current]++; } return distribution; } /** * Find sequence windows with L-scores above the threshold representing * classification as a putative CRM. * * Classification data for all windows is written to the dat file. * For positively-classified sequences, sequence data and substring * indexes of the window in the source sequence are written to the * fa and xtr files. * * @param threshold L-score boundary between classification as positive * or negative * @param scores feature scores for each window in this sequence * @param s Sequence to classify * @param dat PrintWriter for recognition profile of this Sequence */ private void sequenceClassify(double threshold, double[] scores, Sequence s, PrintWriter dat) { long pos = 0, posStart = 0; boolean inCRM = false; for (int i = 0; i < scores.length; i++) { pos++; dat.println(scores[i]); // found the end of a CRM; save it if (inCRM && scores[i] < threshold) { inCRM = false; spProfilePrint(posStart, pos, s); } // found the beginning of a CRM; save its start index if ((! inCRM) && scores[i] > threshold) { inCRM = true; posStart = pos; } } // print final data if we ended in a CRM if (inCRM) { spProfilePrint(posStart, pos, s); } } /** * write sequence data to output files * * @param end start of this window * @param begin end of this window * @param s sequence to write */ private void spProfilePrint(long begin, long end, Sequence s) { // CRMs must be wider than params.peakWidth if (end - begin > this.params.peakWidth) { this.xtr.println(s.data +" "+ begin + "-"+ end); String fasta = Util.fasta(s.sequence.substring(((int)begin) - 1,(int) end)); this.fa.println(">"+ s.data +" "+ begin +"-"+ end +"\n"+ fasta); } } /** * Write GNUPlot config file for displaying a sequence profile. * * @param list */ private void plot(List list) throws IOException { String plotFile = this.params.scanoutdir + "/" +this.params.seqfname + ".plt"; PrintWriter out = new PrintWriter(new BufferedWriter(new FileWriter(plotFile))); double[] profile = null; Iterator i = list.iterator(); while (i.hasNext()) { profile = (double[]) i.next(); for (int j = 0; j < profile.length; j++) out.print(profile[j] + " "); out.println(); } out.flush(); out.close(); out = new PrintWriter(new BufferedWriter(new FileWriter("gnuplot.ini"))); out.println("plot '"+ plotFile +"' with lines"); out.flush(); out.close(); if (this.params.plot) { // system("wgnuplot") if $p->{'plotit'}; } } /** * Calculate running time and processing speed (KB/second). * * @param start when processing started * @param end when processing ended * @param s the Sequence that was processed. */ private String runtimeDisplay(Date start, Date end, Sequence s) { // running time in seconds long time = (end.getTime() - start.getTime()) / 1000; int speed = (int) ((s.length / (10 * time)) / 100); int hour = (int) (time / (60 * 60)); int min = (int) ((time / 60) % 60); int sec = (int) (time % 60); return ("runtime: " +(hour > 9 ? ""+hour : "0" + hour)+":" +(min > 9 ? ""+min : "0" + min)+":" +(sec > 9 ? ""+sec : "0" + sec)+" ("+speed+" KB/s)"); } } Sequence.java package edu.dartmouth.bglab.crm.lwf; import java.io.PrintStream; import edu.dartmouth.bglab.crm.lwf.train.TrainSequenceParams; /** * LWF helper class; basically a bit-bucket for sequence-related data. * * This code is based on Perl written by Dmitri Papatsenko. The source * paper is Nazina, A. and Papatsenko, D. (2003). Statistical extraction * of Drosophila cis-regulatory modules using exhaustive assessment of * local word frequency, BMC Bioinformatics 4:65. * * * @author: $Author: zburke $ * @version: $Revision: 1.3 $ * @update: $Date: 2006/02/07 23:10:23 $ * */ public class Sequence { /** length of sequence in chars */ public long length = 0; /** file header data */ public String data = null; /** */ public String input = null; /** the sequence characters */ public String sequence = null; /** */ public long maxPosition = 0; /** source filename */ public String filename = null; /** * simple constructor requires the filename this sequence was * created from * */ public Sequence(String filename) { this.filename = filename; this.input = new String(); this.sequence = new String(); } /** * simple constructor sets some params * * @param data * @param length * @param sequence * @param input */ public Sequence(String data, long length, String sequence, String input) { this.data = data; this.length = length; this.sequence = sequence; this.input = input; } /** * print information about this seq to the specified stream * @param sequenceName * @param p * @param out */ public void headerPrint(String sequenceName, TrainSequenceParams p, PrintStream out) { out.println("Data: " + sequenceName); out.println("Word: " + p.length); out.println("Mism: " + p.mismatchCount); out.println("Win: " + p.window); out.println("Exp: " + p.expr); out.println("Chann: " + p.channelCount); out.println("Limits: " + p.channelsString()); out.println(""); out.println("Statistics: "); out.println(""); } } Train.java package edu.dartmouth.bglab.crm.lwf; import java.io.BufferedInputStream; import java.io.FileInputStream; import java.io.FileNotFoundException; import java.io.IOException; import java.io.BufferedWriter; import java.io.PrintWriter; import java.io.FileWriter; import java.util.ArrayList; import java.util.List; import java.util.Iterator; import java.util.Date; import java.util.Properties; import java.text.NumberFormat; import java.text.DecimalFormat; import edu.dartmouth.bglab.crm.lwf.train.TrainProgressBar; import edu.dartmouth.bglab.crm.lwf.train.TrainSequenceParams; import edu.dartmouth.bglab.crm.lwf.train.TrainingSet; /** * Generates statistical model for a training set based on word * frequency distributions. * * input * * files containing input sequences in fasta format * * parameters * * word length * * string divergence (number of mismatches) * * sliding window size * * number of frequency channels * * output * * .sts files containing corresponding word frequency statistics * for each training set * * * This code is based on Perl written by Dmitri Papatsenko. The source * paper is Nazina, A. and Papatsenko, D. (2003). Statistical extraction * of Drosophila cis-regulatory modules using exhaustive assessment of * local word frequency, BMC Bioinformatics 4:65. * * * @author: $Author: zburke $ * @version: $Revision: 1.9 $ * @update: $Date: 2005/05/16 01:20:05 $ * */ public class Train { /** * collect training parameters, then parse the input file into * a training set file. * * @param args */ public static void main(String[] args) { try { Train t = null; if (args != null && args.length > 0) t = new Train(args[0]); else t = new Train(); // get and display calculated params // user properties if availablet; otherwise prompt user if (t.props != null) { t.params.set( (new Integer(t.props.getProperty("word-length"))).intValue(), (new Integer(t.props.getProperty("mismatch-count"))).intValue(), (new Integer(t.props.getProperty("window-length"))).intValue(), (new Integer(t.props.getProperty("channel-count"))).intValue(), t.props.getProperty("files") ); } else { t.params.get(); } t.params.show(); Date startTime = new Date(); // parse the training files for (int i = 0; i < t.params.trainsets.length; i++) t.parseFile(t.params.trainsets[i]); Date endTime = new Date(); long time = endTime.getTime() - startTime.getTime(); System.out.println("runtime: "+(time/1000)+" seconds\n\n"); } catch (FileNotFoundException e) { e.printStackTrace(); } catch (IOException e) { e.printStackTrace(); } catch (LwfException e) { System.err.println(e.getMessage()); } } // runtime params public TrainSequenceParams params = null; // model properties file public Properties props = null; // for pretty printing private NumberFormat nf = null; /** * simple constructor initializes instance variables * */ public Train() { init(); } /** * Constructs a new Trainer with runtime parameters set in a * properties file. * * @param pfile name of the properties file * * @throws IOException if the properties file cannot be read */ public Train(String pfile) throws IOException { try { init(); Properties runtimeProps = new Properties(); BufferedInputStream b = null; b = new BufferedInputStream(new FileInputStream(pfile)); if (b != null) runtimeProps.load(b); else System.err.println("ERROR: could not read "+pfile); this.props = runtimeProps; } // error reading properties file; set it to null and continue // as if we never had a file to read. catch (IOException e) { System.err.println("couldn't read "+pfile); this.props = null; } } /** * Initialize instance variables. * */ private void init() { this.params = new TrainSequenceParams(); this.nf = NumberFormat.getInstance(); if (this.nf instanceof DecimalFormat) ((DecimalFormat) this.nf).applyPattern("#.##############################"); } /** * parse a FASTA file into a training set file * * @param filename FASTA sequence file to read * * @throws FileNotFoundException if the sequence file cannot be found * @throws IOException if the sequence file cannot be read */ private void parseFile(String filename) throws FileNotFoundException, IOException { // distribution matrix of words in channels long[][] distribution = Util.matrix(this.params.windowFull, this.params.channelCount, 0); int shiftRows = 0; // output filename + header String outfile = this.params.length +"_"+this.params.mismatchCount+"_"+this.params.window+"_"+this.params.channelCount+"_"+filename+".sts"; // read sequence file, then prepare outfile TrainingSet ts = new TrainingSet(filename, this.params); PrintWriter out = new PrintWriter(new BufferedWriter(new FileWriter(outfile))); headerPrint(filename, out); TrainProgressBar progress = new TrainProgressBar(filename, ts.length); // parse the sequences in the input file int count = 0; java.util.Iterator i = ts.sequences.iterator(); while (i.hasNext()) count = parseSequence(distribution, shiftRows, (Sequence) i.next(), progress, count); shiftRows++; distributionMatrixWrite(distribution, out); out.flush(); out.close(); System.out.println(""); System.out.println(""); } /** * build distribution matrix for a sequence * * @param matrix * @param shiftrows * @param s sequence to parse * @param progress progress meeter * */ private int parseSequence(long[][] matrix, long shiftrows, Sequence s, TrainProgressBar progress, int i) { // TODO: why initialize to windowSmall instead of 0? long pos = this.params.windowSmall; List detectChann = new ArrayList(); List channelNSum = new ArrayList(); while (pos < s.maxPosition) { // detect channels for words in this window channelNSum = Util.channelDetect(s, pos, detectChann, this.params); if (channelNSum != null && channelNSum.size() > 1) { // System.out.println("channelNSum.size: "+channelNSum.size()); long countChanSum = 0; Iterator iter = channelNSum.iterator(); while (iter.hasNext()) { int test1 = (int) (this.params.channelCount * shiftrows + countChanSum); int test2 = ((Long) iter.next()).intValue(); // TODO: figure out why test2 is ever successful if (test1 >= matrix.length || test2 >= matrix[test1].length) { System.out.println("DANGER: "+test1+" >= " +matrix.length + " || "+test2 + " >= "+ matrix[test1].length); continue; } // matrix[ new Long(this.params.channelCount * shiftrows + countChanSum).intValue() ][ ((Long) iter.next()).intValue() ]++; matrix[ new Long(this.params.channelCount * shiftrows + countChanSum).intValue() ][ test2 ]++; countChanSum++; } } // if (i == progress.progress) // progress.update(i); i++; pos++; } return i; } /** * write the distribution matrix for each channel to the output file * * @param channels 2D array of windows in each channel * @param out output file */ private void distributionMatrixWrite(long[][] channels, PrintWriter out) { int distrmrows = 0; for (int i = 0; i < channels.length; i++) { long[] current = channels[i]; long unique = current[0]; int count0 = 0; for (int j = 1; j < current.length; j++) { if (0 == current[j]) count0++; } // smooth values and shift unique value back onto list double[] smoothMatrix = smoothed(current, count0, unique); double sumSmoothMatrix = Util.sum(smoothMatrix); for (int j = 0; j < smoothMatrix.length; j++) smoothMatrix[j] /= sumSmoothMatrix; // print it if (distrmrows++ == this.params.channelCount) out.println(); for (int j = 0; j < smoothMatrix.length; j++) out.print(this.nf.format(smoothMatrix[j]) + " "); out.println(); } } /** * Smooth values on the list, and return a list with the unique * value shifted back onto the start of th list. The original list * is not modfied. * * @param list list of values to smooth * @param count0 * @param unique list value which is exempt from smoothing * * @return smoothed list */ private double[] smoothed(long[] list, int count0, long unique) { long[] currentShift = new long[list.length - 1]; System.arraycopy(list, 1, currentShift, 0, currentShift.length); int smoothRate = StrictMath.round((this.params.windowFull - count0) / 10); double[] smoothMatrixTemp = Util.smooth(smoothRate, currentShift); double[] smoothMatrix = new double[smoothMatrixTemp.length + 1]; smoothMatrix[0] = unique; System.arraycopy(smoothMatrixTemp, 0, smoothMatrix, 1, smoothMatrixTemp.length); return smoothMatrix; } /** * print output file header * * @param sequenceName * @param out */ public void headerPrint(String sequenceName, PrintWriter out) { out.print("Data: " + sequenceName + "\n"); out.print("Word: " + this.params.length + "\n"); out.print("Mism: " + this.params.mismatchCount + "\n"); out.print("Win: " + this.params.window + "\n"); out.print("Exp: " + this.params.expr + "\n"); out.print("Chann: " + this.params.channelCount + "\n"); out.print("Limits: " + this.params.channelsString() +"\n\n"); out.print("Statistics:\n\n"); } } Util.java package edu.dartmouth.bglab.crm.lwf; import java.util.ArrayList; import java.util.List; import java.util.Iterator; import java.util.Collection; import java.util.Arrays; import edu.dartmouth.bglab.crm.lwf.Word; /** * helper methods for train.pl, including input parsing and mathematical * functions * * * This code is based on Perl written by Dmitri Papatsenko. The source * paper is Nazina, A. and Papatsenko, D. (2003). Statistical extraction * of Drosophila cis-regulatory modules using exhaustive assessment of * local word frequency, BMC Bioinformatics 4:65. * * * @author: $Author: zburke $ * @version: $Revision: 1.12 $ * @update: $Date: 2005/05/16 02:33:15 $ * */ public class Util { /** * return sum of input values * * @param c list of values to sum * * @return sum of values on list * */ public static double sum(Collection c) { long total = 0; Iterator i = c.iterator(); while (i.hasNext()) total += ((Long) i.next()).longValue(); return total; } /** * return sum of values on list * * @param c list of values to sum * * @return sum of values on list */ public static double sum(double c[]) { double total = 0; for(int i = 0; i < c.length; i++) total += c[i]; return total; } /** * compute probability distribution for a binomial distribution * allowing $mism mismatches over a sliding window of length $winSize, * * * @param n - sample size (length of string) * @param mismatch - number of mismatches permitted * @param winSize - sliding window size * * @return */ public static double binomial(long n, long mismatch, long winSize, double p) { double prob = 0.0; for (long x = 0; x <= mismatch; x++) { prob += winSize * ( StrictMath.pow(p,x) ) * ( StrictMath.pow((1.0 - p),( n - x )) ) * (factorial(n) / (factorial(x) * factorial(n - x))); } prob = 1 + 2 * prob; return prob; } /** * compute factorial of a number * * @param f number to compute factorial of * * @return factorial of f * */ public static double factorial(double f) { double value = 1; while (f > 0) { value *= f--; } return value; } /** * create frequency channels based on the number of desired channels * (channelCount) and their expected normal distribution around * expr. The value for each item on the returned list is the * maximum frequency contained by that channel. * * @param channelCount number of channels * @param expr * * @return */ public static long[] getChannels(long channelCount, double expr) { List channels = new ArrayList(); long limit_previous = 0; channels.add(new Long(0)); for (int i = 1; i < channelCount; i++) { int temp = new Double(0.5 + channelCount / 2).intValue(); double limChann = StrictMath.pow((i * (StrictMath.sqrt(expr) / temp)), 2); long limChannInt = StrictMath.round(StrictMath.floor(limChann)); if (limChannInt != limit_previous) { channels.add(new Long(limChannInt)); limit_previous = StrictMath.round(limChann); } } // convert the List channels into an array long[] channelArray = new long[channels.size() + 1]; int j = 0; Iterator i = channels.iterator(); while (i.hasNext()) channelArray[j++] = ((Long) i.next()).longValue(); return channelArray; } /** * create a 2D matrix * * @param cols * @param rows * @param value * * @return a new rows-by-cols matrix with value in every cell */ public static long [][] matrix(long cols, long rows, long value) { long matrix[][] = new long[new Long(rows).intValue()][new Long(cols).intValue()]; for (int i = 0; i < rows; i++) Arrays.fill(matrix[i], value); return matrix; } /** * smooth values on a list * * @param rate how much to smooth values * @param list values to smooth * * @return a copy of the list with its values smoothed */ public static double[] smooth(double rate, long list[]) { double smoothList[] = new double[list.length]; double item = 0.00; // current item double sum = 0.00; // current sum double count = 0; // item count int intCount = 0; for (int i = 0; i < list.length; i++) { if (i < rate) item = 1 + 2 * i; else if (list.length - i - 1 < rate) item = 2 * (list.length - i) - 1; else item = 2 * rate + 1; sum = 0; count = (1 - item) / 2; while (count <= (item - 1)/2) { intCount = (int) count; sum += list[i + intCount]; count++; } smoothList[i] = sum/item; } return smoothList; } /** * smooth values on a list * * @param rate how much to smooth values * @param list values to smooth * * @return a copy of the list with its values smoothed */ public static double[] smooth(double rate, double list[]) { double smoothList[] = new double[list.length]; double item = 0.00; // current item double sum = 0.00; // current sum double count = 0; // item count int intCount = 0; for (int i = 0; i < list.length; i++) { if (i < rate) item = 1 + 2 * i; else if (list.length - i - 1 < rate) item = 2 * (list.length - i) - 1; else item = 2 * rate + 1; sum = 0; count = (1 - item) / 2; while (count <= (item - 1)/2) { intCount = (int) count; sum += list[i + intCount]; count++; } smoothList[i] = sum/item; } return smoothList; } /** * return the largest value on a list * * @param list * @return */ public static long max(long[] list) { long max = list[0]; for (int i = 0; i < list.length; i++) if (list[i] > max) max = list[i]; return max; } /** * return the largest value on a list * * @param list * @return */ public static long max(List list) { Long max = (Long) list.get(0); Iterator i = list.iterator(); while (i.hasNext()) { Long l = (Long) i.next(); if (l.compareTo(max) > 0) max = l; } return max.longValue(); } /** * return input in FASTA format * * @param s string to format * @return FASTA-formatted string */ public static String fasta(String s) { int width = 70; int l = s.length(); StringBuffer fasta = new StringBuffer(); for (int i = 0; i < l; i += width) { if (i + width > l) fasta.append(s.substring(i) + "\n"); else fasta.append(s.substring(i, i + width) + "\n"); } return fasta.toString(); } /** * in a window from 0 - p.window_full in sequence, count the frequency * of word. * * @param p * @param sequence * @param globbackwin * @param word */ public static void freqCount(Model m, String sequence, long globbackwin, Word word) { int position = 0; while (position < m.windowFull + 1) { int offsetStart = (int)(position + globbackwin); int offsetEnd = (int)(offsetStart + m.length); char[] string = sequence.substring(offsetStart, offsetEnd).toLowerCase().toCharArray(); int[] mismatches = mismatchCount(string, word, m); if (mismatches[0] <= m.mismatchCount) word.frequency++; if (mismatches[1] <= m.mismatchCount) word.frequency++; position++; } } /** * count mismatches between two words, allowing p.mismatch_count * mismatches between the two strings. * * @param string * @param word * @param p * @return */ public static int[] mismatchCount(char string[], Word word, Model m) { int j = 0; int[] mismatches = new int[2]; for (int i = 0; i < string.length; i++) { if (string[i] != word.aWord[j]) mismatches[0]++; if (string[i] != word.aWordR[j]) mismatches[1]++; if (mismatches[0] > m.mismatchCount && mismatches[1] > m.mismatchCount) break; j++; } return mismatches; } /** * loop over the list of channels in order of frequency and stop when the * channel frequency is higher than the word frequency, i.e. stop at the * channel that contains the current word based on its frequency. return * this index int to the array. * * @param word - word with a frequency count * @param list - list of frequency channels * * @return array index of the frequency channel containing this word */ public static int limitCount(Word word, long[] list) { // TODO: throw an error here instead? if (list == null) return 0; for (int i = 0; i < list.length; i++) { if (word.frequency <= list[i]) return i; } // TODO: throw an error here instead? System.err.println("word frequency "+word.frequency +" is not contained in any channel"); for (int i = 0; i < list.length; i++) System.err.print (list[i] + " "); System.err.println(""); return 0; } /** * count the number of members that each frequency channel has in a given * window. channels is a list of channels for each word in the window. the * channels themselves are just values on an array, so an array with the * same dimsensions (sum) is used to store the member-counts for each channel. * * @param channels - list of channels for each word in a window * @param m.channelCount - how many channels to count * * @return */ public static List channelMemberCount(List channels, Model m) { ArrayList sum = new ArrayList(); int count = 1; Iterator i = null; while (count <= m.channelCount) { int members = 0; i = channels.iterator(); while (i.hasNext()) { if (((Long) i.next()).longValue() == count) members++; } // System.err.println("members: "+members); sum.add(new Long(members)); count++; } // System.err.println(""); return sum; } /** * detect channels for words in a window. * * @param s - DNA string to parse * @param pos - current position in the string * @param channels - list of channels detected so far * @param m - model for this sequence * * @return list of member-counts, per-channel */ public static List channelDetect(Sequence s, long pos, List channels, Model m) { List sum = null; int offsetStart = (int) pos; int offsetEnd = offsetStart + (int) m.length; Word word = new Word(s.sequence.substring(offsetStart, offsetEnd)); long globBackWin = pos - m.windowSmall; // count word's frequency in this sequence window. // store count in word.frequency Util.freqCount(m, s.sequence, globBackWin, word); // save the index of the channel containing this word int limitCount = Util.limitCount(word, m.channels); channels.add(new Long(limitCount)); // we've determined channels for all words; now store each channel's // member-count in sum. if (channels.size() == m.windowFull) { sum = Util.channelMemberCount(channels, m); channels.remove(0); } return sum; } /** * return sum of values in list between from and to, inclusive. * * @param list list of values to sum * @param from begin index * @param to end index * @return sum of values in list between from and to, inclusive. */ public static double sumArray(double[] list, int from, int to) { // TODO: throw an exception here instead if (to >= list.length) to = list.length - 1; double total = 0; for (int i = from; i <= to; i++) total += list[i]; return total; } /** * descending in-place sort of a list * * @param list to be reverse-sorted */ public static void reverseSort(double [] list) { // ascending sort Arrays.sort(list); // reverse it int max = list.length / 2; double j; for (int i = 0; i < max; i++ ) { j = list[i]; list[i] = list[list.length - i - 1]; list[list.length - i - 1] = j; } } /** * Plot outcome of a continuous random variable x. * * // TODO: figure out why sigma is outside the squareroot * compute a normal distribution (density function) * * @param x - * @param mu - mean * @param sigma - standard deviation * * @return */ public static double normalDistribution(double x, double mu, double sigma) { // return ((1 / (Math.sqrt(2 * Math.PI * sigma))) * // Math.pow(Math.exp(1), ((-1) * Math.pow((x - mu), 2) / (2 * Math.pow(sigma, 2))))); return ((1 / (sigma * Math.sqrt(4 * Math.atan2(1, 0)))) * Math.pow(Math.exp(1), ((-1) * Math.pow((x - mu), 2) / (2 * Math.pow(sigma, 2))))); } } Word.java package edu.dartmouth.bglab.crm.lwf; import java.lang.StringBuffer; /** * LWF helper class stores information about a substring * (word) within a sequence * * * This code is based on Perl written by Dmitri Papatsenko. The source * paper is Nazina, A. and Papatsenko, D. (2003). Statistical extraction * of Drosophila cis-regulatory modules using exhaustive assessment of * local word frequency, BMC Bioinformatics 4:65. * * * @author: $Author: zburke $ * @version: $Revision: 1.2 $ * @update: $Date: 2006/02/07 23:10:23 $ * */ public class Word { public String word = null; public String wordR = null; public char aWord[] = null; public char aWordR[] = null; public long frequency = 0; public Word(String w) { this.word = w.toLowerCase(); this.wordR = new StringBuffer(this.word).reverse().toString(); this.aWordR = this.wordR.toCharArray(); for (int i = 0; i < this.aWordR.length; i++) { if (this.aWordR[i] == 'a') this.aWordR[i] = 't'; else if (this.aWordR[i] == 't') this.aWordR[i] = 'a'; else if (this.aWordR[i] == 'g') this.aWordR[i] = 'c'; else if (this.aWordR[i] == 'c') this.aWordR[i] = 'g'; } this.wordR = new String(this.aWordR); this.aWord = this.word.toCharArray(); } } TrainProgressBar.java package edu.dartmouth.bglab.crm.lwf.train; import java.text.DecimalFormat; import java.text.NumberFormat; /** * LWF helper class has methods for displaying a progress * thermometer as parsing progresses. * * * This code is based on Perl written by Dmitri Papatsenko. The source * paper is Nazina, A. and Papatsenko, D. (2003). Statistical extraction * of Drosophila cis-regulatory modules using exhaustive assessment of * local word frequency, BMC Bioinformatics 4:65. * * * @author: $Author: zburke $ * @version: $Revision: 1.3 $ * @update: $Date: 2006/02/11 01:37:34 $ * */ public class TrainProgressBar { private String seqName; private long seqLen; private long inc; public long progress; private StringBuffer bar; private NumberFormat nf; // thermometer progress formatter // width of thermometer private final int WIDTH = 70; public TrainProgressBar() { this.bar = new StringBuffer(); } public TrainProgressBar(String seqName, long length) { this.bar = new StringBuffer(); init(seqName, length); } public void init(String seqName, long length) { this.seqName = seqName; this.seqLen = length; this.inc = length / 70; this.progress = this.inc; this.bar.append("."); this.nf = NumberFormat.getInstance(); if (this.nf instanceof DecimalFormat) ((DecimalFormat) this.nf).applyPattern("##.##"); System.out.println(""); System.out.println("Building statistics for "+seqName+"..."); } /** * update the progress thermometer * * print a thermometer that looks like this: * 4.25 [........... ] * the percentage and progress of the dots between the [ ] brackets is * regularly updated by appending \r instead of \n to the end of each * line * * @param i */ public void update(long i) { this.bar.append("."); this.progress += this.inc; System.out.print("\r" + this.bar + "\r" + ( 100 * i / this.seqLen )); StringBuffer pad = new StringBuffer(); int max = this.WIDTH - this.bar.length(); for (int j = 0; j <= max; j++) pad.append(" "); System.out.print( "\r ["+this.bar + pad +"]\r "+this.nf.format(100 * i / (double)this.seqLen ) +"\r"); // ^ ^ ^ ^ // ^ thermometer ^ percent complete so far // return to|start of line return to|start of line } } TrainSequenceParams.java package edu.dartmouth.bglab.crm.lwf.train; import java.lang.StrictMath; import java.io.InputStreamReader; import java.io.BufferedReader; import java.io.IOException; import edu.dartmouth.bglab.crm.lwf.Util; import edu.dartmouth.bglab.crm.lwf.Model; import edu.dartmouth.bglab.crm.lwf.LwfException; /** * LWF helper class; bit-bucket for runtime model-creation parameters. * * * This code is based on Perl written by Dmitri Papatsenko. The source * paper is Nazina, A. and Papatsenko, D. (2003). Statistical extraction * of Drosophila cis-regulatory modules using exhaustive assessment of * local word frequency, BMC Bioinformatics 4:65. * * * @author: $Author: zburke $ * @version: $Revision: 1.10 $ * @update: $Date: 2006/02/03 06:56:28 $ * */ public class TrainSequenceParams extends Model { /** * default constructor * */ public TrainSequenceParams() { } /** * constructor with some params * @param length * @param mmc * @param win * @param winSmall * @param cCount * @param files */ public TrainSequenceParams(long length, long mmc, long win, long winSmall, long cCount, String files) throws LwfException { set(length, mmc, win, cCount, files); this.windowSmall = winSmall; } /** * show information about this run * */ public void show() { System.out.println(""); System.out.println("---------------------------------------------------"); System.out.println("Expected frequency: "+this.expr); System.out.println("Channels created: "+this.channelCount); System.out.println("Channel limits: "+channelsString()); System.out.println("---------------------------------------------------"); System.out.println(""); } /** * return channel list as a string * * @return */ public String channelsString() { StringBuffer s = new StringBuffer(); for (int i = 0; i < this.channels.length; i++) s.append(this.channels[i] + " "); return s.toString(); } /** * gather input params * * @throws IOException if there is trouble reading params from command line */ public void get() throws IOException, LwfException { BufferedReader in = new BufferedReader(new InputStreamReader(System.in)); String line = null; long length = 0; long mmc = 0; long win = 0; long cCount = 0; String files = null; while (files == null || files.equals("")) { // word length System.out.print("\nSelect string length (1-12) [5]: "); line = in.readLine().trim(); length = line.equals("") ? 5 : new Integer(line).intValue(); // mismatch count long suggestmism = StrictMath.round(length / 3 - 0.1); System.out.print("\nSelect number of mismatches (0-4) ["+suggestmism+"]: "); line = in.readLine().trim(); mmc = line.equals("") ? suggestmism : new Integer(line).intValue(); // window size long suggestwin = StrictMath.round(64 * StrictMath.pow(2.0, suggestmism)); System.out.print("\nSelect detection window size (64-512) ["+suggestwin+"]: "); line = in.readLine().trim(); win = line.equals("") ? suggestwin : new Integer(line).intValue(); // channel count long suggestccount = 10; System.out.print("\nSelect maximal number of channels (2-24) ["+suggestccount+"]: "); line = in.readLine().trim(); cCount = line.equals("") ? suggestccount : new Integer(line).intValue(); // sequence file names System.out.print("\nEnter training sequence filenames (FASTA format): "); files = in.readLine().trim(); } set(length, mmc, win, cCount, files); } public void set(long length, long mmc, long win, long cCount, String files) throws LwfException { if (mmc >= length) throw new LwfException("mismatch count ("+mmc+") >= word-length ("+length+")"); // word length this.length = length; // mismatch count this.mismatchCount = mmc; // window size this.window = win; this.windowSmall = StrictMath.round(StrictMath.floor(this.window / 2)); this.windowFull = this.windowSmall * 2; // channel count this.channelCount = cCount; // sequence file names this.fileList = files; // sequence filenames this.trainsets = this.fileList.split("\\s"); this.expr = Util.binomial(this.length, this.mismatchCount, this.window, this.probability); System.err.println("expr: "+expr + "based on len/mism/win/prob: " + this.length +"/"+ this.mismatchCount+"/"+ this.window+"/"+ this.probability); // get channels this.channels = Util.getChannels( this.channelCount, this.expr ); this.channelCount = this.channels.length - 1; // add a final channel wide enough for every word on each // strand to be unique. this.channels[this.channels.length - 1] = 2 * this.window; } } TrainingSet.java package edu.dartmouth.bglab.crm.lwf.train; import java.io.BufferedReader; import java.io.FileNotFoundException; import java.io.FileReader; import java.io.IOException; import java.io.PrintStream; import java.util.Collection; import java.util.LinkedList; import java.util.List; import edu.dartmouth.bglab.crm.lwf.Sequence; import edu.dartmouth.bglab.crm.lwf.Util; import edu.dartmouth.bglab.crm.lwf.Model; /** * Stores run-time parameters about a set of paired positive and negative * training sequences that will later be used to scan an unknown sequence * to determine whether it more closely resembles the positive or negative * training set. * * * This code is based on Perl written by Dmitri Papatsenko. The source * paper is Nazina, A. and Papatsenko, D. (2003). Statistical extraction * of Drosophila cis-regulatory modules using exhaustive assessment of * local word frequency, BMC Bioinformatics 4:65. * * @author: $Author: zburke $ * @version: $Revision: 1.7 $ * @update: $Date: 2006/02/13 01:17:13 $ * */ public class TrainingSet { /** number of Sequences in this set */ public long count = 0; /** combined length of Sequences in this set */ public long length = 0; /** name of the file containing the raw sequence data */ public String input = null; /** length of each sequence */ public List lengths = null; /** Sequence objects in this training set */ public List sequences = null; /** * initializes lists * */ public TrainingSet() { init(); } /** * initializes lists and reads a sequence file * * @param filename FASTA sequence file to read * @param p * * @throws FileNotFoundException if sequence file cannot be found * @throws IOException if sequence file cannot be read */ public TrainingSet(String filename, Model m) throws FileNotFoundException, IOException { init(); sequenceFileParse(filename, m); } /** * initialize instance variables * */ private void init() { this.sequences = new LinkedList(); this.lengths = new LinkedList(); } /** * print sequence info to requested stream * * @param sequenceName * @param p * @param out */ public void headerPrint(String sequenceName, TrainSequenceParams p, PrintStream out) { out.print("Data: " + sequenceName + "\n"); out.print("Word: " + p.length + "\n"); out.print("Mism: " + p.mismatchCount + "\n"); out.print("Win: " + p.window + "\n"); out.print("Exp: " + p.expr + "\n"); out.print("Chann: " + p.channelCount + "\n"); out.print("Limits: " + p.channelsString() +"\n\n"); out.print("Statistics:\n\n"); } /** * parse a file containing FASTA data into a training set * * @param filename * @param p * * @throws FileNotFoundException if filename cannot be found * @throws IOException if filename cannot be read */ public void sequenceFileParse(String filename, Model m) throws FileNotFoundException, IOException { BufferedReader br = null; String line = null; int seqtotal = 0; Sequence sequence = new Sequence(filename); br = new BufferedReader(new FileReader(filename)); while ((line = br.readLine()) != null) { if (line.trim().equals("")) continue; // get new seq's header data if (line.matches(">.*")) { if (sequence.length > 0) { this.lengths.add(new Long(sequence.length)); this.sequences.add(sequence); } line.replaceAll(" ", "_"); sequence = new Sequence(filename); sequence.data = line; seqtotal++; } else if (line.matches("\\^\\^.*")) { if (sequence.length > 0) { this.lengths.add(new Long(sequence.length)); this.sequences.add(sequence); } sequence = new Sequence(filename); } else { sequence.sequence += line; sequence.length += line.length(); } } if (sequence.data == null) sequence.data = filename; this.lengths.add(new Long(sequence.length)); this.sequences.add(sequence); if (seqtotal == 0) seqtotal++; Sequence s = null; java.util.Iterator i = this.sequences.iterator(); while (i.hasNext()) { s = (Sequence) i.next(); s.maxPosition = s.length - m.windowSmall - m.length + 1; } this.count = seqtotal; this.input = filename; this.length = StrictMath.round(Util.sum((Collection) this.lengths)); } } ProgressBar.java package edu.dartmouth.bglab.crm.lwf.scan; import java.text.DecimalFormat; import java.text.NumberFormat; /** * LWF helper class has methods for displaying a progress thermometer * as parsing progresses. * * * This code is based on Perl written by Dmitri Papatsenko. The source * paper is Nazina, A. and Papatsenko, D. (2003). Statistical extraction * of Drosophila cis-regulatory modules using exhaustive assessment of * local word frequency, BMC Bioinformatics 4:65. * * * @author: $Author: zburke $ * @version: $Revision: 1.6 $ * @update: $Date: 2005/04/12 12:55:21 $ * */ public class ProgressBar { /** */ public long progress; private long seqLen; // total length of sequence private long win; // size of window in sequence private long inc; // unit of progress (progress proceeds in inc-sized increments) private StringBuffer bar; // the thermometer private NumberFormat nf; // thermometer progress formatter // width of thermometer private final int WIDTH = 70; /** * initialize the thermometer and display the sequence and model being * used for parsing. * * @param seqName name of sequence to parse * @param length * @param win * @param modelName name of model to parse */ public ProgressBar(String seqName, long length, long win, String modelName) { this.bar = new StringBuffer(); this.seqLen = length; this.win = win * 2; this.inc = (length - this.win) / this.WIDTH; this.progress = this.inc; this.bar.append("."); this.nf = NumberFormat.getInstance(); if (this.nf instanceof DecimalFormat) ((DecimalFormat) this.nf).applyPattern("##.##"); System.out.println(""); System.out.println(" processing "+seqName+" using model "+modelName+"..."); } /** * update the progress thermometer * * print a thermometer that looks like this: * 4.25 [........... ] * the percentage and progress of the dots between the [ ] brackets is * regularly updated by appending \r instead of \n to the end of each * line * * @param i */ public void update(long i) { this.bar.append("."); this.progress += this.inc; StringBuffer pad = new StringBuffer(); int max = this.WIDTH - this.bar.length(); for (int j = 0; j <= max; j++) pad.append(" "); System.out.print( "\r ["+this.bar + pad +"]\r "+this.nf.format(100 * i / ((double)this.seqLen - (double)this.win)) +"\r"); // ^ ^ ^ ^ // ^ thermometer ^ percent complete so far // return to|start of line return to|start of line } } ScanSequenceParams.java package edu.dartmouth.bglab.crm.lwf.scan; import java.io.BufferedReader; import java.io.BufferedWriter; import java.io.FileWriter; import java.io.IOException; import java.io.InputStreamReader; import java.io.PrintWriter; import java.util.List; import java.util.ArrayList; import java.util.Hashtable; import java.util.Iterator; import java.io.File; import edu.dartmouth.bglab.crm.lwf.LwfException; import edu.dartmouth.bglab.crm.lwf.Util; import edu.dartmouth.bglab.crm.lwf.Model; /** * LWF helper class; bit-bucket for runtime sequence scanning params. * * This code is based on Perl written by Dmitri Papatsenko. The source * paper is Nazina, A. and Papatsenko, D. (2003). Statistical extraction * of Drosophila cis-regulatory modules using exhaustive assessment of * local word frequency, BMC Bioinformatics 4:65. * * @author: $Author: zburke $ * @version: $Revision: 1.11 $ * @update: $Date: 2006/02/03 06:56:29 $ * */ public class ScanSequenceParams extends Model { /** filename of sequence to scan */ public String seqfname = null; /** output directory where parse files will be saved */ public String scanoutdir = null; /** string containing list of model files */ // public String fileList = null; /** list of +/- training model pairs */ public List stsList = null; public long stsmaxcount = 0; public List mdlwins = null; public List models = null; public List modelNames = null; public List channels = null; /** whether to write GNU plot data for a sequence recognition profile */ public boolean plot = false; /** * run-time param, freely enterable regardless of model-data * range is 0-100; default is 10 * passed as the standard-deviation to the normalDistribution method */ public long supression = 10; /** * threshold above which a window is classified positively. * Range: 0 - 1; 0 = stringent, 1 = relaxed */ public double profileCutoff = 0.0; /** * model-based param * minimum width of a positively-classified span that will be * considered a CRM */ public long peakWidth = 0; /** * model-based param * smoothing window */ public long inputSmooth = 0; /** * default constructor initializes lists * */ public ScanSequenceParams() { this.stsList = new ArrayList(); this.mdlwins = new ArrayList(); this.models = new ArrayList(); this.modelNames = new ArrayList(); this.channels = new ArrayList(); } /** * collect runtime params from the user for model parsing * * @throws IOException */ public void trainGet() throws IOException, LwfException { BufferedReader in = new BufferedReader(new InputStreamReader(System.in)); String line = null; // sequence name System.out.print("\nEnter sequence name: "); String sequenceName = in.readLine().trim(); // model names System.out.print("\nEnter +/- model pairs: "); String models = in.readLine().trim(); // equalize channels System.out.print("\nEqualize channels (1-100) [10]: "); line = in.readLine().trim(); int supression = line.equals("") ? 10 : new Integer(line).intValue(); trainSet(sequenceName, models, supression); } /** * save runtime * @param seqfname * @param fileList * @param supression * @throws LwfException */ public void trainSet(String seqfname, String fileList, int supression) throws LwfException { // sequence name this.seqfname = seqfname; // make sure we got pairs of models String[] stsList = fileList.split("\\s"); if (stsList.length % 2 != 0) throw new LwfException("you must specify positive and negative training model pairs"); // create the output directory outDirCreate(); for (int i = 0; i < stsList.length; i += 2) { Hashtable hash = new Hashtable(); hash.put("pos", stsList[i]); hash.put("neg", stsList[i + 1]); this.stsList.add(hash); } // equalize channels this.supression = supression; } /** * collect runtime params for scanning * * @throws IOException */ public void scanGet() throws IOException, LwfException { BufferedReader in = new BufferedReader(new InputStreamReader(System.in)); String line = null; // profile cutoff double cutoffSuggest = profileCutoffSuggest(); System.out.print("\nSelect profile cutoff (0-1) [" + cutoffSuggest + "]: "); line = in.readLine().trim(); double cutoff = line.equals("") ? cutoffSuggest : new Double(line).doubleValue(); // peak width cutoff (min CRM width) // suggestion is 1/2 size of biggest channel long peakSuggest = peakWidthSuggest(); System.out.print("\nSelect peak width cutoff [" + peakSuggest + "]: "); line = in.readLine().trim(); long width = line.equals("") ? peakSuggest : new Integer(line).intValue(); // how much to smooth data // suggestion is 1/2 size of biggest channel long smoothSuggest = smoothSuggest(); System.out.print("\nSelect smoothing window [" + smoothSuggest + "]: "); line = in.readLine().trim(); long smooth = line.equals("") ? smoothSuggest : new Integer(line).intValue(); // whether to write GNUPlot files System.out.print("\nSelect plot option (y/n) [" + plotSuggest() + "]: "); line = in.readLine().trim(); boolean gnuplot = line.equals("y") ? true : false; scanSet(cutoff, width, smooth, gnuplot); } public double profileCutoffSuggest() { return 0.3; } public long peakWidthSuggest() { return Util.max(this.mdlwins) / 2; } public long smoothSuggest() { return peakWidthSuggest(); } public boolean plotSuggest() { return false; } /** * set run-time scanning parameters and stash a README file in the * output directory that records them. * * @param profile * @param peak * @param smooth * @param gnuplot * * @throws LwfException in case the README can't be written */ public void scanSet(double profile, long peak, long smooth, boolean gnuplot) throws LwfException { // profile cutoff this.profileCutoff = profile; // peak width cutoff (max channel width???) this.peakWidth = peak; // how much to smooth data this.inputSmooth = smooth; // whether to write GNUPlot files this.plot = gnuplot; // put a readme in the outdir that specifies the models outDirReadme(); } /** * Create a uniquely named output directory in the pattern *
$number |
$number | \n";
print FILE "|
| phi: $round->{'phi'} $round->{'motifs'} |
EOT
}
print GENE "
| ". (join " | \n\t", @{ $summary{$row} } ) . " | \n
| $block->{'density'} |