<p>I don't know how fast you need it to be, but this R solution runs on 20 million SNPs in under 2 min on my laptop:</p>

<pre><code>allele_counts &lt;-matrix(rep(c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0, 0, 0,
                     0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 0, 0,
                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                     2, 2, 0, 0, 2, 1, 0, 0, 0, 0, 2, 4, 0, 0, 0, 2, 0), 10000000), 
                   ncol=17, byrow=TRUE)
allele1&lt;-allele_counts[ c(TRUE,FALSE), ]
allele2&lt;-allele_counts[ c(FALSE,TRUE), ]
total_counts &lt;- allele1+allele2
allele_counts[c(TRUE,FALSE), ] &lt;-allele1/total_counts
allele_counts[c(FALSE,TRUE), ] &lt;-allele2/total_counts
</code></pre>

<p><em>Edit</em></p>

<p>I decided to a head-to-head test of this solution vs. <a href="https://bioinformatics.stackexchange.com/a/2795/581">bli's pandas solution</a> and your awk solution</p>

<p>I used readr for file IO:</p>

<pre><code>library(readr)
allele_counts&lt;-as.matrix(read_tsv("allele_counts.txt", col_names=F))
...
write_tsv(as.data.frame(allele_counts), "allele_freq_R.txt", col_names=F)
</code></pre>

<p>Here are the timings:</p>

<pre><code>cat Benchmarks/R.txt 
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
103.7556        0:01:43 2211.16 2518.59 2209.11 2209.15 0.00    27.58   73.62

cat Benchmarks/Python.txt 
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
132.6425        0:02:12 4872.45 6042.28 4870.89 4870.92 2.67    46.59   99.92
cat Benchmarks/awk.txt 
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
126.4455        0:02:06 4.22    413.55  0.93    0.97    0.00    130.68  91.28
</code></pre>

<p>The R solution is slightly faster than the pandas solution and uses half as much memory. Your awk one-liner is about the same speed as the pandas script, but requires a tiny fraction of the memory footprint.</p>

