I have a Gene-Expression file contain 18,000 genes (RPKM) for 150 patients of Breast Cancer (Normalized).

I would like to calculate the first 4 principle component (PCA1, PCA2, PCA3,PCA4) of this data and then insert them in the model along with "gender" and "age" for adjusting in Limma. Could you please advise how to do that.

design=model.matrix(~0+age+gender+PCA1+PCA2+PCA3+PCA4+treatment)

Let's assume your matrix of RPKM is set up with genes in rows and samples in columns, which is how these things usually go. Further, assume this matrix is called 'my.rpkms'.

pca <- prcomp(t(my.rpkms))
design <- model.matrix(~0+age+gender+treatment)
design <- cbind(design, pca$x[,1:4])