☰

如何将每行中的值除以两行的总和（将 bash 翻译为 perl 或 python）

tag-icon • tag-icon python perl
如何将每行中的值除以两行的总和（将 bash 翻译为 perl 或 python）

我有一个人口等位基因计数数据，看起来像这样

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

列是群体，行是 SNP。每个 SNP 有 2 行（一行表示每个群体中等位基因“1”的拷贝数，一行表示等位基因“2”的拷贝数）。在上面的示例中，第一行和第二行是 SNP1 的等位基因 1 和 2 的数量，第三行和第四行是 SNP 2 的等位基因数量，整个数据集也是如此。我想计算所有群体的每个 SNP 的群体等位基因频率： 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这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 意味着对于每个 SNP，每个等位基因的拷贝数应除以每个群体的等位基因“1”和“2”的拷贝数之和。这是我想要的输出：

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

我有 R 和 bash 解决方案，但是有什么方法可以在 perl 或 python 中进行此估计吗？任何帮助表示赞赏。

我这里有 bash 解决方案

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'

但不知道如何翻译成perl或python>

答案1

这是一个基本的 Python 解决方案（没有花哨的第 3 方包），应该非常接近您正在寻找的内容：

#!/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]))

那对你有用吗？如果格式有一点偏差，应该很容易调整。

相关内容