<p>The following python script takes a population name as argument, searches for all files for that population (assuming some regularity in the file names), and records the counts of "0" and "1" for all the genes encountered:</p>

<pre class="lang-python prettyprint-override"><code>#!/usr/bin/env python3

from sys import argv
from glob import glob
from collections import defaultdict

filenames = glob("HI*Index_*_{pop}-*_all_PA".format(pop=argv[1]))

gene_counts = defaultdict(lambda : defaultdict(int))

for filename in filenames:
    with open(filename, "r") as genes_file:
        for line in genes_file:
            (gene, num) = line.strip().split()
            gene_counts[gene][num] += 1

for (gene, counts) in gene_counts.items():
    print(gene, counts["1"], counts["0"], sep="\t")
</code></pre>

<p>With the two example files you give being present in the current directory:</p>

<pre><code>$ ./count_gene.py MK4
Ha8_00040788    2   0
Ha4_00024045    2   0
Ha4_00025366    2   0
Ha16_00022130   1   1
Ha16_00023451   1   1
Ha8_00040789    2   0
Ha4_00025367    2   0
Ha4_00024046    0   2
Ha16_00022131   2   0
Ha16_00023452   2   0
</code></pre>

<p>I haven't tested with python 2. If you use python 2, you should at least add <code>from __future__ import print_function</code> as first import, or adapt the <code>print</code> part.</p>

<p>Then, to run on all populations:</p>

<pre><code>for pop in $(cat pop.txt)
do
    ./count_gene.py ${pop} &gt; ${pop}
done
</code></pre>

CC BY-SA 3.0
2018-04-04T17:10:01.943Z
2018-04-04T17:10:01.943Z
bli
<p>I can show you how I would approach it in <code>R</code>. You'll have to adjust and customize it to your own data of course. I only had two files to play with, from your example. Especially for the <code>merge</code> step you should build a <code>for</code> loop around it so that it will merge every next file automatically, since you have 500 files.</p>

<p>Here an example for the KS2W files:</p>

<pre><code>myFiles &lt;- dir(pattern="KS2W")

myFiles
[1] "HI.3577.005.Index_25.QPG4_KS2W-35_all_PA.txt"
[2] "HI.3577.005.Index_9.QPA2_KS2W-8_all_PA.txt"  

fileNames &lt;- list()

for(i in 1:2){
    fileNames[[i]] &lt;- read.table(myFiles[i], sep = " ")
}  

merged &lt;- merge(fileNames[[1]], fileNames[[2]], by.x = "V1", by.y = "V1", all = F, sort = F)

present &lt;- rowSums(merged[,2:3])

zeros &lt;- rowSums(merged[,2:3] == 0)

output &lt;- cbind(present, zeros)

rownames(output) &lt;- merged[,1]

output
              present zeros
Ha8_00040788        2     0
Ha4_00024045        2     0
Ha4_00025366        2     0
Ha16_00022130       1     1
Ha16_00023451       1     1
Ha8_00040789        2     0
Ha4_00025367        2     0
Ha4_00024046        0     2
Ha16_00022131       2     0
Ha16_00023452       2     0
</code></pre>

<p>To add more files with <code>merge</code>, you could add something like this:</p>

<pre><code>merged &lt;- merge(fileNames[[1]], fileNames[[2]], by.x = "V1", by.y = "V1", all = F, sort = F)

for(i in 3:30){
  merged &lt;- merge(merged, fileNames[[i]], by.x = "V1", by.y = "V1", all = F, sort = F)
}
</code></pre>

<p>You'll have to adjust all the numbers and coordinates to fit your data.</p>

CC BY-SA 3.0
2018-04-03T14:37:46.340Z
2018-04-03T15:37:16.977Z
2018-04-03T15:37:16.977Z
benn
<p>I have a folder containing 500 files from different individuals:</p>

<pre class="lang-bash prettyprint-override"><code>$ ls
HI.3577.005.Index_25.QPG4_KS2W-35_all_PA
HI.3577.005.Index_27.QPH4_SD1W-13_all_PA
HI.3577.005.Index_9.QPA2_KS2W-8_all_PA
HI.3577.007.Index_1.Q004B_SD1A-1_all_PA
HI.3577.007.Index_21.Q076C_SD2W-23_all_PA
HI.3577.007.Index_23.Q079_KS1A-15_all_PA
</code></pre>

