I know that $$\int_0^T f(B_t, t)dB_t=\lim_{n\to \infty }\sum_{i=1}^n f(B_{t_i^{(n)}},t_i^{(n)})(B_{t_{i+1}^{(n)}}-B_{t_i^{(n)}}),\quad \text{in }L^2,$$ where $\{t_i^{(n)}\}_{i=1}^n$ is a sequence of partition of $[0,T]$ s.t. $$\sup|t_{i+1}^{(n)}-t_i^{(n)}|\to 0,\quad \text{when }n\to \infty .$$
Now, there is a subsequence $n_k$ s.t. $$\lim_{k\to \infty }\sum_{i=1}^{n_k}f(B_{t_i^{(n_k)}},t_i^{(n_k)})(B_{t_{i+1}^{(n_k)}}-B_{{t_i}^{(n_k)}})=\int_0^T f(B_t,t)dB_t\quad a.s.$$
So, up to a subsequence, can the stochastic integral be seen as a Stiljes-Riemann integral ? A sort of week Stiljes-Riemann integral in the sense that indeed $$\lim_{n\to \infty }\sum_{i=1}^{n}f(B_{t_i^{(n)}},t_i^{(n)})(B_{t_{i+1}^{(n)}}-B_{{t_i}^{(n)}}),$$ may not exist, but if it exist, then it converges to $\int_0^T f(B_t,t)dB_t$ a.s. Does this makes sense ?
The convergence of the subsequence holds almost surely, i.e. there is an exceptional null set where convergence fails to hold. This null set depends on the partitioning sequence $(t^{(n)})_{n \geq 1}$. Since there are uncountably many sequences, this is pretty bad - the union of the null exceptional sets is going to be quite huge, in general, and therefore we cannot expect to have a "uniform" exceptional null set for all sequences $(t^{(n)})_{n \in \mathbb{N}}$ with mesh size converging to zero. In fact, that's exactly the reason why the Itô integral is defined as an $L^2$-limit of the Riemann sums and not as a pointwise limit.
The phenomena which you are observing is a very general one: If a sequence of random variables $Y_n$ converges to a random variable $Y$ in probability (or in $L^2$), then we can choose a subsquence which converges almost surely to $Y$. Nevertheless, convergence in probability is a much weaker notion of convergence than pointwise convergence; this means, in particular, that the pointwise convergence of a subsequence of Riemann sums is far from giving a notion for a pointwise integration: If we want to get a pointwise notion for a stochastic integral, then we would like to fix $\omega \in \Omega$ and then define the stochastic integral, say, as a pointwise limit of Riemann sums along a suitable partition. That's, however, not what happens if we use the subsequence procedure from your question. If we take a sequence $(t^{(n)})_n$ then we get pointwise convergence with probability $1$ but we have no control about the null set. In particular, we don't have a clue how to choose a sequence $(t^{(n)})_n$ such that the Riemann sums convergence for our fixed $\omega$.
Let me give one further remark. In order to get the convergence
$$\int_0^T f(s) \, dB_s = \lim_{n \to \infty} \sum_{i=1}^n f(t_i^{(n)}) (B_{t_{i+1}^{(n)}}-B_{t_i^{(n)}}) \quad \text{in $L^2$}$$
(and hence the pointwise convergence of the subsequence of Riemann sums) you will typically need some continuity assumptions on $f$. If $f$ is a general progressively measurable function with $\mathbb{E}\int_0^t f(s)^2 \, ds < \infty$ for all $T>0$, then there exists some sequence of approximating simple functions $(f_n)_{n \in \mathbb{N}}$, i.e. a sequence of functions such that$$\mathbb{E}\int_0^t |f_n(s)-f(s)|^2 \ ds \to 0 \quad \text{and} \quad \int_0^t f(s) \, dB_s = L^2-\lim_{n \to \infty} \int_0^t f_n(s) \, dB_s;$$
however, the approximating functions $f_n$ will be, in general, not of the form
$$f_n(s) := \sum_{i=1}^n f(t_i^{(n)}) 1_{[t_i^{(n)},t_{i+1}^{(n)})}(s)$$
(which would give rise to the Riemann sums you are stating at the very beginning of the question); see e.g. Proposition 15.16 and Lemma 15.19/Theorem 15.20 in the book by Schilling & Partzsch (2nd edition) for more information.