bgzip, tabix and sorting

[carolinehey]

Hi @wdecoster

I have 2 nanopolish output tsv files (calls and frequencies). I would like to use them as input for methplotlib. I have tried to sort as you wrote in an example but when i want to use tabix it gives an unsorted error.

cat <(head -n1 sample8_methylation_calls_2.tsv) <(tail -n +2 sample8_methylation_calls_2.tsv | sort -k2,2 -k3,3) | bgzip > sample8_sort_calls.tsv.gz

tabix -S1 -s1 -b3 -e4 sample8_methylation_calls_2.tsv.gz

Looking at the sorted and decompressed tsv file i see that the problem is between the plus and minus strand:

Should it only be sorted by the 3 column and how can this be done?

[carolinehey]

When i sort like this and use tabix afterwards the file does not work in methplotlib
cat <(head -n1 sample8_methylation_calls_2.tsv) <(tail -n +2 sample8_methylation_calls_2.tsv | sort -k 3) | bgzip > sample8_sort_calls.tsv.gz

(methplotlib_caroline_python3) root@srvodeappgsq02:/longread/G128-2021-Imprinting-disorders/20220315_1228_1G_PAI63326_dd13a3c5# methplotlib -m sample8_sort2_calls.tsv.gz -n calls -w chr15:24,500,000-25,500,000 -f /longread/resources/reference_genomes/GRCh38_Hg38/hg38.fa
Traceback (most recent call last):
File "/root/miniconda3/envs/methplotlib_caroline_python3/lib/python3.10/site-packages/pandas/core/indexes/base.py", line 3621, in get_loc
return self._engine.get_loc(casted_key)
File "pandas/_libs/index.pyx", line 136, in pandas._libs.index.IndexEngine.get_loc
File "pandas/_libs/index.pyx", line 142, in pandas._libs.index.IndexEngine.get_loc
TypeError: '85285 False
85303 True
85322 True
85341 True
85358 True
...
13851 True
13873 True
13895 False
13917 False
13939 True
Name: log_lik_ratio, Length: 166259, dtype: bool' is an invalid key

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "/root/miniconda3/envs/methplotlib_caroline_python3/lib/python3.10/site-packages/pandas/core/series.py", line 1086, in setitem
self._set_with_engine(key, value)
File "/root/miniconda3/envs/methplotlib_caroline_python3/lib/python3.10/site-packages/pandas/core/series.py", line 1147, in _set_with_engine
loc = self.index.get_loc(key)
File "/root/miniconda3/envs/methplotlib_caroline_python3/lib/python3.10/site-packages/pandas/core/indexes/base.py", line 3628, in get_loc
self._check_indexing_error(key)
File "/root/miniconda3/envs/methplotlib_caroline_python3/lib/python3.10/site-packages/pandas/core/indexes/base.py", line 5637, in _check_indexing_error
raise InvalidIndexError(key)
pandas.errors.InvalidIndexError: 85285 False
85303 True
85322 True
85341 True
85358 True
...
13851 True
13873 True
13895 False
13917 False
13939 True
Name: log_lik_ratio, Length: 166259, dtype: bool

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "/root/miniconda3/envs/methplotlib_caroline_python3/bin/methplotlib", line 8, in
sys.exit(main())
File "/root/miniconda3/envs/methplotlib_caroline_python3/lib/python3.10/site-packages/methplotlib/methplotlib.py", line 26, in main
meth_browser(meth_data=meth_data,
File "/root/miniconda3/envs/methplotlib_caroline_python3/lib/python3.10/site-packages/methplotlib/methplotlib.py", line 53, in meth_browser
meth_traces = plots.methylation(meth_data, dotsize=dotsize, binary=binary, minqual=minqual)
File "/root/miniconda3/envs/methplotlib_caroline_python3/lib/python3.10/site-packages/methplotlib/plots.py", line 106, in methylation
make_per_read_meth_traces_llr(table=meth.table,
File "/root/miniconda3/envs/methplotlib_caroline_python3/lib/python3.10/site-packages/methplotlib/plots.py", line 183, in make_per_read_meth_traces_llr
table.loc[:, "llr_scaled"] = rescale_log_likelihood_ratio(table["log_lik_ratio"].copy())
File "/root/miniconda3/envs/methplotlib_caroline_python3/lib/python3.10/site-packages/methplotlib/plots.py", line 269, in rescale_log_likelihood_ratio
llr[llr > 0] = scaler.fit_transform(llr[llr > 0].values.reshape(-1, 1)).tolist()
File "/root/miniconda3/envs/methplotlib_caroline_python3/lib/python3.10/site-packages/pandas/core/series.py", line 1128, in setitem
self._set_values(indexer, value)
File "/root/miniconda3/envs/methplotlib_caroline_python3/lib/python3.10/site-packages/pandas/core/series.py", line 1192, in _set_values
self._mgr = self._mgr.setitem(indexer=key, value=value)
File "/root/miniconda3/envs/methplotlib_caroline_python3/lib/python3.10/site-packages/pandas/core/internals/managers.py", line 337, in setitem
return self.apply("setitem", indexer=indexer, value=value)
File "/root/miniconda3/envs/methplotlib_caroline_python3/lib/python3.10/site-packages/pandas/core/internals/managers.py", line 304, in apply
applied = getattr(b, f)(**kwargs)
File "/root/miniconda3/envs/methplotlib_caroline_python3/lib/python3.10/site-packages/pandas/core/internals/blocks.py", line 943, in setitem
values[indexer] = value
ValueError: shape mismatch: value array of shape (96245,1) could not be broadcast to indexing result of shape (96245,)

[wdecoster]

For the calls file you should sort on the first (chromosome), third (start), and fourth (end) columns.

[carolinehey]

Sorry for my ignorance, but how is this done?

And should it also be done on the frequency file?

[wdecoster]

I think

cat <(head -n1 sample8_methylation_calls_2.tsv) <(tail -n +2 sample8_methylation_calls_2.tsv | sort -k1,1n -k3,3 -k4,4) | bgzip > sample8_sort_calls.tsv.gz

The frequency file has chromosome, start and end in columns 1,2,3 so that would probably look like

cat <(head -n1 sample8_methylation_freqs.tsv) <(tail -n +2 sample8_methylation_freqs.tsv | sort -k1,1n -k2,2 -k3,3) | bgzip > sample8_sort_freqs.tsv.gz

[wdecoster]

Please create a separate issue for this and also provide the version of methplotlib that you are using.