Please help me to calculate the value of resistor R1 so that the current through R2 was 100uA. The only thing I noticed is that the output current is pretty much sensitive to changes of saturation current Is of both diode and transistor.
If you assume that the junctions are the same (patently false, and this dubious assumption must be stated), then the current through R2 ~= current through R1 since \$\beta>>1\$ . You know the voltage at the base of the transistor within a reasonable margin, so you can calculate the current through R1. The value of R2 is not important unless you take 2nd-order effects such as the Early voltage into account. To see how far off the initial assumption is in reality, let's simulate it with a 100K resistor:
The current through R1 is 145uA and the current through R2 is 360uA according to the simulation. So the error is about 2:1 and there's really no point bothering about the 0.3% or whatever due to hFE.
This error originates from the difference in saturation currents of the two junctions and also the ideality factor of the diode junction, which is 1.45 in the above simulation. If I replace the diode with a diode-connected 2N2222, the currents will be much closer:
Now the currents are 145uA and 169ua so the error is only about 15%, but hFE is still an effect of little significance. The error at this point is mostly due to the Early voltage.
Note that you have to use a reasonably accurate transistor model to get reasonably accurate results. In this case Circuitlab (and also the free LTSpice) use good models (Ebers-Moll) with good default parameters and give consistent results showing about 15-20% error.
What you have there is a circuit that, when you increase the current through R1, the current through R2 will increase.
However, that's the only thing that can be said about it for certain. The actual ratio of R2 to R1 current is poorly controlled. As you say, it depends on the detail of the diode and transistor base-emitter conduction.
It's easy enough to arrive at a current through R1 to give any R2 current, if you make assumptions about D1 and Q1 parameters, and the temperature. It will be dependent on which model the simulator uses for the transistor. Be aware that this mode is almost never used, well, certainly never trusted in this region, so the models may well be quite inaccurate, as they don't need to be accurate, as nobody would trust them here anyway. If you give an answer based on your simulator, make sure you quote the transistor model, and its parameters, in your answer.
You might well, for the sake of your answer, assume the diode conduction curves for D1 and Q1 are identical if you like, as long as you state the assumption.
To make a proper (robust with temperature, and various device types) current source, place a resistor in series with the diode, and a resistor in series with Q1 emitter. Aim to drop around one volt at the design output current. This will be sufficient to swamp the differences in diode voltage for most 'biassing' purposes, though a precision current source it still isn't.
For 100uA output, put 10k in Q1 emitter. Now you can choose the D1 series resistor to scale the output current if you like. A 100k resistor would result in an output current of 10x I(R1), a 1k resistor of 0.1xI(R1).
I(R1) will equal I(R2) when the same transistor is used for the diode within 2% unless Hfe of Q2 is very large (>500) Here simulated with hFE=100 for Q2