1
$\begingroup$

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'
$\endgroup$
  • 1
    $\begingroup$ 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? $\endgroup$ – llrs Nov 7 '17 at 7:54
  • $\begingroup$ File is 8GB. I included bash solution but do not know exactly how to translate it into perl or python! $\endgroup$ – Anna1364 Nov 7 '17 at 15:17
  • 1
    $\begingroup$ 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. $\endgroup$ – Bioathlete Nov 7 '17 at 15:27
  • 1
    $\begingroup$ 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 $\endgroup$ – llrs Nov 7 '17 at 15:28
  • 1
    $\begingroup$ Your "bash" code is an awk solution $\endgroup$ – gringer Nov 7 '17 at 17:24
4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Intelligent use of vector reusing! I would recommend scale though $\endgroup$ – llrs Nov 7 '17 at 16:05
  • $\begingroup$ Thanks. I nicked the idea from here $\endgroup$ – heathobrien Nov 7 '17 at 16:08
  • 1
    $\begingroup$ Thanks for the benchmarks. Maybe you can answer this question of mine, about the meaning of the benchmarking values given by psutil (via snakemake): stackoverflow.com/q/46813371/1878788 $\endgroup$ – bli Nov 9 '17 at 9:42
  • $\begingroup$ 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 $\endgroup$ – heathobrien Nov 9 '17 at 11:59
0
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ Interesting, though I actually used the example times 10 million. I added some proper benchmarking to my anwser $\endgroup$ – heathobrien Nov 8 '17 at 16:00
  • 1
    $\begingroup$ @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. $\endgroup$ – bli Nov 9 '17 at 9:51
0
$\begingroup$

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
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Anna1364 Jan 28 '18 at 19:19
  • $\begingroup$ @Anna1364 Are you sure you are using python3? $\endgroup$ – holmrenser Feb 5 '18 at 14:46
-1
$\begingroup$

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:

  1. The previous line (at least for even-numbered lines) needs to be stored
  2. 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$?

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – llrs Nov 7 '17 at 8:24
  • $\begingroup$ @gringer, thanks but it does not give the output that I want. I added bash solution.Yes, LIopis is right $\endgroup$ – Anna1364 Nov 7 '17 at 15:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.