Just as the paper says
$$L_i(\theta_i)= E_{(s,a)\sim p}[(y_i-Q(s,a;\theta_i))^2]$$
where
$$y_i = E_{s' \sim \mathcal{E}}[r+\gamma \max_{a'}Q(s',a';\theta_{i+1})\mid s,a]$$
Then in the Background section of the paper, it says
Differentiating the loss function with respect to the weights we arrive at the following gradient:
$$\nabla_{\theta_i} L_i(\theta_i)\\= E_{(s,a)\sim p,s'\sim\mathcal{E}}\left[\left(r+\gamma \max_{a'}Q(s',a';\theta_{i+1})-Q(s,a;\theta_i)\right)\nabla_{\theta_i}Q(s,a;\theta_i)\right]\tag{1}$$
Rather than computing the full expectations in the above gradient, it is often computationally expedient to optimize the loss function by stochastic gradient descent.
...
and the expectations are replaced by single samples from the behavior distribution $ρ$ and the emulator $\mathcal{E}$ respectively.
If you're familiar with SGD and Stochastic Optimization then you know what happens here:
The expression inside the expectation of (1), i.e. $\left(r+\gamma \max_{a'}Q(s',a';\theta_{i+1})-Q(s,a;\theta_i)\right)\nabla_{\theta_i}Q(s,a;\theta_i)$ , is an unbiased estimation of the real gradient $\nabla_{\theta_i}L_i(\theta_i)$ - its expectation is the real gradient. In other words,
$$\widehat{\nabla L}=\left(r+\gamma \max_{a'}Q(s',a';\theta_{i+1})-Q(s,a;\theta_i)\right)\nabla_{\theta_i}Q(s,a;\theta_i).$$
Then by theory of Stochastic Optimization we can optimize $L$ by $\theta\leftarrow \theta - \alpha\widehat{\nabla L}$ , which is how SGD works.
The unbiased estimation $\left(r+\gamma \max_{a'}Q(s',a';\theta_{i+1})-Q(s,a;\theta_i)\right)\nabla_{\theta_i}Q(s,a;\theta_i)$ can be sampled and calculated directly - you can run the emulator $\mathcal{E}$ and current behavior policy to collect $r, s', a'$ and calculate the gradient of your Q network $\nabla Q$ using TensorFlow. So the gradient to $L$ can be approximated.
The $y_j$ in the pseudocode of the paper is also an estimation. It'll be more proper to denote it as $\hat{y_j}$.
(I'm not an native English speaker so forgive my poor expression.)