以下の人口対立遺伝子数データがあります。
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行は各母集団の対立遺伝子「1」のコピー数を表し、1行は対立遺伝子「2」のコピー数を表します)。上の例では、1行目と2行目はSNP1の対立遺伝子1と2の数、3行目と4行目はSNP2の対立遺伝子の数であり、データセット全体でも同様です。すべての母集団について、各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ソリューション(クールなサードパーティパッケージなし)です。
#!/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]))
それはあなたに効果がありますか?フォーマットが少しずれても調整しやすいでしょう。