One way to do it would be to center your data on (x_0, y_0), run linear regression (specifying no intercept), and then transform the predictions back to the original scale.
Here's an example:
import numpy as np
import matplotlib.pyplot as plt
x = np.random.randn(100)
y = np.random.randn(100) + x
plt.plot(x, y, "o")
plt.plot(x[0], y[0], "o") # x_0, y_0 is the orange dot

Next, use linear regression without an intercept to fit a line that goes through the transformed data.
from sklearn.linear_model import LinearRegression
lm = LinearRegression(fit_intercept = False)
# center data on x_0, y_0
y2 = y - y[0]
x2 = x - x[0]
# fit model
lm.fit(x2.reshape(-1, 1), y2)
Last, predict the line and plot it back on the original scale
# predict line
preds = lm.predict(np.arange(-5, 5, 0.1).reshape(-1,1))
# plot on original scale
plt.plot(x, y, "o")
plt.plot(x[0], y[0], "o")
# add x_0 and y_0 back to the predictions
plt.plot(np.arange(-5, 5, 0.1) + x[0], preds + y[0])
