【问题标题】:awk code to filter lines in one file according to matching conditions in another fileawk 代码根据另一个文件中的匹配条件过滤一个文件中的行
【发布时间】:2020-05-14 03:59:29
【问题描述】:

我有两个如下所示的文件:

file1 有 511 行:

chr start end
1 1227897 2779043
1 6644723 8832944
1 11067792 11372913
1 17287414 17342924
1 23254576 23590651

file2 有 12594903 行:

CHR POS REF A1 OBS_CT BETA SE P ALT_FREQS RSID_UKB
1 58396 T C 382851 0.0882097 0.0677502 0.192923 0.000249012 rs570371753
1 91588 G A 382852 0.265908 0.0879796 0.00250811 0.000148375 rs554639997
1 713979 C G 382837 0.00630607 0.0925289 0.945664 0.000138059 rs117217250
1 715265 C T 377557 0.00260829 0.00617561 0.672768 0.0331599 rs12184267
1 715367 A G 377954 0.00212642 0.00615857 0.729886 0.0333038 rs12184277
1 717485 C A 377980 0.00449142 0.00615965 0.465899 0.0332908 rs12184279
1 6702159 G T 378749 0.00305772 0.00604916 0.613223 0.0345562 rs116801199
1 9902231 G C 378573 0.00216983 0.00607117 0.720793 0.0342995 rs12565286
1 23364524 C G 377155 0.00505093 0.00588132 0.390447 0.0368034 rs2977670

我正在尝试删除 file2 中的所有行,其中 file2 中的 POS(即位置)位于 file1 的 startend 之间,并且 file2 中的 CHR 等于 file1 中的 chr .我可以用 R 和下面的代码做到这一点,但它需要大约。以这种方式运行它需要 1 小时:

for (row in 1:nrow(file1)) {
  file2 <- file2[!(file2$CHR == file1$chr[row] &
                      file2$POS >= file1$start[row] &
                      file2$POS <= file1$end[row])]
}

输出应如下所示,删除了满足条件的 2 行(第 7 行和第 9 行):

CHR POS REF A1 OBS_CT BETA SE P ALT_FREQS RSID_UKB
1 58396 T C 382851 0.0882097 0.0677502 0.192923 0.000249012 rs570371753
1 91588 G A 382852 0.265908 0.0879796 0.00250811 0.000148375 rs554639997
1 713979 C G 382837 0.00630607 0.0925289 0.945664 0.000138059 rs117217250
1 715265 C T 377557 0.00260829 0.00617561 0.672768 0.0331599 rs12184267
1 715367 A G 377954 0.00212642 0.00615857 0.729886 0.0333038 rs12184277
1 717485 C A 377980 0.00449142 0.00615965 0.465899 0.0332908 rs12184279
1 9902231 G C 378573 0.00216983 0.00607117 0.720793 0.0342995 rs12565286

必须有一个简洁快速的 awk 单行代码来执行此操作。不过,我不知道怎么做。任何帮助将不胜感激!

【问题讨论】:

  • 给定样本数据的预期输出是什么?
  • 所以输出应该和file2一模一样,但是那些符合条件的行应该被删除。
  • 好的,请将预期的输出添加到问题中。我想要一些东西来比较我的结果。谢谢!
  • 啊,当然。已添加预期输出!

标签: r awk


【解决方案1】:

这是一种非常快速而丑陋的方法。 At 只是创建一个由所有可能性索引的列表:

awk '(NR==FNR) {for(i=$2;i<=$3;++i) a[$1,i]; next}(FNR==1);!(($1,$2) in a)' file1 file2

问题是,虽然写得很快,但速度很慢;-)

一种更快的方法可能是:

awk '(NR==FNR) { a[$1]=a[$1] FS $2 FS $3; next}
     { t=$0;c=$1;p=$2; $0=a[c] }
     { for(i=1;i<=NF;i+=2) if( $i<=p && p<=$(i+1) ) next; print t }' file1 file2

在后者中,我们创建如下列表:

a[1] = "1227897 2779043 6644723 8832944 11067792 11372913"

每次读取file2的新记录时,我们将记录存储在t下并跟踪CHR值(c=$1)。通过将a[c] 分配给$0,我们重新计算了所有字段,所以$1=1227897$2=2779043 等。所以我们要做的就是在所有字段中以两个为一组跳转并检查当前记录是否在各自的范围与否。如果它在范围内,我们跳到下一条记录。

