I have file called pangenome.txt which has only one column containing the names of the genes

head pangenome.txt

I have another bunch of files in a tabular blast format (330 in total: individual1.blast.txt,individual2.blast.txt, ...,individual330.blast.txt ) which are the outputs of Diamond protein aligner. Below is part of individual1.blast.txt file

head blast.individual1.txt
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha16_00021664   77.4    31  7   0   94  2   49  79  7.3e-159.3
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha8_00042762    74.2    31  8   0   94  2   1   31  9.6e-158.9
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha4_00025807    77.4    31  7   0   94  2   1   31  1.3e-058.5
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha16_00022939   77.4    31  7   0   94  2   295 325 2.1e-057.8
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha12_00032339   71.0    31  9   0   94  2   2   32  2.8e-057.4
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha16_00021553   71.0    

I want to know which genes from the pangemome.txt is present in my individuals file and which are absent. for the example above this is my desired output:

head output.individual1.txt
    Ha1_00044161   0
    Ha11_00031248  0
    Ha16_00021664  1
    Ha8_00042762   1
    Ha4_00025807   1

My desired out put is a file with 2 columns (tabular format), the first column the list of the all genes form pangemome.txt pan, and the second column is the presence and absence of these genes in the individuals files if the gene is present present in the blast.individual1.txt, I want to assign 1 and if it is absent I want to assign 0. Can anyone help me to do this? Thanks a lot


I wouldn't do this in Python, myself. This is a very simple text parsing problem and the standard *nix tools will be able to do it very easily. For example with awk:

$ awk 'NR==FNR{a[$1]=0; next}
        if($2 in a){
        for(i in a){
            print i,a[i]
       }' pangenome.txt blast.individual1.txt 
Ha1_00044161 0
Ha4_00025807 1
Ha11_00031248 0
Ha16_00021664 1
Ha8_00042762 1

The above is just a one-liner, I expanded it for clarity. You can run it on as many files as you like:

awk 'NR==FNR{a[$1]=0; next}{if($2 in a){a[$2]=1}}END{for(i in a){print i,a[i]}}' pangenome.txt blast.individual*.txt 

For many large files, the fastest approach would probably be to just grep the names (unlike the awk, this will also match if the target gene is in the 1st field, but that's probably not an issue for you):

$ grep -wof pangenome.txt blast.individual1.txt.orig 

That will give you the list of genes found in the target blast files. To get the 1 and 0, you can do:

$ sort <(grep -owf pangenome.txt blast.individual1.txt) pangenome.txt | 
    uniq -c | awk '{print $2,$1-1}'
Ha1_00044161 0
Ha11_00031248 0
Ha16_00021664 1
Ha4_00025807 1
Ha8_00042762 1

For multiple files, just run it on all of them:

$ sort <(grep -owf pangenome.txt blast.individual*.txt) pangenome.txt | 
    uniq -c | awk '{print $2,$1-1}'

Using python it should be relatively trivial, using sets. The following snippet works with your pangenome file and one blast file; adding more blast files should be trivial, making use of the union method of python's sets.

pangenome = {x.rstrip() for x in open('pangenome.txt')}
blast = {x.rstrip().split()[1] for x in open('blast.txt')}
for x in pangenome:
    if x in blast:
        print('%s\t1' % x)
        print('%s\t0' % x)
  • 1
    $\begingroup$ Good approach but I would use a list comprehension [] for pangenome to preserve the original order. Also you'd be safer to use the with statement to ensure the file handles are explicitly closed. There is no need to use r.strip() when defining blast. The latest python string formatting uses .format() $\endgroup$ – Chris_Rands Feb 12 '18 at 16:39
  • 1
    $\begingroup$ Actually a generator expression for pangenome would be best and most memory efficient $\endgroup$ – Chris_Rands Feb 12 '18 at 16:48

Here is a way to do it using a python dictionary:

#!/usr/bin/env python3

import sys

gene_filename = sys.argv[1]
blast_filename = sys.argv[2]

# Set to 0 all genes in the gene file
genes = {}
with open(gene_filename, "r") as gene_file:
    for line in gene_file:
        genes[line.strip()] = 0

# Set to 1 those genes found in the blast file
with open(blast_filename, "r") as blast_file:
    for line in blast_file:
        gene = line.strip().split()[1]
        if gene in genes:
            genes[gene] = 1

# Print the results
for (gene, presence) in genes.items():
    print(gene, presence, sep="\t")


Executing it:

$ ./find_genes_in_blast.py pangenome.txt blast.individual1.txt 
Ha1_00044161    0
Ha11_00031248   0
Ha16_00021664   1
Ha8_00042762    1
Ha4_00025807    1
  • $\begingroup$ thanks for the code. I tried it, but I get an error message: File "./find_genes_in_blast.py", line 23 print(gene, presence, sep="\t") ^ SyntaxError: invalid syntax $\endgroup$ – Anna1364 Feb 16 '18 at 23:44
  • $\begingroup$ If you are using python 2, you need to put from __future__ import print_function at the beginning of the code. Or you need to use a correct syntax for print in python 2, for instance print "%s\t%d" % (gene, presence). $\endgroup$ – bli Feb 19 '18 at 9:25

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.