Here is some code that might do the job. I am assuming that your image is truecolor, so the data have dimensions M x N x 3 for some M and N, and the colors are 24-bit, so the data type in R is uint8. If your images are different, you'll have to adjust the code:
n = length(dirlist);
R = imread([pathname, dirlist(1).name]);
RR = zeros([size(Rav) n]);
% Put the images into an M x N x 3 x 60 array
for m2 = 2:n
RR(:,:,:,m2) = imread([pathname, dirlist(m2).name]);
end
Rav = mean(RR,4); % average along the 4th dimension of the data
RR = bsxfun(@minus,RR,Rav); % subtract the average from each image
You can examine the kth image with
R = RR(:,:,:,k);
image(R)