Why does sklearn's permutation_test_score return a different ROC AUC score than when I calculate manually using predict_proba and roc_auc_score?

I can't reproduce the ROC AUC score I get from permutation_test_score when I try to calculate the score manually using predict_proba and roc_auc_score. It matters because this could be the difference between sig./non-sig. results for a project.

Visual of the output:

• (In yellow) Score from permutation_test_score = 0.5256
• (In green) Score from roc_auc_score using predict_proba = 0.5416
• (In red) 97.5 percentile line representing p=.05 sig. threshold

Here is the code that produced that visual, which I adapted from an sklearn example in their permutation_test_score documentation. The grid.best_estimator_ object is a RandomForestClassifier that is the result of a randomized grid search with the same exact cross-validation you see below — I can include the code for that if helpful. Also, if it would help to have a self-contained reproducible example ripped straight from that sklearn example, I can provide that too (didn't for sake of space):

rskf = StratifiedKFold(n_splits=5)

n_permutations = 300
###### Use best estimator and run it on the Validation set. Validation set targets are permuted n times.
score_ofc, perm_scores_ofc, pvalue_ofc = permutation_test_score(
    grid.best_estimator_, 
    X_val, 
    y_val, 
    scoring="roc_auc", 
    cv=rskf, 
    n_permutations=n_permutations, 
    n_jobs=6, 
    random_state=42,
    verbose=1
)

###### manual calculation of roc_auc score
y_pred_val = grid.best_estimator_.predict_proba(X_val)[:,1]
roc_auc_val = roc_auc_score(y_val, y_pred_val)
p_val_man = (np.sum(perm_scores_ofc >= roc_auc_val) + 1.0) / (n_permutations + 1)


##### Plot permutations 
fig, ax = plt.subplots()
plt.figure(figsize = (5,5))

ax.hist(perm_scores_ofc, bins=20, density=True)
###### Compare roc_auc_val score to score_ofc 
ax.axvline(roc_auc_val, ls="--", color="g", lw=3)
ax.axvline(score_ofc, ls="--", color="y", lw=3)
###### Include line showing the p=.05 significance level
ax.axvline(np.percentile(perm_scores_ofc, 97.5), ls="-", color="r", lw=3)
ax.set_xlabel("ROC AUC score")
_ = ax.set_ylabel("Probability")     # copy-pasted all this code from the sklearn documentation, and I'm not sure why they called this "probability"

print('''
Green = Score on original data using "manual" predict_proba method
      = {}
p-val = {}

Yellow = Score on original data using "automatic" grid.score_ method
      = {}
p-val = {}

97.5 Percentile value: {}
'''.format(roc_auc_val, p_val_man, score_ofc, pvalue_ofc, np.percentile(perm_scores_ofc, 97.5)))

plt.show()

I've seen one or two other related questions (ex: here) that had something to do with the difference between a scorer using decision_function vs. predict_proba, but that shouldn't be an issue here since RandomForestClassifier has no decision_function attribute. So permutation_test_score has to be using predict_proba, right? But then why am I getting different results?

Thank you for any and all help! I've been trying to figure this out for days.

Edit

Including my code for the original pipeline and grid search, for completeness.

pca = PCA()
pipe = Pipeline(
    [
        ('scaler', MaxAbsScaler()),
        ('pca', pca),
        ('classifier', RandomForestClassifier()),
    ]
)

param_grid = [
    {
        'classifier': [RandomForestClassifier(random_state=42, n_jobs=-1)],
        'classifier__max_depth' : [i for i in range(1, 8, 2)],
        'scaler': [RobustScaler()],
        'pca__n_components': [33],
        'classifier__n_estimators' : [250],
        'classifier__criterion' : ['gini'],
        'classifier__max_features' : [0.3],
        'classifier__min_samples_split': [12],
        'classifier__min_samples_leaf': [9]
    }
]
###################################################
### USE GRID SEARCH TO FIND BEST HYPERPARAMETERS ###
# SCORING = ROC AUC

rskf = StratifiedKFold(n_splits=5)

grid = RandomizedSearchCV(pipe, 
                          param_grid,
                          n_iter=60,                           # Seemed like the right balance between computation time and exhaustiveness
                          random_state=42,
                          scoring='roc_auc',
                          cv=rskf,
                          refit=True,
                          return_train_score=True,
                          verbose=1,
                          n_jobs=6
                         ).fit(X_train, y_train)               # Fit all iterations on training data