The Matlab function beamr.m takes as input the observed gene expression levels of a single gene in two experiments under different conditions, and outputs the BEAM estimate of the ratio of the transcript levels, along with the probability that a ratio at least that extreme would be observed if the true transcript levels were identical.

Required data files:

  • rat_val_est.mat
  • rat_val_std.mat
  • x_range.mat
  • p.mat
  • p_range.mat

    function [ratio, p_value, ratio_std_factor] = beamr(X, Y)
    % Input is two values X and Y which represent two readings of the same gene
    % from experiments under different conditions.  RATIO returns the BEAM
    % estimate of the transcript level ratio of X to Y.  P is the significance
    % value showing the probability that the same or a more extreme ratio would
    % be observed if the underlying transcript levels were the same.
    % STD_FACTOR is the standard deviation of the ratio estimate, reported as
    % a factor (i.e., the 1-sigma interval for the true ratio is
    % X and Y may be vectors or matrices of the same size, in which case
    % RATIO, P, and STD_FACTOR are also vectors of the same size and return 
    % the BEAM ratio and p-value for corresponding elements.
    if any(size(X) ~= size(Y))
      error('The two input matrices must have the same size');
    % Load lookup tables
    load rat_val_est
    load rat_val_std
    load rat_val_range
    [r_x, r_y] = meshgrid(rat_val_range, rat_val_range);
    % Estimate ratios.  For reasons of symmetry and numerical accuracy, our
    % lookup table stores Bayes least squares estimates of log ratios in
    % rat_val_est.  To obtain ratio estimates, we exponentiate these values,
    % after performing the appropriate interpolation.
    ratio = 10.^(interp2(r_x, r_y, rat_val_est, X, Y, 'linear'));
    % Find standard deviations of ratio distributions.  Again, for reasons of
    % symmetry and numerical accuracy, rat_val_std stores standard deviations of
    % the log ratio distributions.  The interpolated values are stored in
    % log_ratio_std.  We exponentiate these values to compute
    % ratio_std_factor.
    log_ratio_std = interp2(r_x, r_y, rat_val_std, X, Y, 'linear');
    ratio_std_factor = 10.^(log_ratio_std);
    load p_val;
    load p_range;
    % bounds checking
    min_range = max(min(rat_val_range), min(p_range));
    max_range = min(max(rat_val_range), max(p_range));
    if any(min(X(:),Y(:)) < min_range) | any(max(X(:),Y(:)) > max_range)
      warning(sprintf('One or more input values is outside the range \n of the precomputed lookup tables (i.e., less than %1.0f \n or greater than %1.0f).  The corresponding output \n values will be stored as ''NaN'' (Not a Number)',min_range,max_range));
    [p_x, p_y] = meshgrid(p_range, p_range);
    p_value = interp2(p_x, p_y, p_val, X, Y, 'linear');