counting individuals with the present gene and absent gene for each population

Question

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

Answer 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.

Answer 2

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