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

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

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


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)

• 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() Feb 12 '18 at 16:39
• Actually a generator expression for pangenome would be best and most memory efficient 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")

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

• 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 Feb 16 '18 at 23:44
• 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).
– bli
Feb 19 '18 at 9:25