assigning different characters to present and absent genes from blast output using perl or python

Question

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

head pangenome.txt
    Ha1_00044161
    Ha11_00031248
    Ha16_00021664
    Ha8_00042762
    Ha4_00025807

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

Answer 1

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){
            a[$2]=1
        }
       }END{
        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 
Ha16_00021664
Ha8_00042762
Ha4_00025807

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}'

Answer 2

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)
    else:
        print('%s\t0' % x)

Answer 3

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")

exit(0)

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