CC BY-SA 3.0
2017-11-07T16:05:56.720Z
Intelligent use of vector reusing! I would recommend scale  though
llrs
CC BY-SA 3.0
2017-11-07T16:08:35.760Z
Thanks. I nicked the idea from [here](https://stackoverflow.com/questions/24440258/selecting-multiple-odd-or-even-columns-rows-for-dataframe)
heathobrien
CC BY-SA 3.0
2017-11-09T09:42:59.883Z
Thanks for the benchmarks. Maybe you can answer this question of mine, about the meaning of the benchmarking values given by psutil (via snakemake): https://stackoverflow.com/q/46813371/1878788
bli
CC BY-SA 3.0
2017-11-09T11:59:33.597Z
I'm new new to benchmarking with Snakmake and also couldn't find any useful documentation for it. I did my best to answer with the info I've found, but I guess it's a case of the blind leading the blind
heathobrien
CC BY-SA 3.0
2017-11-07T16:01:29.857Z
2017-11-08T16:21:06.937Z
2017-11-08T16:21:06.937Z
heathobrien
<p>Here is an implementation in pure python. Probably not the fastest, but you should not encounter any memory issues with this.</p>

<p>Note: this is python3, use as follows: <code>cat snp_counts.tsv | python allele_frequencies.py</code></p>

<pre><code>#!/usr/bin/env python
"""
allele_frequencies.py
"""
import sys

is_odd = True
for line in sys.stdin:
    parts = line.strip().split()
    if not parts:
        continue
    allele_counts = [ int(p) for p in parts ]
    if is_odd:
        allele1 = allele_counts
        is_odd = False
    else:
        allele2 = allele_counts
        allele_sum = [ sum([a1,a2]) for a1,a2 in zip(allele1,allele2) ]
        allele1_freq = [ a / tot if tot else 0.0 for a,tot in zip(allele1, allele_sum) ]
        allele2_freq = [ a / tot if tot else 0.0 for a,tot in zip(allele2, allele_sum) ]
        print(*allele1_freq, sep = '\t')
        print(*allele2_freq, sep = '\t')
        is_odd = True
</code></pre>

CC BY-SA 3.0
2018-01-28T19:19:42.027Z
thanks for the code. I tried your code but I get an error message:  File "./allele_frequencies.py", line 21
    print(*allele1_freq, sep = '\t')
          ^
SyntaxError: invalid syntax.
Anna1364
CC BY-SA 3.0
2018-02-05T14:46:56.847Z
@Anna1364 Are you sure you are using python3?
holmrenser
CC BY-SA 3.0
2017-12-15T09:10:32.840Z
2017-12-15T09:10:32.840Z
holmrenser
<p>In python, if your files can fit in RAM, you might try using pandas and cytoolz.</p>

<pre class="lang-python prettyprint-override"><code>import pandas as pd
from cytoolz import interleave

def odd(num):
    return num % 2

def even(num):
    return 1 - (num % 2)

# Read in two separate dataframes
alleles_1 = pd.read_table("snps.txt", skiprows=odd, sep="\s+", header=None)
alleles_2 = pd.read_table("snps.txt", skiprows=even, sep="\s+", header=None)
sums = alleles_1 + alleles_2
# https://stackoverflow.com/a/45565745/1878788
pd.DataFrame(list(interleave([
    alleles_1.div(sums).fillna(0).values,
    alleles_2.div(sums).fillna(0).values]))).to_csv(
        "freqs.txt", sep="\t", index=False, header=False)
</code></pre>

<p>I don't know how efficient this is. On my workstation it runs in 1m27 on data obtained by taking 1 million times the example given in the question. See <a href="https://bioinformatics.stackexchange.com/a/2794/292">heathobrien's answer</a> for benchmarks on larger data.</p>

CC BY-SA 3.0
2017-11-08T16:00:21.867Z
Interesting, though I actually used the example times 10 million. I added some proper benchmarking to my anwser
heathobrien
CC BY-SA 3.0
2017-11-09T09:51:50.180Z
@heathobrien Indeed, I had misread your code. With 20 million snps, on my workstation where Firefox is eating one third of the 16GB RAM, my script fails on a MemoryError after 4 minutes.
bli
CC BY-SA 3.0
2017-11-07T18:36:03.367Z
2017-11-09T09:54:43.703Z
2017-11-09T09:54:43.703Z
bli
<p>Here's my quick Perl attempt after misreading the question, which divides each row elements by the row total:</p>

<pre><code>cat popCountData.txt | perl -lane '
  my $total=0;
  grep {$total += $_} @F; ## create per-row total
  foreach (@F) {printf("%0.2f ", ($total) ? $_ / $total : 0)} @F; ## divide elements by total
  print ""; ## print a line break (note: line will end with a space)
'
0.20 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.20 0.40 0.20 0.00 0.00 0.00 0.00 
0.00 0.25 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.50 0.00 0.25 0.00 0.00 0.00 
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
0.13 0.13 0.00 0.00 0.13 0.07 0.00 0.00 0.00 0.00 0.13 0.27 0.00 0.00 0.00 0.13 0.00
</code></pre>

<p>The challenge here is that the processing is in chunks of two lines, which means a couple of things:</p>

<ol>
<li>The previous line (at least for even-numbered lines) needs to be stored</li>
<li>The sum of counts (or whatever statistic is necessary) will have different processing on each line</li>
</ol>

<p>Unfortunately, I can't quite see how to get from your example input to your example output. For example, how does a count of 2 and 4 turn into a frequency of $1/3$ and $2/3$?</p>

CC BY-SA 3.0
2017-11-07T08:24:01.640Z
Because in that colum there is a 2+4 = 6 So 2/6 = 1/3 and 4/6 = 2/3, The scale is by each two rows and column wise.
llrs
CC BY-SA 3.0
2017-11-07T15:20:12.217Z
@gringer, thanks but it does not give the output that I want. I added bash solution.Yes, LIopis is right
Anna1364
CC BY-SA 3.0
2017-11-07T08:17:52.723Z
2017-11-07T08:17:52.723Z
gringer
<p>I already know the Bash AWK solution and also R code for the question that I am asking. But, my file is so big and with R it takes very long time and I am afraid of AWK's mistake for this computation, so I wonder if there is a Perl or Python solution for my problem.</p>

<p>I have a population allele count data, and it looks like this:</p>

<pre><code>1   0   0   0   0   0   0   0   0   0   1   2   1   0   0   0   0
0   2   0   0   0   0   0   0   0   0   0   4   0   2   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
2   2   0   0   2   1   0   0   0   0   2   4   0   0   0   2   0
</code></pre>

<p>Columns are populations and rows are SNPs. I have two rows per SNP (one row for the number of copies of allele "A" in each population and one row for the number of copies of allele "a").</p>

<p>In the example above, the first and second rows are number of alleles for the SNP1, the third and fourth rows are alleles for SNP 2, and like this until SNP20000000.</p>

<p>I want to calculate population allele frequency per SNP for all populations: frequency of allele A at SNP1 in population 1 = number of copies of allele A in population/total number of A/a gene copies in population and frequency of allele a at SNP1 in population 1 = number of copies of allele a in population/total number of A/a gene copies in the population. This means for each SNP number of copies of each allele should be divided to sum of the number of copies of allele "A" and "a" for each population. This is my desired output:</p>

<pre><code>1   0   0   0   0   0   0   0   0   0   1   0.333333    1   0   0   0   0
0   1   0   0   0   0   0   0   0   0   0   0.666667    0   1   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0           0   0   0   0   0
1   1   0   0   1   1   0   0   0   0   1   1           0   0   0   1   0
</code></pre>

<p>As I said, I have an R and Bash solution, but is there a way to do this estimation in Perl or Python?</p>

<p>Here is the Bash solution, but I do not know how to translate it into Perl or Python? I am learning Perl...</p>

<pre><code>awk '{for(i = 1; i &lt;= NF; i++)
    tally[i] += $i};
    (NR%2) == 1
    {for(i=1; i&lt;=NF; i++)
        allele1[i] = $i};
        (NR%2) == 0
        {for(i = 1; i &lt;= NF; i++)
        allele2[i] = $i;
            for(i = 1; i &lt;= NF; i++)
                if(tally[i] == 0)
                    tally[i] = 1;
                for(i = 1; i &lt;= NF; i++)
                    printf allele1[i]/tally[i]"\t";
                    printf "\n";
                for(i = 1; i &lt;= NF; i++)
                    printf allele2[i]/tally[i]"\t";
                    printf "\n";
                for(i = 1; i &lt;= NF; i++)
                    tally[i]=0}' MyData |
