How can I perform a DC simulation of a two-stage op-amp design in NGSpice to determine if my transistors are operating in saturation? Is there a direct method to identify this, and what are the steps to achieve this?

I don't think there is a direct method in this. What I did is I just measured the Vds, Vgs, and Vth then did a saturation test from Vds >/ |Vgs - Vth| manually. Also curious if there are other method to do this.

Thats great, How did you print those value inside schematic window? I am new to Ngspice and xschem