# beamr.m

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)
% [RATIO, P, 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
% [RATIO/STD_FACTOR RATIO*STD_FACTOR].
%
% 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');
end

[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);

% 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));
end

[p_x, p_y] = meshgrid(p_range, p_range);

p_value = interp2(p_x, p_y, p_val, X, Y, 'linear');

```

• rif@alum.mit.edu