0

所以我有两个文件，一个 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 列的所有行进行比较？

最好的，塞尔格

4

1 回答 1

2

我希望有所帮助（编辑我更改代码以获得更有效的算法）

gawk '
  #read input.genes and create list of limits (min, max)
  NR == FNR {
    #without header in input
    if(NR>1) {
      for(i=$2; i<=$3; i++){
        limits[i]=limits[i]","$2"-"$3;
      }
    };
    next
  }
  #read input.vcf, if column 3 is range of limits then print
  {
    if($3 in limits){
      print $0, "between("limits[$3]")"
    }
  }' input.genes input.vcf

你得到：

88  Chr1    25  C   -   3   2   1   1 between(,11-155)
88  Chr1    88  A   T   7   2   1   1 between(,11-155)
88  Chr1    92  A   C   16  4   1   1 between(,11-155)

python中的这个算法使用字典针对非常大的文件进行了优化

limits = [line.strip().split() for line in open("input.genes")]
limits.pop(0) #remove the header
limits = [map(int,v[1:]) for v in limits]

dict_limits = {}
for start, finish in limits:
  for i in xrange(start, finish+1):
    if i not in dict_limits:
      dict_limits[i] = []
    dict_limits[i].append((start,finish))

OUTPUT = open("my_output.txt", "w")
for reg in open("input.vcf"):
  v_reg = reg.strip().split()
  if int(v_reg[2]) in dict_limits:
    OUTPUT.write(reg.strip() + "\tbetween({})\n".format(str(dict_limits[int(v_reg[2])])))

OUTPUT.close()

你得到：

88 Chr1 25 C - 3 2 1 1 介于([(11, 155)])
88 Chr1 88 AT 7 2 1 1 介于([(11, 155)])
88 Chr1 92 AC 16 4 1 1 介于([(11, 155)])
于 2015-05-12T13:46:57.070 回答