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 % 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)); end [p_x, p_y] = meshgrid(p_range, p_range); p_value = interp2(p_x, p_y, p_val, X, Y, 'linear');