<p>Each of these individuals belongs to a population (in the example above individuals belong to KS2W, SD1W, KS2W, SD1A, SD2W, and KS1A respectively). Each of my files looks as follows (this sample belongs to from "MK4" population):</p>

<pre class="lang-bash prettyprint-override"><code>$ head HI.3579.001.Index_9.Q089_MK4-4R_all_PA
Ha8_00040788 1
Ha4_00024045 1
Ha4_00025366 1
Ha16_00022130 0
Ha16_00023451 0
Ha8_00040789 1
Ha4_00025367 1
Ha4_00024046 0
Ha16_00022131 1
Ha16_00023452 1
</code></pre>

<p>or another sample from MK4 population</p>

<pre class="lang-bash prettyprint-override"><code>$ head HI.221.009.Index_9.P089_MK4-4R_all_PA
Ha8_00040788 1
Ha4_00024045 1
Ha4_00025366 1
Ha16_00022130 1
Ha16_00023451 1
Ha8_00040789 1
Ha4_00025367 1
Ha4_00024046 0
Ha16_00022131 1
Ha16_00023452 1
</code></pre>

<p>the first column is genes and the second columns refers to the presence of that gene (coded as 1) or absence of gene (coded as 0).</p>

<p>For all individuals belonging to each of the populations, I would like to count the total number individuals that have the gene and individuals that do not have the gene.</p>

<p>For example, for the 2 examples above, as they both belong to MK4 population my desired output is going to looke like this:</p>

<pre class="lang-bash prettyprint-override"><code>$ head MK4
Ha8_00040788   2  0
Ha4_00024045   2  0
Ha4_00025366   2  0
Ha16_00022130  1  1 
Ha16_00023451  1  1
Ha8_00040789   2  0
Ha4_00025367   2  0
Ha4_00024046   0  2
</code></pre>

<p>The first column is the gene, the second column is the number of individuals in the MK4 population that have the gene and the third column is the number of individuals that do not have that gene.</p>

<p>I have made a file containing population IDs</p>

<pre class="lang-bash prettyprint-override"><code>$ head pop.txt
KS2W
SD1W
SD1A
KS1A
</code></pre>

<p>Can anyone suggest me how I could do this? thanks</p>

CC BY-SA 3.0
2018-04-02T20:10:59.853Z
It's unlikely that *genes* are absent in individuals; do you mean SNPs?
gringer
CC BY-SA 3.0
2018-04-03T13:50:02.740Z
Why not use R??
benn
CC BY-SA 3.0
2018-04-04T16:31:42.323Z
@b.nota Some people have programming-language preferences.
bli
<python><perl>
counting individuals with the present gene and absent gene for each population
# counting individuals with the present gene and absent gene for each population

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

I have a folder containing 500 files from different individuals:
$ ls
HI.3577.005.Index_25.QPG4_KS2W-35_all_PA
HI.3577.005.Index_27.QPH4_SD1W-13_all_PA
HI.3577.005.Index_9.QPA2_KS2W-8_all_PA
HI.3577.007.Index_1.Q004B_SD1A-1_all_PA
HI.3577.007.Index_21.Q076C_SD2W-23_all_PA
HI.3577.007.Index_23.Q079_KS1A-15_all_PA

Each of these individuals belongs to a population (in the example above individuals belong to KS2W, SD1W, KS2W, SD1A, SD2W, and KS1A respectively). Each of my files looks as follows (this sample belongs to from "MK4" population):
$ head HI.3579.001.Index_9.Q089_MK4-4R_all_PA
Ha8_00040788 1
Ha4_00024045 1
Ha4_00025366 1
Ha16_00022130 0
Ha16_00023451 0
Ha8_00040789 1
Ha4_00025367 1
Ha4_00024046 0
Ha16_00022131 1
Ha16_00023452 1

or another sample from MK4 population
$ head HI.221.009.Index_9.P089_MK4-4R_all_PA
Ha8_00040788 1
Ha4_00024045 1
Ha4_00025366 1
Ha16_00022130 1
Ha16_00023451 1
Ha8_00040789 1
Ha4_00025367 1
Ha4_00024046 0
Ha16_00022131 1
Ha16_00023452 1

