In the investigation of fatigue properties of metals, the microstructure-based modelling has shown its powerful applicability in predicting the microcrack initiation as well as the fatigue life. However, proper treatment of the inclusions, which are the major fatigue crack trigger especially for the very high cycle fatigue regime, is still missing. It is emphasised that in addition to the geometrical representation and the basic mechanical properties assignment of the inclusions, the residual stresses developed between the steel matrix and inclusions during the cooling processes due to their distinct thermal expansion coefficients play a non-negligible role in determining the fatigue properties. Therefore, it is aimed, in this study, to propose a microstructure-based modelling approach to account for the effects of residual stresses induced by the rapid cooling process on the fatigue crack initiation behaviour of a martensitic steel, for which the majority of the fatigue crack is formed around the calcium aluminate inclusions in experiments. The entire approach is decomposed into two processes: i) simulation of the cooling process to obtain the residual stress profile around the inclusion and ii) conducting the fatigue simulation using a crystal plasticity model including the mapped residual stress profile from the previous step. It is shown that the proposed approach accurately predicts the fatigue crack initiation sites around the inclusions corresponding to the experimental findings, while the modelling approach without the residual stresses fails to predict the correct locations of the crack initiation, revealing the necessity to consider the residual stresses for the future fatigue modelling and assessment.