Genotyping 1-2x Sample Resulting in nan

[Pkaps25]

Hello,

I am using danbing-tk to genotype a set of 1-2x bacterial samples. When running length prediction, my output is entirely nan for all loci in all samples. Is this level of coverage too low for genotyping SRS?

I made the following modifications to script/kmc2length.py to fix Numpy and Python errors:
line 14: q1s = np.nanquantile(data[m], 0.25)
line 15: q3s = np.nanquantile(data[m], 0.75)
line 54: nloci0 = df.shape[0]

The changes on lines 14 and 15 were made because of the following error:

Loading metadata and precomputed TR/NTR LSB
Computing NTR LSB in current dataset
../danbing-tk/script/kmc2length.py:30: RuntimeWarning: invalid value encountered in true_divide
normcovmat = covmat / (covmat@ctrlsize / np.sum(ctrlsize))[:,None]
Traceback (most recent call last):
File "../danbing-tk/script/kmc2length.py", line 127, in
pbamcov, _, ntrbiasmat = processCtrlBamCov(rawcovmat[:,2:].astype(float))
File "../danbing-tk/script/kmc2length.py", line 34, in processCtrlBamCov
badmask = np.logical_or(badmask, get1DIQRmask(stds))
File "../danbing-tk/script/kmc2length.py", line 14, in get1DIQRmask
q1s = np.quantile(data[m], 0.25)
File "<array_function internals>", line 5, in quantile
File "/home/pk445/miniconda3/envs/danbing/lib/python3.8/site-packages/numpy/lib/function_base.py", line 3844, in quantile
return _quantile_unchecked(
File "/home/pk445/miniconda3/envs/danbing/lib/python3.8/site-packages/numpy/lib/function_base.py", line 3851, in _quantile_unchecked
r, k = _ureduce(a, func=_quantile_ureduce_func, q=q, axis=axis, out=out,
File "/home/pk445/miniconda3/envs/danbing/lib/python3.8/site-packages/numpy/lib/function_base.py", line 3429, in _ureduce
r = func(a, **kwargs)
File "/home/pk445/miniconda3/envs/danbing/lib/python3.8/site-packages/numpy/lib/function_base.py", line 3967, in _quantile_ureduce_func
x1 = take(ap, indices_below, axis=axis) * weights_below
File "<array_function internals>", line 5, in take
File "/home/pk445/miniconda3/envs/danbing/lib/python3.8/site-packages/numpy/core/fromnumeric.py", line 191, in take
return _wrapfunc(a, 'take', indices, axis=axis, out=out, mode=mode)
File "/home/pk445/miniconda3/envs/danbing/lib/python3.8/site-packages/numpy/core/fromnumeric.py", line 58, in _wrapfunc
return bound(*args, **kwds)
IndexError: cannot do a non-empty take from an empty axes.

The change on line 54 was made in response to this error:

Loading metadata and precomputed TR/NTR LSB
Computing NTR LSB in current dataset
Traceback (most recent call last):
File "../danbing-tk/script/kmc2length.py", line 128, in
pbamcov, _, ntrbiasmat = processCtrlBamCov(rawcovmat[:,2:].astype(float))
File "../danbing-tk/script/kmc2length.py", line 26, in processCtrlBamCov
cov = covmat@ctrlsize / np.sum(ctrlsize)
ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 3460 is different from 3467)

[mchaisso]

We would recommend ~30X. At 1-2X, the counts of coverage will be too low to estimate length.

…

On Mon, Feb 1, 2021 at 2:08 PM Pkaps25 ***@***.***> wrote: Hello, I am using danbing-tk to genotype a set of 1-2x bacterial samples. When running length prediction, my output is entirely nan for all loci in all samples. Is this level of coverage too low for genotyping SRS? I made the following modifications to script/kmc2length.py to fix Numpy and Python errors: line 14: q1s = np.nanquantile(data[m], 0.25) line 15: q3s = np.nanquantile(data[m], 0.75) line 54: nloci0 = df.shape[0] The changes on lines 14 and 15 were made because of the following error: Loading metadata and precomputed TR/NTR LSB Computing NTR LSB in current dataset ../danbing-tk/script/kmc2length.py:30: RuntimeWarning: invalid value encountered in true_divide normcovmat = covmat / ***@***.*** / np.sum(ctrlsize))[:,None] Traceback (most recent call last): File "../danbing-tk/script/kmc2length.py", line 127, in pbamcov, _, ntrbiasmat = processCtrlBamCov(rawcovmat[:,2:].astype(float)) File "../danbing-tk/script/kmc2length.py", line 34, in processCtrlBamCov badmask = np.logical_or(badmask, get1DIQRmask(stds)) File "../danbing-tk/script/kmc2length.py", line 14, in get1DIQRmask q1s = np.quantile(data[m], 0.25) File "<*array_function* internals>", line 5, in quantile File "/home/pk445/miniconda3/envs/danbing/lib/python3.8/site-packages/numpy/lib/function_base.py", line 3844, in quantile return _quantile_unchecked( File "/home/pk445/miniconda3/envs/danbing/lib/python3.8/site-packages/numpy/lib/function_base.py", line 3851, in _quantile_unchecked r, k = _ureduce(a, func=_quantile_ureduce_func, q=q, axis=axis, out=out, File "/home/pk445/miniconda3/envs/danbing/lib/python3.8/site-packages/numpy/lib/function_base.py", line 3429, in _ureduce r = func(a, **kwargs) File "/home/pk445/miniconda3/envs/danbing/lib/python3.8/site-packages/numpy/lib/function_base.py", line 3967, in _quantile_ureduce_func x1 = take(ap, indices_below, axis=axis) * weights_below File "<*array_function* internals>", line 5, in take File "/home/pk445/miniconda3/envs/danbing/lib/python3.8/site-packages/numpy/core/fromnumeric.py", line 191, in take return _wrapfunc(a, 'take', indices, axis=axis, out=out, mode=mode) File "/home/pk445/miniconda3/envs/danbing/lib/python3.8/site-packages/numpy/core/fromnumeric.py", line 58, in _wrapfunc return bound(*args, **kwds) IndexError: cannot do a non-empty take from an empty axes. The change on line 54 was made in response to this error: Loading metadata and precomputed TR/NTR LSB Computing NTR LSB in current dataset Traceback (most recent call last): File "../danbing-tk/script/kmc2length.py", line 128, in pbamcov, _, ntrbiasmat = processCtrlBamCov(rawcovmat[:,2:].astype(float)) File "../danbing-tk/script/kmc2length.py", line 26, in processCtrlBamCov cov = ***@***.*** / np.sum(ctrlsize) ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 3460 is different from 3467) — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#8>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAKPGPSKPSBYMFT2DSRUVTTS44Q47ANCNFSM4W5RQ7CA> .

[joyeuxnoel8]

Hi @Pkaps25 ,
As I recall, you have a different set of VNTR loci of interest. The model would not be able to estimate lengths unless pre-trained with some ground truths (SRS + VNTR length in assembly) from other samples. The model works by correcting VNTR-specific LSB using prior observations. If there's no available data to compute VNTR-specific LSB, cross-sample comparison over kmer dosages is still valid as long as LSB analysis at non-repetitive regions indicate similar bias pattern.

[Pkaps25]

Hi, where do the ground truths come into play in LSB calculation? By ground truths I am assuming knowledge of true VNTR lengths in assemblies. I thought that to calculate LSB, all we need is the SRS data and the set of control regions. Please let me know if I am not understanding LSB calculation correctly.

[joyeuxnoel8]

So LSB for each locus is computed as (kmer_dosage / global_coverage) / VNTR_length_in_assembly. Although it might be possible to train a model to predict the LSB in VNTR regions based on the LSB in control regions so as to remove the dependency on ground truths, this hasn't been implemented yet. Hope this clarifies your question.

[Pkaps25]

Yes, thank you very much!

[Pkaps25]

Hi, I have one more question about finding LSBs for custom RPGG. I constructed my RPGG without pruning, so without matching SRS for assemblies. I have a set of ~100 samples I would like to genotype using this RPGG. Should I compute LSB for a subset of these samples and run length prediction on the remaining samples using these LSBs or compute LSB for all samples? Thank you.

Also, to output kmer dosages in estimated_TR_len.tsv for intra-sample comparison, is it sufficient to remove the bias[m] term in line 83 of script/kmc2length.py?

I know I've been asking many questions, so thank you for all of your time and help.

[joyeuxnoel8]

Just checking if I understand your dataset correctly, so you have a few TB assemblies but none of those match (coming from the same genome/strain) the ~100 SRS samples you have, and you don't have any LRS samples that match the SRS samples either. If this is the case, then the current implementation does not support computing LSB directly from SRS samples.

I probably wasn't clear enough on this, so let me rewrite how LSB is computed:
(kmer_dosage[i,j] / global_coverage[i]) / VNTR_length_in_assembly[i,j]
, where i and j denote sample i and locus j; kmer_dosage is obtained from danbing-tk align; VNTR_length_in_assembly is measured from the long-read assembly. You can see that kmer_dosage and VNTR_length_in_assembly have to be from the same sample.

Just a note, it might be helpful to visualize the LSBs at non-repetitive regions to better understand the bias pattern of your samples. As shown in LSB_analysis.ipynb.

Also, to output kmer dosages in estimated_TR_len.tsv for intra-sample comparison, is it sufficient to remove the bias[m] term in line 83 of script/kmc2length.py?

I don't think it's safe to remove the bias term for intra-sample comparison, or comparing across loci, unless you can show that the sequencing depth is almost constant regardless of GC content or kmer compositions. It is generally okay to ignore the bias term for inter-sample comparison as long as they are from the same sequencing run with good consistency.

I know I've been asking many questions, so thank you for all of your time and help.

Not a problem at all, glad to have your feedback! Hope your research goes well!