0

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

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 调用这些位置，但没有成功。提前致谢。

4

3 回答 3

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 might be better for such tabular data processing, and it's python.

于 2015-04-06T17:09:09.157 回答
1

Here's some python code to do what you are asking for. It assumes that your data lives in a text file called '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)

Sample Output:

0.42
于 2015-04-06T17:09:22.177 回答
1

这是一个 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_randhigh_r数组中查找匹配项。如果没有匹配，则跳过任何其他处理
  • 根据low和high测试、递增total和cnt匹配范围检查包含范围。
  • 在 处END，打印匹配时的简单平均值

当脚本（称为script.awk）可执行时，它可以像这样运行：

$ ./script.awk ranges data
0.42

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

于 2015-04-06T18:58:30.157 回答