I have a final report methylation data for ~800k rows for CpG sites for 50 samples. There are many missing-values (NA) in this dataset which I have a problem when I want use Normalization and also differential methylation analysis in Limma. I just provide you a small sample of my data as follows: (Values of each cell is Beta-values)

My question is: How to handle these missing-data (NA) in this huge dataset (~800k rows + 50 columns)? Is there any package in R to consider missing data? Is there a fast program to impute missing data in R? Thanks in advance for any advise

Well, is your matrix numeric? This works fine for me:

CpG <- matrix(rnorm(1000), ncol=10, nrow=100)
CpG[sample(length(CpG), 100)] <- NA # Adding some missing values.
design <- model.matrix(~gl(2,5))
fit <- lmFit(CpG, design)

You'll have to be more precise about the problems you're having with normalization. As far as I am aware, the normalization methods in limma can deal with missing values. (Assuming, of course, that it makes sense to apply them to methylation data - I'm not familiar enough with this to be sure.) In any case, you shouldn't be using beta values for linear modelling, see https://dx.doi.org/10.1186/1471-2105-11-587.