【问题标题】:Comparing one column value to all columns in linux enviroment在linux环境中将一列值与所有列进行比较
【发布时间】: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


【解决方案1】:

希望对您有所帮助(编辑我更改代码以获得更高效的算法)

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 A T 7 2 1 1 介于([(11, 155)]) 88 Chr1 92 A C 16 4 1 1 介于([(11, 155)])

【讨论】:

  • 我试过了,根本没有输出:/尝试重定向到一个文件,结果是空的。你能解释一下为什么你在第五行有 [NR -1,1] = $2 吗? -1,1 在做什么?
  • @srx 因为你的input.genes 文件有标题,我必须删除它
  • 处理 for(i=1; i
  • @srx 是的for (i = 1; i in a ;i++ ) 不起作用,.....因为i in a 返回false
  • 如果我将其保留在原始命令中,它会显示第 11 行:对数组 a 的非法引用。似乎我在这台机器上缺少 GNU awk (GAWK) 它可以工作!非常感谢您的帮助!
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 2019-09-11
  • 1970-01-01
  • 2022-08-05
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多