Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-12T12:45:14.613Z
Answers.LastActivityDate: 2018-02-12T12:45:14.613Z
Answers.UserDisplayName: terdon
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-12T16:39:27.907Z
Answers.Text: 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()`
Answers.UserDisplayName: Chris_Rands
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-12T16:48:18.863Z
Answers.Text: Actually a generator expression for pangenome would be best and most memory efficient
Answers.UserDisplayName: Chris_Rands
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-12T10:41:18.547Z
Answers.LastActivityDate: 2018-02-12T10:41:18.547Z
Answers.UserDisplayName: mgalardini
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-16T23:44:26.090Z
Answers.Text: 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
Answers.UserDisplayName: Anna1364
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-19T09:25:23.133Z
Answers.Text: 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)`.
Answers.UserDisplayName: bli
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-12T16:08:52.723Z
Answers.LastActivityDate: 2018-02-12T16:08:52.723Z
Answers.UserDisplayName: bli
Body: <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>
  
Tags: <python><blast><perl>
Title: assigning different characters to present and absent genes from blast output using perl or python
ThreadText: # 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`](https://docs.python.org/3.6/library/stdtypes.html#frozenset.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)`.  
  
  ---
  
  