将来自 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 列的平均值。应确定前五行，并将其第 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()

pandas可能更适合这样的表格数据处理，python.

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[] = ; high_r[] = ; next }

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

 >= l &&  <= h { total+=; 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。

将来自 VCF 测序数据的等位基因频率合并在一起

Binning Together Allele Frequencies From VCF Sequencing Data

bash

awk

vcf-vcard

sequencing

fastq

示例输出：