sed 's/\t$//g'
</code></pre>

CC BY-SA 3.0
2017-11-07T07:54:59.500Z
It is strange that with R it takes so much time. Usually with matrices it is quite fast. Maybe you can post the code in R to translate it to python or perl, or improve it. perl is usually better for text editing, as with numbers it is a bit more difficult but it is quite fast usually. What have you tried in those languages? Did you translate your R code? Did you start from zero?
llrs
CC BY-SA 3.0
2017-11-07T15:17:17.693Z
File is 8GB. I included bash solution but do not know exactly how to translate it into perl or python!
Anna1364
CC BY-SA 3.0
2017-11-07T15:27:32.547Z
If you have a solution that works in on language then I dont' think this is the right forum.  It is no longer a bioinformatics question it is a programming question that happens to use biological data.
Bioathlete
CC BY-SA 3.0
2017-11-07T15:28:57.930Z
It seems like you do 7 loops in your *awk* code. I am not an expert of awk, but this calls for an optimization, I think it could be reduced to 4 or so. How is your code in R? (I am more confident of my R skills, and perhaps I could help you) And as Bioathlete says this might suit better to code Reviw. Where they could help you speed and improve your code
llrs
CC BY-SA 3.0
2017-11-07T15:30:32.647Z
Could you also share the R solution for for those of us that find awk one liners to be completely impenetrable?
heathobrien
CC BY-SA 3.0
2017-11-07T15:48:27.587Z
@ LIopis, OK. what is code review? is that what you meant? is that a forum?
Anna1364
CC BY-SA 3.0
2017-11-07T16:05:28.180Z
@Anna, code review is when someone looks your code to improve it. See [this site](https://codereview.stackexchange.com) in the stackexchange network.
llrs
CC BY-SA 3.0
2017-11-07T17:24:22.013Z
Your "bash" code is an `awk` solution
gringer
<python><perl>
Convert population allele count into population allele frequency by using Perl or Python (translate available Bash into Perl or Python)
# Convert population allele count into population allele frequency by using Perl or Python (translate available Bash into Perl or Python)

(Asked by: username_0 on 2017-11-07)

I already know the Bash AWK solution and also R code for the question that I am asking. But, my file is so big and with R it takes very long time and I am afraid of AWK's mistake for this computation, so I wonder if there is a Perl or Python solution for my problem.
I have a population allele count data, and it looks like this:
1   0   0   0   0   0   0   0   0   0   1   2   1   0   0   0   0
0   2   0   0   0   0   0   0   0   0   0   4   0   2   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
2   2   0   0   2   1   0   0   0   0   2   4   0   0   0   2   0

Columns are populations and rows are SNPs. I have two rows per SNP (one row for the number of copies of allele "A" in each population and one row for the number of copies of allele "a").
In the example above, the first and second rows are number of alleles for the SNP1, the third and fourth rows are alleles for SNP 2, and like this until SNP20000000.
I want to calculate population allele frequency per SNP for all populations: frequency of allele A at SNP1 in population 1 = number of copies of allele A in population/total number of A/a gene copies in population and frequency of allele a at SNP1 in population 1 = number of copies of allele a in population/total number of A/a gene copies in the population. This means for each SNP number of copies of each allele should be divided to sum of the number of copies of allele "A" and "a" for each population. This is my desired output:
1   0   0   0   0   0   0   0   0   0   1   0.333333    1   0   0   0   0
0   1   0   0   0   0   0   0   0   0   0   0.666667    0   1   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0           0   0   0   0   0
1   1   0   0   1   1   0   0   0   0   1   1           0   0   0   1   0

As I said, I have an R and Bash solution, but is there a way to do this estimation in Perl or Python?
Here is the Bash solution, but I do not know how to translate it into Perl or Python? I am learning Perl...
awk '{for(i = 1; i <= NF; i++)
    tally[i] += $i};
    (NR%2) == 1
    {for(i=1; i<=NF; i++)
        allele1[i] = $i};
        (NR%2) == 0
        {for(i = 1; i <= NF; i++)
        allele2[i] = $i;
            for(i = 1; i <= NF; i++)
                if(tally[i] == 0)
                    tally[i] = 1;
                for(i = 1; i <= NF; i++)
                    printf allele1[i]/tally[i]"\t";
                    printf "\n";
                for(i = 1; i <= NF; i++)
                    printf allele2[i]/tally[i]"\t";
                    printf "\n";
                for(i = 1; i <= NF; i++)
                    tally[i]=0}' MyData |
