beam2.m

The Matlab function beam2.m takes as input the observed gene expression levels of a single gene in two repeated experiments, and outputs the BEAM estimate of the true transcript level, and the standard deviation of this estimate.

Required data files:

two_val_est.mat

two_val_std.mat

x_range.mat

function [est, std] = beam2(X1, X2)
% [EST, STD] = BEAM2(X1, X2)
% Input is two values X1 and X2 which represent readings of the same gene
% from two repeated experiments.  EST returns the BEAM estimate of the
% underlying transcript level.  STD is the standard deviation of the estimate.
% X1 and X2 may be vectors or matrices of the same size, in which case
% EST and STD are also the same size and return the estimates for the
% corresponding elements.

if any(size(X1) ~= size(X2))
  error('The two input matrices must have the same size');
end

% Load lookup table
load -ascii two_val_est;
load -ascii two_val_std;
load -ascii two_val_range;

min_range = min(two_val_range);
max_range = max(two_val_range);
if any(min(X1(:),X2(:)) < min_range) | any(max(X1(:),X2(:)) > max_range)
  warning(sprintf('One or more input values is outside the range \n of the precomputed lookup table (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

[m_x, m_y] = meshgrid(two_val_range, two_val_range);

est = interp2(m_x, m_y, two_val_est, X1, X2, 'linear');
std = interp2(m_x, m_y, two_val_std, X1, X2, 'linear');