1
$\begingroup$

This is an add-on to the previous question I asked here.

I was able to deploy a "chunking" strategy to my kmer matrix and was able to get a file with the distance measurements! From what I could tell by using the head command was that the data looks like what it did before when I was using this approach on a smaller dataset so now I'm going to use this data to generate a tree using MEGAX.

Normally when I'd convert the tsv file to a meg file I'd manually add in the meg headers/requirements for conversion. However, I can't open the distance tsv file (its ~23 GB) to make these changes.

Is there an easy solution to this problem? How do you convert a tsv file to a meg file without opening the file?

EDIT:

For more context here is the format for each file.

Currently, my tsv file looks like this:

0 1 2 3 4 5 6 7 8 9 10 ... 178
1 0 # # # # # # # # #  ... #
2 # 0 # # # # # # # #  ... #
3 # # 0 # # # # # # #  ... #
4 # # # 0 # # # # # #  ... #
5 # # # # 0 # # # # #  ... #
6 # # # # # 0 # # # #  ... #
7 # # # # # # 0 # # #  ... #
8 # # # # # # # 0 # #  ... #
9 # # # # # # # # 0 #  ... #
10 # # # # # # # # # 0 ... #
...
178 # # # # # # # # # #  ... 0

This is a jaccard distance matrix of 178 samples. The # represent the distance values for each sample compared to the other. I was able to calcualte the distance using pandas and scipy: matrix = pd.DataFrame(squareform(pdist(data, 'jaccard')))

And I would like it to look like this:

#mega
TITLE: Jaccard distance of 178 samples

sample_1
sample_2
sample_3
...
sample_178

distance_1
distance_1 distance_2
distance_1 distance_2 distance_3
distance_1 distance_2 distance_3 distance_4
...
distance_1 distance_2 distance_3 distance_4...distance_178
$\endgroup$
2
  • 1
    $\begingroup$ I added some extra information to give more context to this issue! Let me know if you'd like me to clarify something further or if you'd like more information. I appreciate the help! $\endgroup$
    – rimo
    Jan 12 at 20:47
  • $\begingroup$ Thanks @rimo .. okay this ain't the data set I imagined. I'll need to reproduce it. Are you using SciPy, i.e. where is the pdist() method coming from? $\endgroup$
    – M__
    Jan 12 at 23:41

1 Answer 1

1
$\begingroup$

Note, its a while since I've used scipy and I don't fancy loading an deleting for this test ... so in the absence of pdist() the solution would be ...

import numpy as np
import pandas as pd
from pathlib import Path

data = [[1,2,3],[0,5,6],[7,0,9],[10,11,0]]
arr = np.tril([data]).astype('float')
arr[arr == 0] = 'nan' 
arr = np.where(np.isnan(arr), '', arr)
columns = ['a','b','c', 'd']
file = Path('/Volumes/Data2/test/test1.meg')
if file.is_file():
    file.unlink()
with open(file, 'a') as fout:
    fout.write("#MEGA\n")
    for col in columns:
        fout.write(col + "\n")
    fout.write("\n0\n")
    np.savetxt(fout, arr.reshape(3,-1), fmt='%s')

output

#MEGA
a
b
c
d

0
1.0 0.0
5.0 7.0 0.0
9.0 10.0 11.0 0.0

That works. Keep in mind where it says arr.reshape(3,-1) the 3 should be the number of columns -1 in your matrix: in your case this appears to be 177 (178 - 1).

The gist is to use np.tril, which doesn't exists in pandas. Hence don't use pandas.

At a guess with scipy then it would be ...

import numpy as np
import pandas as pd
from pathlib import Path
from scipy import pdist

tmp = [[1,2,3],[4,5,6],[7,8,9],[10,11,12]]

data = pdist(tmp, 'jaccard')
arr = np.tril([data]).astype('float')
arr[arr == 0] = 'nan' 
arr = np.where(np.isnan(arr), '', arr)
columns = ['a','b','c', 'd']
file =Path('/Volumes/Data2/test/test1.meg')
if file.is_file():
    file.unlink()
with open(file, 'a') as fout:
    fout.write("#MEGA\n")
    for col in columns:
        fout.write(col + "\n")
    fout.write("\n0\n")
    np.savetxt(fout, arr.reshape(3,-1), fmt='%s')

You might need a quick perl -p -i -e 's/[" "]+/ /g' test1.meg or even perl -p -i -e 's/[" "]+/\t/g' test1.meg to set the file for MEGA. I dunno why but \s or " " in the transform part of perl pie isn't working for me.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.