Demosaicing and noise reduction should be done simultaneously. It's best to use the DCRaw program to extract the raw pixel data and then to tackle the problem of reconstructing the image from first principles instead of using the standard algorithms which won't yield the best possible reconstruction of the image.
In such a first principle treatment, you frame the problem as a so-called inverse problem where you calculate the most likely picture given incomplete information in the form of the sensor data. You then need to build a model based on measurements that gives the probability P(Y|X) for getting sensor data Y given that the real image data is X. Once you have a good model for this function, you want to find the probability P(X|Y) for the probability that the image is X given that the sensor data is Y. The relation between the two is given by Bayes' theorem and it involves the prior probability P(X) for the image. This prior can be described using a simple model appropriate for case at hand.
If you choose your models simple enough e.g. Gaussian noise, Gaussian priors, with only local correlations ( e.g. only correlations between between nearest and next nearest neighbors), then the model will be exactly solvable. It is then possible to write down the image X that maximizes P(X|Y) in a pixelwise fashion, the solution for the optimal pixel value at some point is then given as a weighted summation over the sensor data. Such models are too simplistic to yield very good results, but they already are quite good compared to what the raw processors have to offer.
E.g. I've implemented a simple noise filter that is based on defining P(X) on the basis of how least squarest deviations of 3rd degree polynomial fits to 6 by 6 blocks. So, given some arbitrary picture defined by pixels values p_{n,m} chop it up into 6 by 6 blocks and write down the formal expression of the sum of all the least squarest deviations when fitting the polynomials to each of the blocks in terms of p_{n,m}. The expression S you get is a quadratic function of all the p_{n,m}. I then assume that the probability P(X) is proportional to exp(-a S) where a is some parameter that can be fixed later. Because S is a quadratic expression, it is a Gaussian model and can be solved exactly (when imposing periodic boundary conditions, translational invariance allows you to solve the model using Fourier transforms).
The result is then a filter that ends up reducing the noise averaging over neigboring pixels, but it does so using both positive and negative weights. The reason the negative weight enter is to undo the effects of blurring that you would get if all the weights were positive (it's similar to how the unsharp mask works using negative weights). But this comes out of the math, you don't put anything in here on an ad hoc basis, other than your assumption of the prior P(X) and a noise model.
But you can do much better by writing down realistic models tuned to whatever scene you are photographing and the sharpness. The picture that maximizes P(X|Y) then won't have an analytic solution but you can then solve it using numerical methods involving successive approximations that get better and better.
