1
$\begingroup$

I was wondering if anyone could help with my little issue in manipulating tsv files in order to calculate estimated coverage across many tsv files.

I have tsv files that look something like this:

ID        Info    Info2  x
M80428    info    info   1000 
T83241    info    info   8400
E81292    info    info   10200

I also have another tsv file with all possible IDs and a corresponding number like so:

ID        y
M80428    27311
E11123    46531
P38163    14213
A34197    12892
T83241    37416
Q21345    27321
E81292    13101

x in the first tsv are estimated read counts and the y in the second is the genome length of each ID.

I would like to divide the "x" by the "y" in the second tsv that has the corresponding ID in order to get the estimated coverage of each hit. e.g. the first row would be 1000/27311, second row 8400/37416. I would then want this to be written to the fifth column of the tsv.

I am working in bash shell mostly here, and think that the way I would probably do this is with awk and by making the second tsv an associative array but I can't work out how exactly to do so. The python package dict seems like another option but my Python skills are practically zero.

Thanks for the help in advance!

$\endgroup$
  • 1
    $\begingroup$ I'm voting to close this question as off-topic because there is nothing related to bioinformatics in this post $\endgroup$ – Ram RS Feb 25 at 16:06
  • 1
    $\begingroup$ Hi @stanleyHo, could you provide more detail, what have you tried and what biological context is your problem? $\endgroup$ – Michael Feb 25 at 16:38
  • $\begingroup$ x in the first tsv are estimated read counts and the y in the second is the genome length of each ID. Trying to divide them to get estimated coverage. I have tried to set the second tsv as variables and then tried to awk out the columns with divisions but i'm stuck at how to link the ID and y together. Something like this: cat genome-size.txt | while read eachline; do echo $eachline; vname=echo $eachline| awk '{print $1}'; div=echo $eachline| awk '{print $2}' cat .tsv | awk -v div="$div" '{FS ="\t"} {OFS ="\t"} {print $1,$2,$3,$4,$4/div}' > new.tsv; done $\endgroup$ – Stanley Ho Feb 25 at 17:10
  • $\begingroup$ Hi @StanleyHo, note you can edit your questions to make them clearer. $\endgroup$ – Kamil S Jaron Feb 25 at 17:21
  • 1
    $\begingroup$ Use Python's pandas ... $\endgroup$ – Michael Feb 25 at 19:48
2
$\begingroup$

Python dictionaries are indeed a good tool for that.

Here is a script that takes 2 arguments (the file containing the counts and the file containing the genome length info), and prints the input data, with an added column containing coverage:

#!/usr/bin/env python3
import sys

counts_fname = sys.argv[1]
genome_lengths_fname = sys.argv[2]

genome_length_dict = {}
with open(genome_lengths_fname) as genome_lengths_f:
    # Checking the header line, for extra safety
    header = genome_lengths_f.readline().strip().split("\t")
    assert header[0] == "ID"
    assert header[1] == "y"
    # Reading the data lines
    for line in genome_lengths_f:
        [genome_id, genome_length] = line.strip().split("\t")
        genome_length_dict[genome_id] = int(genome_length)

with open(counts_fname) as counts_f:
    print("ID\tInfo\tInfo2\tx\tx/y")
    # Checking the header line, for extra safety
    header = counts_f.readline().strip().split("\t")
    assert header[0] == "ID"
    assert header[1] == "Info"
    assert header[2] == "Info2"
    assert header[3] == "x"
    # Reading the data lines
    for line in counts_f:
        [genome_id, info1, info2, counts] = line.strip().split("\t")
        # /!\ It would be different in python 2:
        # different print and division behaviours
        print(
            genome_id, info1, info2, int(counts),
            int(counts) / genome_length_dict[genome_id],
            sep="\t")

sys.exit(0)

Testing it:

$ ./coverage.py x.tsv y.tsv
ID  Info    Info2   x   x/y
M80428  info    info    1000    0.0366152832192157
T83241  info    info    8400    0.22450288646568314
E81292  info    info    10200   0.7785665216395695
| improve this answer | |
$\endgroup$
  • $\begingroup$ This worked perfectly thanks so much! $\endgroup$ – Stanley Ho Feb 26 at 11:55
1
$\begingroup$

I would do it in R:

tab1 <- read.table('file1.tsv', header = T, stringsAsFactors = F)
tab2 <- read.table('file2.tsv', header = T, stringsAsFactors = F)
merged_tav <- merge(reads_tab, genome_len_tab)
merged_tav$z <- merged_tav$x / merged_tav$y
merged_tav

returns

      ID Info Info2     x     y          z
1 E81292 info  info 10200 13101 0.77856652
2 M80428 info  info  1000 27311 0.03661528
3 T83241 info  info  8400 37416 0.22450289

note that only IDs that are in both tables will be preserved.

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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