import numpy as np import plotly.graph_objs as go from sklearn.covariance import MinCovDet from sklearn.neighbors import LocalOutlierFactor from scipy.stats import chi2 from statsmodels.robust.scale import mad # ---------------------------- # Config toggles # ---------------------------- SHOW_BOOTSTRAP_CENTER = True N_BOOT = 500 OUTLIER_ALPHA = 0.95 ELLIPSOID_ALPHA = 0.95 MCD_SUPPORT = 0.75 AXIS_Z_THRESHOLD = 2.5 LOF_NEIGHBORS = 7 RANDOM_SEED = 42 ALLOWED_3D = {'circle','circle-open','cross','diamond','diamond-open','square','square-open','x'} def symbol3d(s): return s if s in ALLOWED_3D else 'cross' # ---------------------------- # Your data # ---------------------------- # ---------------------------- #[add here] # ---------------------------- # Helpers # ---------------------------- def make_ellipsoid_surface(center, cov, alpha=0.95, n_u=60, n_v=30): df = 3 chi2_val = chi2.ppf(alpha, df) eigvals, eigvecs = np.linalg.eigh((cov + cov.T) / 2.0) eigvals = np.clip(eigvals, 1e-9, None) radii = np.sqrt(chi2_val) * np.sqrt(eigvals) u = np.linspace(0, 2*np.pi, n_u) v = np.linspace(0, np.pi, n_v) xs = np.outer(np.cos(u), np.sin(v)) ys = np.outer(np.sin(u), np.sin(v)) zs = np.outer(np.ones_like(u), np.cos(v)) E = np.stack([radii[0]*xs, radii[1]*ys, radii[2]*zs], axis=0).reshape(3, -1) R = eigvecs E_rot = R @ E X = (E_rot[0, :] + center[0]).reshape(xs.shape) Y = (E_rot[1, :] + center[1]).reshape(xs.shape) Z = (E_rot[2, :] + center[2]).reshape(xs.shape) return X, Y, Z def robust_mcd(X, support_fraction=MCD_SUPPORT, random_state=RANDOM_SEED): mcd = MinCovDet(support_fraction=support_fraction, random_state=random_state).fit(X) center = mcd.location_ cov = (mcd.covariance_ + mcd.covariance_.T) / 2.0 w, V = np.linalg.eigh(cov) w = np.clip(w, 1e-9, None) precision = V @ np.diag(1.0 / w) @ V.T diffs = X - center D2 = np.einsum('ij,jk,ik->i', diffs, precision, diffs) return mcd, center, cov, precision, D2 def bootstrap_centers(X, n_boot=500, random_state=RANDOM_SEED, support_fraction=MCD_SUPPORT): rng = np.random.default_rng(random_state) centers = [] for _ in range(n_boot): idx = rng.integers(0, X.shape[0], size=X.shape[0]) Xm = X[idx, :] try: mcd = MinCovDet(support_fraction=support_fraction, random_state=rng.integers(0, 1_000_000)).fit(Xm) centers.append(mcd.location_) except Exception: centers.append(np.median(Xm, axis=0)) return np.array(centers) # ---------------------------- # Build matrices and robust fit # ---------------------------- X = np.column_stack([or_genes, gr_genes, ir_genes]).astype(float) mcd, center, cov, precision, D2 = robust_mcd(X) df = 3 thr = chi2.ppf(OUTLIER_ALPHA, df) thr975 = chi2.ppf(0.975, df) order = np.argsort(D2)[::-1] print("Top Mahalanobis D² (largest first):") for idx in order[:10]: print(f"{species[idx]:22s} D²={D2[idx]:.3f} | >95%? {int(D2[idx] > thr)} | >97.5%? {int(D2[idx] > thr975)}") is_outlier = D2 > thr inlier_idx = np.where(~is_outlier)[0] out_idx = np.where(is_outlier)[0] vmin = np.min(D2) vmax = np.max(D2) med = np.median(X, axis=0) mad_vals = np.array([mad(X[:, j], center=med[j]) for j in range(3)]) robust_z = (X - med) / (mad_vals + 1e-12) axis_outlier = np.any(np.abs(robust_z) > AXIS_Z_THRESHOLD, axis=1) axis_only = axis_outlier & (~is_outlier) # ------------------------------------- # Recompute robust MCD on inliers only # ------------------------------------- X_inliers = X[~is_outlier] mcd_in, center_in, cov_in, _, _ = robust_mcd( X_inliers, support_fraction=MCD_SUPPORT, random_state=RANDOM_SEED ) # ---------------------------- # Diagnostic: Axis-only outliers check # ---------------------------- ax_idx = np.where(axis_only)[0] print(f"\n[INFO] Axis-only outliers (|z| > {AXIS_Z_THRESHOLD}) that are NOT Mahalanobis outliers: {len(ax_idx)}") if len(ax_idx) > 0: for i in ax_idx: print(f" - {species[i]:20s} | Z-scores: {robust_z[i]}") else: print(" - None found.") # ---------------------------- # Plotly traces # ---------------------------- traces = [] lof = LocalOutlierFactor(n_neighbors=min(LOF_NEIGHBORS, X.shape[0] - 1), contamination='auto') lof_labels = lof.fit_predict(X) lof_scores = lof.negative_outlier_factor_ # Print LOF outliers and scores print("\n📌 LOF Outliers (density-based):") lof_out = (lof_labels == -1) lof_idx = np.where(lof_out)[0] for i in lof_idx: print(f" {species[i]:<25} LOF score: {lof_scores[i]:.3f}") lof_out = (lof_labels == -1) lof_idx = np.where(lof_out)[0] if lof_idx.size > 0: traces.append(go.Scatter3d( x=X[lof_idx, 0], y=X[lof_idx, 1], z=X[lof_idx, 2], mode='markers', marker=dict(size=6, symbol=symbol3d('x'), line=dict(width=2), opacity=1.0, color='gray'), name='LOF outliers (density-based)' )) Ex, Ey, Ez = make_ellipsoid_surface(center_in, cov_in, alpha=ELLIPSOID_ALPHA, n_u=80, n_v=40) boot_centers = bootstrap_centers(X, n_boot=N_BOOT, random_state=RANDOM_SEED) if SHOW_BOOTSTRAP_CENTER else None # Inliers colored by D² traces.append(go.Scatter3d( x=X[inlier_idx, 0], y=X[inlier_idx, 1], z=X[inlier_idx, 2], mode='markers', marker=dict( size=[18 if species[i] == 'D. enhydrobia' else 14 for i in inlier_idx], color=[D2[i] for i in inlier_idx], colorscale='Viridis', cmin=vmin, cmax=vmax, colorbar=dict(title='Mahalanobis D²', thickness=20, len=0.9), opacity=0.85 ), name='Inliers (Mahalanobis)' )) # Outliers colored by D² traces.append(go.Scatter3d( x=X[out_idx, 0], y=X[out_idx, 1], z=X[out_idx, 2], mode='markers', marker=dict( size=10, symbol=symbol3d('diamond'), color=[D2[i] for i in out_idx], colorscale='Viridis', cmin=vmin, cmax=vmax, showscale=False, # Hide duplicate colorbar opacity=1.0 ), name=f'Outliers (D² > χ²₍{OUTLIER_ALPHA:.2f}₎)' )) # Axis outliers ax_idx = np.where(axis_only)[0] if ax_idx.size > 0: traces.append(go.Scatter3d( x=X[ax_idx, 0], y=X[ax_idx, 1], z=X[ax_idx, 2], mode='markers', marker=dict(size=10, symbol=symbol3d('diamond-open'), line=dict(width=2), color='gray', opacity=1.0), name=f'Axis outliers (|z|>{AXIS_Z_THRESHOLD} MAD)' )) z_base = float(np.min(X[:, 2]) - 1.0) drop_lines = [ go.Scatter3d( x=[X[i, 0], X[i, 0]], y=[X[i, 1], X[i, 1]], z=[X[i, 2], z_base], mode='lines', line=dict(color='rgba(128,128,128,0.8)', width=1), showlegend=False ) for i in range(X.shape[0]) ] # Main ellipsoid from robust MCD on full data Ex_full, Ey_full, Ez_full = make_ellipsoid_surface(center, cov, alpha=ELLIPSOID_ALPHA, n_u=80, n_v=40) ellipsoid_full_trace = go.Surface( x=Ex_full, y=Ey_full, z=Ez_full, opacity=0.15, colorscale='Blues', showscale=False, name=f'{int(ELLIPSOID_ALPHA*100)}% robust ellipsoid (MCD, full data)' ) # Inlier-only ellipsoid wireframe (OPTIONAL) # wireframe_traces = [] # # wireframe_traces = [] # # Vertical arcs ("latitude" lines) # for i in range(0, Ex.shape[0], 4): # wireframe_traces.append(go.Scatter3d( # x=Ex[i, :], # y=Ey[i, :], # z=Ez[i, :], # mode='lines', # line=dict(color='gray', width=2, dash='dot'), # showlegend=(i == 0), # name='Inlier-only ellipsoid' # )) # # Horizontal arcs ("longitude" lines) # for j in range(0, Ex.shape[1], 4): # wireframe_traces.append(go.Scatter3d( # x=Ex[:, j], # y=Ey[:, j], # z=Ez[:, j], # mode='lines', # line=dict(color='gray', width=2, dash='dot'), # showlegend=False # )) center_trace = go.Scatter3d( x=[center[0]], y=[center[1]], z=[center[2]], mode='markers', marker=dict(size=12, color='black', opacity=0.95), name='Robust center (MCD)' ) if boot_centers is not None: traces.append(go.Scatter3d( x=boot_centers[:, 0], y=boot_centers[:, 1], z=boot_centers[:, 2], mode='markers', marker=dict(size=3, opacity=0.15, color='black'), name=f'Bootstrap centers (n={N_BOOT})' )) traces = traces + [center_trace, ellipsoid_full_trace] + drop_lines layout = go.Layout( title=f'3D OR–GR–IR • {int(ELLIPSOID_ALPHA*100)}% Robust Mahalanobis Ellipsoid (MCD) + Outliers', scene=dict( xaxis=dict(title='OR genes'), yaxis=dict(title='GR genes'), zaxis=dict(title='IR genes'), aspectmode='cube', camera=dict( up=dict(x=0, y=0, z=1), center=dict(x=0, y=0, z=0), eye=dict(x=-2, y=-1.5, z=1.5) ) ), legend=dict(x=0.02, y=0.98) ) fig = go.Figure(data=traces, layout=layout) fig.show()