2

I have a population allele count data, 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. There are 2 rows per SNP (one row for the number of copies of allele "1" in each population and one row for the number of copies of allele "2"). In the example above, the first and second rows are number of alleles 1 and 2 for the SNP1, the third and fourth rows are alleles for SNP 2, and likewise for the entire dataset. I want to calculate population allele frequency per SNP for all populations: frequency of allele 1 at SNP1 in population 1 =Number of copies of allele 1 in population/Total number of 1+2 alleles copies in population and frequency of allele 2 at SNP1 in population 1 =Number of copies of allele a in population/Total number of 1+2 gene copies in population This means for each SNP number of copies of each allele should be divided to sum of the number of copies of allele "1" and "2" 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

I have R and bash solution, but is there any way to do this estimation in perl or python? Any help is appreciated.

I have the bash solution, here

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'

But do not know how to translate it into perl or python>

  • 3
    add the r/bash solution to question... and if you have those solution, you should be able to try something in perl/python.. they are not particularly difficult language to pick up – Sundeep Nov 7 '17 at 5:42
  • thanks @Sundeep. I included bash solution in my post. I am just learning perl and not so skilled yet to be able to translate from other languages. – Anna1364 Nov 7 '17 at 15:16
  • 1
    To translate awk and sed to perl, try 'a2p' and 's2p'. The result should – D McKeon Nov 8 '17 at 5:18
1

Here's a basic Python solution (no fancy 3rd party packages) that should be pretty close to what you're looking for:

#!/usr/bin/env python2

# snp.py

import sys

# Get the name of the data file from the command-line
data_file = sys.argv[1]

# Read and parse the data from the text file
data = []
with open(data_file, 'r') as file_handle:
    for line in file_handle:
        data.append([float(n) for n in line.split()])

# Get the number of rows and columns
rows = len(data)
cols = len(data[0])

# Iterate over adjacent pairs of rows
for r in range(rows//2):

    # Iterate over columns
    for c in range(cols):

        # Compute the sum of the two matching entries the pair of rows
        t = data[2*r][c] + data[2*r+1][c]
        if t:

            # Divide each entry by the sum of the pair
            data[2*r][c] /= t
            data[2*r+1][c] /= t

# Convert the data array back into formatted strings and print the results
for row in data:
    print(' '.join(['{0: <8}'.format(round(x,6)) for x in row]))

Does that work for you? It should be easy to tweak if it's a little off as far as formatting.

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.