the first column is genes and the second columns refers to the presence of that gene (coded as 1) or absence of gene (coded as 0).
For all individuals belonging to each of the populations, I would like to count the total number individuals that have the gene and individuals that do not have the gene.
For example, for the 2 examples above, as they both belong to MK4 population my desired output is going to looke like this:
$ head MK4
Ha8_00040788   2  0
Ha4_00024045   2  0
Ha4_00025366   2  0
Ha16_00022130  1  1 
Ha16_00023451  1  1
Ha8_00040789   2  0
Ha4_00025367   2  0
Ha4_00024046   0  2

The first column is the gene, the second column is the number of individuals in the MK4 population that have the gene and the third column is the number of individuals that do not have that gene.
I have made a file containing population IDs
$ head pop.txt
KS2W
SD1W
SD1A
KS1A

Can anyone suggest me how I could do this? thanks

**Tags:** python, perl
> **Comments:**
> **username_1**: It's unlikely that *genes* are absent in individuals; do you mean SNPs?  
> **username_2**: Why not use R??  
> **username_3**: @b.nota Some people have programming-language preferences.  

---

## 2 answers

### Answer by username_3 on 2018-04-04. Score: 1

The following python script takes a population name as argument, searches for all files for that population (assuming some regularity in the file names), and records the counts of "0" and "1" for all the genes encountered:
#!/usr/bin/env python3

from sys import argv
from glob import glob
from collections import defaultdict

filenames = glob("HI*Index_*_{pop}-*_all_PA".format(pop=argv[1]))

gene_counts = defaultdict(lambda : defaultdict(int))

for filename in filenames:
    with open(filename, "r") as genes_file:
        for line in genes_file:
            (gene, num) = line.strip().split()
            gene_counts[gene][num] += 1

for (gene, counts) in gene_counts.items():
    print(gene, counts["1"], counts["0"], sep="\t")

With the two example files you give being present in the current directory:
$ ./count_gene.py MK4
Ha8_00040788    2   0
Ha4_00024045    2   0
Ha4_00025366    2   0
Ha16_00022130   1   1
Ha16_00023451   1   1
Ha8_00040789    2   0
Ha4_00025367    2   0
Ha4_00024046    0   2
Ha16_00022131   2   0
Ha16_00023452   2   0

I haven't tested with python 2. If you use python 2, you should at least add from __future__ import print_function as first import, or adapt the print part.
Then, to run on all populations:
for pop in $(cat pop.txt)
do
    ./count_gene.py ${pop} > ${pop}
done


---

### Answer by username_2 on 2018-04-03. Score: 1

I can show you how I would approach it in R. You'll have to adjust and customize it to your own data of course. I only had two files to play with, from your example. Especially for the merge step you should build a for loop around it so that it will merge every next file automatically, since you have 500 files.
Here an example for the KS2W files:
myFiles <- dir(pattern="KS2W")

myFiles
[1] "HI.3577.005.Index_25.QPG4_KS2W-35_all_PA.txt"
[2] "HI.3577.005.Index_9.QPA2_KS2W-8_all_PA.txt"  

fileNames <- list()

for(i in 1:2){
    fileNames[[i]] <- read.table(myFiles[i], sep = " ")
}  

merged <- merge(fileNames[[1]], fileNames[[2]], by.x = "V1", by.y = "V1", all = F, sort = F)

present <- rowSums(merged[,2:3])

zeros <- rowSums(merged[,2:3] == 0)

output <- cbind(present, zeros)

rownames(output) <- merged[,1]

output
              present zeros
Ha8_00040788        2     0
Ha4_00024045        2     0
Ha4_00025366        2     0
Ha16_00022130       1     1
Ha16_00023451       1     1
Ha8_00040789        2     0
Ha4_00025367        2     0
Ha4_00024046        0     2
Ha16_00022131       2     0
Ha16_00023452       2     0

To add more files with merge, you could add something like this:
merged <- merge(fileNames[[1]], fileNames[[2]], by.x = "V1", by.y = "V1", all = F, sort = F)

for(i in 3:30){
  merged <- merge(merged, fileNames[[i]], by.x = "V1", by.y = "V1", all = F, sort = F)
}

You'll have to adjust all the numbers and coordinates to fit your data.


---