sed 's/\t$//g'

**Tags:** python, perl
> **Comments:**
> **username_1**: It is strange that with R it takes so much time. Usually with matrices it is quite fast. Maybe you can post the code in R to translate it to python or perl, or improve it. perl is usually better for text editing, as with numbers it is a bit more difficult but it is quite fast usually. What have you tried in those languages? Did you translate your R code? Did you start from zero?  
> **username_0**: File is 8GB. I included bash solution but do not know exactly how to translate it into perl or python!  
> **username_2**: If you have a solution that works in on language then I dont' think this is the right forum.  It is no longer a bioinformatics question it is a programming question that happens to use biological data.  
> **username_1**: It seems like you do 7 loops in your *awk* code. I am not an expert of awk, but this calls for an optimization, I think it could be reduced to 4 or so. How is your code in R? (I am more confident of my R skills, and perhaps I could help you) And as Bioathlete says this might suit better to code Reviw. Where they could help you speed and improve your code  
> **username_3**: Could you also share the R solution for for those of us that find awk one liners to be completely impenetrable?  
> **username_0**: @ LIopis, OK. what is code review? is that what you meant? is that a forum?  
> **username_1**: @Anna, code review is when someone looks your code to improve it. See [this site](https://codereview.stackexchange.com) in the stackexchange network.  
> **username_4**: Your "bash" code is an `awk` solution  

---

## 4 answers

### Answer by username_3 on 2017-11-07. Score: 4

I don't know how fast you need it to be, but this R solution runs on 20 million SNPs in under 2 min on my laptop:
allele_counts <-matrix(rep(c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0, 0, 0,
                     0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 0, 0,
                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                     2, 2, 0, 0, 2, 1, 0, 0, 0, 0, 2, 4, 0, 0, 0, 2, 0), 10000000), 
                   ncol=17, byrow=TRUE)
allele1<-allele_counts[ c(TRUE,FALSE), ]
allele2<-allele_counts[ c(FALSE,TRUE), ]
total_counts <- allele1+allele2
allele_counts[c(TRUE,FALSE), ] <-allele1/total_counts
allele_counts[c(FALSE,TRUE), ] <-allele2/total_counts

Edit
I decided to a head-to-head test of this solution vs. bli's pandas solution and your awk solution
I used readr for file IO:
library(readr)
allele_counts<-as.matrix(read_tsv("allele_counts.txt", col_names=F))
...
write_tsv(as.data.frame(allele_counts), "allele_freq_R.txt", col_names=F)

Here are the timings:
cat Benchmarks/R.txt 
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
103.7556        0:01:43 2211.16 2518.59 2209.11 2209.15 0.00    27.58   73.62

cat Benchmarks/Python.txt 
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
132.6425        0:02:12 4872.45 6042.28 4870.89 4870.92 2.67    46.59   99.92
cat Benchmarks/awk.txt 
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
126.4455        0:02:06 4.22    413.55  0.93    0.97    0.00    130.68  91.28

The R solution is slightly faster than the pandas solution and uses half as much memory. Your awk one-liner is about the same speed as the pandas script, but requires a tiny fraction of the memory footprint.

