【发布时间】:2015-07-23 08:11:04
【问题描述】:
所以我有两个文件,一个看起来像的 VCF
88 Chr1 25 C - 3 2 1 1
88 Chr1 88 A T 7 2 1 1
88 Chr1 92 A C 16 4 1 1
另一个基因看起来像
GENEID Start END
GENE_ID 11 155
GENE_ID 165 999
我想要一个脚本来查看第二个文件的第二个和第三个位置范围内是否存在基因位置(VCF 文件的第三列),然后将其打印出来。
到目前为止我所做的是加入文件并做
awk '{if (3>$12 && $3< $13) print }' > out
我所做的只是比较当前连接文件的行(它只在值在同一行时打印),我怎样才能让它比较第 3 列的所有行和第 12 列和第 13 列的所有行?
最好, 塞尔格
【问题讨论】:
-
当发生这种情况时,您想打印 what 吗?您想将 VCF 文件中的每一行与基因 id 文件中的每一行进行比较吗?没有必要为此加入文件(正如您所见,这并不能帮助您进行交叉匹配,并且实际上会使这变得更加困难)。
-
这些文件有多大?他们有多少行?十,一百万?
-
有一个基因文件和数百个VCF文件。基因一有大约 1000 行,而 VCF 有更多……范围从 5000-20000
-
您需要将基因坐标存储在内存中,然后循环通过 VCF 并根据每个基因的坐标检查每个变体的位置。这对于 awk 可能会很棘手,但对于 python 来说很简单。当然,这将是低效的,因为您必须对每条 VCF 行进行 2000 次比较。你不能使用像 BEDtools 这样的东西来做到这一点吗?
-
@heathobrien 我刚刚开始处理数据,我之前没有使用过 BEDtools。我现在正在查找它,看看它是否对我有用:)
标签: python awk sed bioinformatics vcf-variant-call-format