I can suggest a simple way that is not perfect but is reasonable and not too difficult.
Transaction Costs Will Change the Optimal Portfolio
What you want to do is account for the transactions costs of trading. Note that Engle and Ferstenberg (2007) have shown that including transactions costs shifts the efficient frontier down and to the right (because there is an added mean and variance of transactions costs).
This will lead to a different optimal portfolio than if you ignore transactions costs. A common complaint about Markowitz portfolio optimization (without transactions costs) is that the portfolios underperform "out-of-sample." Including transactions costs (mean and variance) may reduce the underperformance problem.
Approximation for "Small" Trades: Using Only Bid-Ask Spread
For "small" trade amounts, you can approximate the mean and variance of trading as being affected only by the percentage spread $S_i$ for asset $i$
$$
S_i = \log(A_i/B_i)
$$
where $A_i$ is the ask price and $B_i$ is the bid price.
Also, you must account for the cost of entry and exit -- so we multiply the single-trade mean $E(c_i)$ and variance $\textrm{var}(c_i)$ of cost by 2.
If you are patient and wait for a possibly better price before giving up and crossing the spread (buying at the offer or selling at the bid), this would be:
- mean: half the mean bid-ask spread per trade; total = $2\times\bar{S}_i/2=\bar{S}_i$, and
- variance: cf Rademacher random variable; total = $2\times\bar{S}_i^2/4=\bar{S}_i^2/2$.
If you are impatient, your execution price is much more certain (and worse):
- mean: the full spread per trade; total = $2\times\bar{S}_i$, and
- variance: 0.
You then subtract each of these mean costs from your positive expected returns and add it to the negative expected returns and add the cost variances to your return variances. Or, equivalently:
$$
E(r_{i,new}) = E(r_i) - 2\textrm{sign}(r_i)E(c_i), \\
\textrm{var}(r_{i,new}) = \textrm{var}(r_i) + 2\textrm{var}(c_i).
$$
Note that you should not let the expectations cross through 0. If your expected return is 5% and the expected entry+exit cost is 6%, you are better off treating this like it has a 0% expected return rather than a -1% expected return (which might lead to shorting the stock).
Should you propagate the cost variances through to covariances, i.e. create covariances from correlations and (new) variances? That depends on whether you think execution costs are correlated or not. If you want execution costs to be correlated, propagate the uncertainty; if not, only alter your variances and not the covariances.
Approximation for Larger Trades: Using a Price Impact Model
You mention wanting to use bid-ask spreads; however, I will also hint at a better approximation. Basically, we can assume your transactions costs have a few parts and use a price impact model to get a better sense of costs.
While there are lots of price impact models, you are not doing trade scheduling at this step -- so you can use a model that might be unsuitable for trade scheduling but an OK approximation for portfolio optimization.
For trading a quantity $Q_i$ of asset $i$, I would use something like the Torre and Ferrari (1997) model:
$$
E(c_i) = \alpha_i + \beta_{iS} \bar{S}_i + \beta_{iQ} \sigma_i\sqrt{\frac{|Q_i|}{\bar{V}}}
$$
where
- $\alpha_i$ captures fixed costs (like a commission fee),
- $\beta_{iS}$ gets at how much of the log-spread you typically pay,
- $\beta_{iQ}$ estimates the permanent impact (and cumulated decaying impact as per Obizhaeva and Wang (2013)) paid, and
- $\bar{V}$ is the average volume per trading period.
The variance calculation is tricky: you want to incorporate the uncertainty of your model coefficients with the uncertainty of the model terms. If I had to guess at this, I would expect something like
$$
\textrm{var}(c_i) = \sigma_\pi^2 + \hat\beta_{iS}^2 \frac{\bar{S}^2}{4} + \hat\beta_{iQ}^2 \left[
\sigma_i^2\frac{n-1}{2}\left(\frac{\Gamma(n/2)}{\Gamma(\frac{n-1}{2})}\right)^2 \frac{Q_i^2}{\bar{V}^2} +
\sigma_i^2 Q_i^2\frac{\textrm{var}(\bar{V})}{\bar{V}^4}\right]
$$
where $\sigma_\pi$ is the standard deviation from a fitted model prediction interval and $n$ is the number of observations used to estimate $\sigma_i$. The last term comes from the variance of $\sigma_i$, a one-term delta method approximation for $\textrm{var}(1/\bar{V}))$, and the variance of a product of random variables.
With these (hopefully better) estimates of $E(c_i)$ and $\textrm{var}(c_i)$ you can then double them (cost of entry and cost of exit) and add them in as done for the simpler small-trade case.