<p>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 <code>awk</code>:</p>

<pre><code>$ 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
</code></pre>

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

<pre><code>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 
</code></pre>

<p>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):</p>

<pre><code>$ grep -wof pangenome.txt blast.individual1.txt.orig 
Ha16_00021664
Ha8_00042762
Ha4_00025807
</code></pre>

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

<pre><code>$ sort &lt;(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
</code></pre>

<p>For multiple files, just run it on all of them:</p>

<pre><code>$ sort &lt;(grep -owf pangenome.txt blast.individual*.txt) pangenome.txt | 
    uniq -c | awk '{print $2,$1-1}'
</code></pre>

CC BY-SA 3.0
2018-02-12T12:45:14.613Z
2018-02-12T12:45:14.613Z
terdon
<p>Using <code>python</code> it should be relatively trivial, using <code>sets</code>. The following snippet works with your pangenome file and one blast file; adding more blast files should be trivial, making use of the <a href="https://docs.python.org/3.6/library/stdtypes.html#frozenset.union" rel="nofollow noreferrer"><code>union</code></a> method of <code>python</code>'s <code>sets</code>.</p>

<pre><code>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)
</code></pre>

CC BY-SA 3.0
2018-02-12T16:39:27.907Z
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()`
Chris_Rands
CC BY-SA 3.0
2018-02-12T16:48:18.863Z
Actually a generator expression for pangenome would be best and most memory efficient
Chris_Rands
CC BY-SA 3.0
2018-02-12T10:41:18.547Z
2018-02-12T10:41:18.547Z
mgalardini
<p>Here is a way to do it using a python dictionary:</p>

<pre><code>#!/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)
</code></pre>

<p>Executing it:</p>

<pre><code>$ ./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
</code></pre>

CC BY-SA 3.0
2018-02-16T23:44:26.090Z
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
Anna1364
CC BY-SA 3.0
2018-02-19T09:25:23.133Z
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
CC BY-SA 3.0
2018-02-12T16:08:52.723Z
2018-02-12T16:08:52.723Z
bli
<p>I have file called pangenome.txt which has only one column containing the names of the genes</p>

<pre><code>head pangenome.txt
    Ha1_00044161
    Ha11_00031248
    Ha16_00021664
    Ha8_00042762
    Ha4_00025807
</code></pre>

<p>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</p>

<pre><code>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    
</code></pre>

<p>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:</p>

<pre><code>head output.individual1.txt
    Ha1_00044161   0
    Ha11_00031248  0
    Ha16_00021664  1
    Ha8_00042762   1
    Ha4_00025807   1
</code></pre>

<p>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</p>

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

(Asked by: username_0 on 2018-02-12)

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

**Tags:** python, blast, perl

---

## 3 answers

### Answer by username_1 on 2018-02-12. Score: 3 | Accepted answer

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 by username_2 on 2018-02-12. Score: 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)

> **Comments:**
> **username_3**: 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()`  
> **username_3**: Actually a generator expression for pangenome would be best and most memory efficient  

---

### Answer by username_4 on 2018-02-12. Score: 2

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

> **Comments:**
> **username_0**: 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  
> **username_4**: 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)`.  

---