> **Comments:**
> **username_1**: Intelligent use of vector reusing! I would recommend scale  though  
> **username_3**: Thanks. I nicked the idea from [here](https://stackoverflow.com/questions/24440258/selecting-multiple-odd-or-even-columns-rows-for-dataframe)  
> **username_5**: Thanks for the benchmarks. Maybe you can answer this question of mine, about the meaning of the benchmarking values given by psutil (via snakemake): https://stackoverflow.com/q/46813371/1878788  
> **username_3**: I'm new new to benchmarking with Snakmake and also couldn't find any useful documentation for it. I did my best to answer with the info I've found, but I guess it's a case of the blind leading the blind  

---

### Answer by username_6 on 2017-12-15. Score: 0

Here is an implementation in pure python. Probably not the fastest, but you should not encounter any memory issues with this.
Note: this is python3, use as follows: cat snp_counts.tsv | python allele_frequencies.py
#!/usr/bin/env python
"""
allele_frequencies.py
"""
import sys

is_odd = True
for line in sys.stdin:
    parts = line.strip().split()
    if not parts:
        continue
    allele_counts = [ int(p) for p in parts ]
    if is_odd:
        allele1 = allele_counts
        is_odd = False
    else:
        allele2 = allele_counts
        allele_sum = [ sum([a1,a2]) for a1,a2 in zip(allele1,allele2) ]
        allele1_freq = [ a / tot if tot else 0.0 for a,tot in zip(allele1, allele_sum) ]
        allele2_freq = [ a / tot if tot else 0.0 for a,tot in zip(allele2, allele_sum) ]
        print(*allele1_freq, sep = '\t')
        print(*allele2_freq, sep = '\t')
        is_odd = True

> **Comments:**
> **username_0**: thanks for the code. I tried your code but I get an error message:  File "./allele_frequencies.py", line 21     print(*allele1_freq, sep = '\t')           ^ SyntaxError: invalid syntax.  
> **username_6**: @Anna1364 Are you sure you are using python3?  

---

### Answer by username_5 on 2017-11-07. Score: 0

In python, if your files can fit in RAM, you might try using pandas and cytoolz.
import pandas as pd
from cytoolz import interleave

def odd(num):
    return num % 2

def even(num):
    return 1 - (num % 2)

# Read in two separate dataframes
alleles_1 = pd.read_table("snps.txt", skiprows=odd, sep="\s+", header=None)
alleles_2 = pd.read_table("snps.txt", skiprows=even, sep="\s+", header=None)
sums = alleles_1 + alleles_2
# https://stackoverflow.com/a/45565745/1878788
pd.DataFrame(list(interleave([
    alleles_1.div(sums).fillna(0).values,
    alleles_2.div(sums).fillna(0).values]))).to_csv(
        "freqs.txt", sep="\t", index=False, header=False)

I don't know how efficient this is. On my workstation it runs in 1m27 on data obtained by taking 1 million times the example given in the question. See heathobrien's answer for benchmarks on larger data.

> **Comments:**
> **username_3**: Interesting, though I actually used the example times 10 million. I added some proper benchmarking to my anwser  
> **username_5**: @heathobrien Indeed, I had misread your code. With 20 million snps, on my workstation where Firefox is eating one third of the 16GB RAM, my script fails on a MemoryError after 4 minutes.  

---

### Answer by username_4 on 2017-11-07. Score: -1

Here's my quick Perl attempt after misreading the question, which divides each row elements by the row total:
cat popCountData.txt | perl -lane '
  my $total=0;
  grep {$total += $_} @F; ## create per-row total
  foreach (@F) {printf("%0.2f ", ($total) ? $_ / $total : 0)} @F; ## divide elements by total
  print ""; ## print a line break (note: line will end with a space)
'
0.20 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.20 0.40 0.20 0.00 0.00 0.00 0.00 
0.00 0.25 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.50 0.00 0.25 0.00 0.00 0.00 
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
0.13 0.13 0.00 0.00 0.13 0.07 0.00 0.00 0.00 0.00 0.13 0.27 0.00 0.00 0.00 0.13 0.00

The challenge here is that the processing is in chunks of two lines, which means a couple of things:

The previous line (at least for even-numbered lines) needs to be stored
The sum of counts (or whatever statistic is necessary) will have different processing on each line

Unfortunately, I can't quite see how to get from your example input to your example output. For example, how does a count of 2 and 4 turn into a frequency of $1/3$ and $2/3$?

> **Comments:**
> **username_1**: Because in that colum there is a 2+4 = 6 So 2/6 = 1/3 and 4/6 = 2/3, The scale is by each two rows and column wise.  
> **username_0**: @gringer, thanks but it does not give the output that I want. I added bash solution.Yes, LIopis is right  

---

