将VCF测序数据中的等位基因频率结合起来

Question

我有一个测序数据文件，其中包含来自基因组的碱基对位置，如下所示：

chr1 814 G A 0.5
chr1 815 T A 0.3
chr1 816 C G 0.2
chr2 315 A T 0.3
chr2 319 T C 0.8
chr2 340 G C 0.3
chr4 514 A G 0.5

我想比较一下第2栏中bp的位置所定义的某些组，然后我想要匹配区域第5栏中数字的平均值。

因此，使用上面的示例，假设我正在寻找跨越chr1 810-820和chr2 310-330的所有样本的第5列的平均值。前五行应该被识别，它们的第五列编号应该是平均的，等于0.42。

我尝试创建一个范围数组，然后使用awk调用这些位置，但是没有成功。提前谢谢。

Answer 1

import pandas as pd
from StringIO import StringIO

s = """chr1 814 G A 0.5
chr1 815 T A 0.3
chr1 816 C G 0.2
chr2 315 A T 0.3
chr2 319 T C 0.8
chr2 340 G C 0.3
chr4 514 A G 0.5"""

sio = StringIO(s)
df = pd.read_table(sio, sep=" ", header=None)
df.columns=["a", "b", "c", "d", "e"]

# The query expression is intuitive 
r = df.query("(a=='chr1' & 810<b<820) | (a=='chr2' & 310<b<330)")
print r["e"].mean()

熊猫可能更适合这样的表格数据处理，而且它是蟒蛇。

Answer 2

这里有一些python代码，可以满足您的要求。它假设数据驻留在一个名为“data.txt”的文本文件中。

#!/usr/bin/env python

data = open('data.txt').readlines()
def avg(keys):
    key_sum = 0
    key_count = 0
    for item in data:
        fields = item.split()
        krange = keys.get(fields[0], None)
        if krange:
            r = int(fields[1])
            if krange[0] <= r and r <= krange[1]:
                key_sum += float(fields[-1])
                key_count += 1
    print key_sum/key_count

keys = {} # Create dict to store keys and ranges of interest
keys['chr1'] = (810, 820)
keys['chr2'] = (310, 330)

avg(keys)

样本输出：

0.42

Answer 3

这是一个awk脚本的答案。为了输入，我创建了第二个文件，我称之为ranges

chr1 810 820
chr2 310 330

脚本本身看起来像：

#!/usr/bin/awk -f

FNR==NR { low_r[$1] = $2; high_r[$1] = $3; next }

{ l = low_r[ $1 ]; h = high_r[$1]; if( l=="" ) next }

$2 >= l && $2 <= h { total+=$5; cnt++ }

END {
        if( cnt > 0 ) print (total/cnt)
        else print "no matched data"
}

细分如下：

• FNR==NR -吸收ranges文件，使low_r和high_r数组按键离开该文件的第一列。
• 然后，对于数据中的每一行，low_r和high_r数组中的查找都匹配。如果没有匹配，则跳过任何其他处理
• 检查基于low和high测试的包含范围，增加匹配范围的total和cnt。
• 在END中，当有匹配时，打印简单的平均值

当脚本(称为script.awk)变成可执行文件时，它可以运行如下：

$ ./script.awk ranges data
0.42

在这里，我调用了数据文件data。