Error summary: IllegalStateException: ARPACK returns non-zero info = 1 Maximum number of iterations taken. (Refer ARPACK user guide for details)
after calling hwe_normalized_pca(mt.GT, k=2) and I have trouble pinpointing the cause of the problem. What’s worse, when I’m calling this function several times in a row on the same input, sometimes it errors out and sometimes it returns without an error. Is that even possible with a proper installation? I have now made myself a jupyter notebook with the following single cell just to prove to myself that I’m not going mad:
hl.init(quiet=True)
hl.set_global_seed(42)
mt = hl.read_matrix_table(my_path)
for i in range(12):
try:
hl.hwe_normalized_pca(
mt.GT,
k = 2
)
print(f"{i} success")
except:
print(f"{i} error")
where mt is a product of a several filtering steps performed on a genuine vcf.
Also, when I do
hl.pca(
hl.int(hl.is_defined(mt.GT)),
k = 2
)
on my mt, I always immediately get another kind of error: Error summary: HailException: Cannot create RowMatrix: filtered entry at row 8111 and col 9,
but I don’t know if it is related to the problem with hwe_normalizd_pca.
Any suggestions would be appreciated:
What qualities of any mt could in general cause the hwe_normalized_pca to not converge?
Should I suspect that there’s something wrong with my installation? I’m running hail locally, on a laptop.
Shape of mt is (18205, 39), its only non-key field is GT, and the col-keys are generic strings, so I’d be happy to share it.
You could make normal PCA work with hl.or_else to set a value for missing entries.
There’s almost certainly no issue with your installation. I can only speculate as to the issues, but I’d guess that your matrix is exhibiting some kind of numerical instability that makes ARPACK’s SVD routines occasionally not converge. Your dataset is quite small. You might try collecting it locally and converting it to a numpy array for use with sklearn.decomposition.PCA. If that succeeds, it’s possible Hail could use whatever parameters sklearn is using to call ARPACK and achieve more consistent results.
Thanks, @danking, for the response!
It’s helpful to hear that reinstalling probably won’t help.
Yup, the dataset is small, partly because I was trying to take a subset of it which still throws the same error as the original
Not sure where I would put hl.or_else in hl.int(hl.is_defined(mt.GT)). But, swapping mt.GT for mt.unfilter_entries().GT made the hl.pca return without an error. Similarly
from sklearn.decomposition import PCA
mt_2 = mt.unfilter_entries()
df = mt_2.select_entries(
n = hl.int(hl.is_defined(mt_2.GT))
).key_rows_by().select_rows().make_table().to_pandas()
PCA(2).fit_transform(df)
returns without an error, but it would complain about NaN without the unfilter_entries (same with say n = hl.or_else(mt_2.GT.n_alt_alleles(), 0)).
Still, unfilter_entries does not help with the hwe_normalized_pca.
Perhaps I’m having a problem with variants that are far from hw equilibrium. E.g. if I filter out low hw p-value before calling hwe_normalized_pca in my jupyter cell:
It sounds to me like hardy-weinberg normalization is producing a matrix that is likely to be poorly conditioned causing ARPACK to not converge.
Let’s compare apples to apples:
import hail as hl
mt = hl.read_mtarix_table('debugging/wierdness.mt')
mt_for_comparison = hl.methods.pca.hwe_normalize(mt.GT)
# ... run mt_for_comparison through sklearn PCA and hl.pca
This is not documented because it’s an internal method. This conversation may end up being an impetus for us to document it and expose it. The usage is similar to what Dan did above, but you’d need to import it the way you did. So:
from hail.methods.pca import hwe_normalize
norm_gt = hwe_normalize(mt.GT)
# note it's not a matrix table, its an entry expression. You can export it for sklearn though using:
norm_gt.export('norm_gt.tsv')
