Trying to divide number in fourth TSV column by a certain number according to its corresponding identifier in the first column

Question

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!

Answer 1

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

Answer 2

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.