正如comment of SiegeX 中提到的,如果file2 已排序,则可以加快速度:

 awk '(NR==FNR) { a[$1]=a[$1] FS $2 FS $3 FS a[$1]; next}
     { t=$0;c=$1;p=$2; $0=a[c] }
     { for(i=NF-1;i>0;i+=2) {
          if ( $(i+1) < p ) { $i=$(i+1)=""; $0=$0; $1=$1; a[c]=$0 }
          else if ( $i<=p && p<=$(i+1) ) { next }
          print t
     }' file1 file2

【讨论】:

  • 您的输出与他正在寻找的相反。您正在打印符合条件的行应该被删除
  • 快得多,非常感谢!现在我只需要启动我的 awk 游戏,直到我明白这里发生了什么。
  • @kvantour 巧妙地制作了以CHR 变量为关键字的记录。获得对数加速的一种方法是,只要其中一条记录的POS 超出一对范围,就剔除范围的哈希映射。例如,一旦到达POS 等于6702159 的记录,您就可以从a[1] 中删除范围对1227897 2779043,这样您就不再需要与它比较进一步的记录。这仅在 file2 按POS 排序时才有效,但从我们给出的示例文件来看,它似乎是这种方式。
  • @SiegeX 谢谢。我已经添加了一个更新。由于字段如何重新定义 $0,我不得不做一些恶作剧来获得所需的结果,但我认为应该这样做。
  • @kvantour @SiegeX 是的,file2 确实按CHRPOS 排序
【解决方案2】:
awk '
  NR==FNR{
    a[$2","$3]=$1
    next
  } 

  {
    for(i in a){
      split(i,b,",")
      if($1 == a[i] && $2 >= b[1] && $2 <= b[2]){next}
    }
  }1
' file1 file2

概念证明

注意:我在 file2 中添加了一个额外的输入行,其中 CHR1 不同,以测试该条件

$ cat file2
CHR POS REF A1 OBS_CT BETA SE P ALT_FREQS RSID_UKB
1 58396 T C 382851 0.0882097 0.0677502 0.192923 0.000249012 rs570371753
1 91588 G A 382852 0.265908 0.0879796 0.00250811 0.000148375 rs554639997
1 713979 C G 382837 0.00630607 0.0925289 0.945664 0.000138059 rs117217250
1 715265 C T 377557 0.00260829 0.00617561 0.672768 0.0331599 rs12184267
1 715367 A G 377954 0.00212642 0.00615857 0.729886 0.0333038 rs12184277
1 717485 C A 377980 0.00449142 0.00615965 0.465899 0.0332908 rs12184279
1 6702159 G T 378749 0.00305772 0.00604916 0.613223 0.0345562 rs116801199
1 9902231 G C 378573 0.00216983 0.00607117 0.720793 0.0342995 rs12565286
1 23364524 C G 377155 0.00505093 0.00588132 0.390447 0.0368034 rs2977670
2 23364524 C G 377155 0.00505093 0.00588132 0.390447 0.0368034 rs2977670

$ awk 'NR==FNR{a[$2","$3]=$1;next}{for(i in a){split(i,b,",");\
         if($1 == a[i] && $2 >= b[1] && $2 <= b[2]){next}}}1' file1 file2
CHR POS REF A1 OBS_CT BETA SE P ALT_FREQS RSID_UKB
1 58396 T C 382851 0.0882097 0.0677502 0.192923 0.000249012 rs570371753
1 91588 G A 382852 0.265908 0.0879796 0.00250811 0.000148375 rs554639997
1 713979 C G 382837 0.00630607 0.0925289 0.945664 0.000138059 rs117217250
1 715265 C T 377557 0.00260829 0.00617561 0.672768 0.0331599 rs12184267
1 715367 A G 377954 0.00212642 0.00615857 0.729886 0.0333038 rs12184277
1 717485 C A 377980 0.00449142 0.00615965 0.465899 0.0332908 rs12184279
1 9902231 G C 378573 0.00216983 0.00607117 0.720793 0.0342995 rs12565286
2 23364524 C G 377155 0.00505093 0.00588132 0.390447 0.0368034 rs2977670

【讨论】:

  • 谢谢,似乎工作正常,但并不比我的 R 代码快。
猜你喜欢
  • 2023-01-20
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2020-05-28
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多