【问题标题】:Filter rows of files conditional on multiple arrays values根据多个数组值过滤文件行
【发布时间】:2018-11-29 14:47:24
【问题描述】:

我有许多带有 qtl 摘要数据的文件 (N>1000),例如让我们假设第一个文件由六行组成(实际上它们都是 GWA/具有 >10M SNP 的估算文件)

cat QTL.1.txt 
Chr Rs  BP  beta    se  pvalue
11  rs11224233  134945522   0.150216    0.736939    0.962375  
11  rs4616056   134945709   0.129518    0.371824    0.910326   
11  rs11823417  134945710   0.103462    0.41737 0.845826  
11  rs80294507  134945765   0.150336    0.735363    0.961403  
11  rs61907173  134946034   0.104531    0.158224    0.884548  
11  rs147621717 134946277   0.105365    0.196168    0.86476

我想根据染色体和基因列表的位置过滤每个数据集(我的列表有 100 个基因,但现在不要假设它有 2 个);因此创建 N_QTL*N_Genes 文件。我想检查每个 QTL 的每个基因/位置。基因的染色体、位置和名称存储在四个数组中,我想迭代地读取这些数组并保存每个基因的每个 qtl 文件的输出。

到目前为止我所做的没有用,我知道 awk 不是最好的方法:

    declare -a array1    
    declare -a array2  
    declare -a array3  
    declare -a array4 

    array1=(11 11) #chromosome   
    array2=(134945709 134945765) #start gene position  
    array3=(134946034 134946277) #end gene position  
    array4=(A B) # gene name  


    for qtl in 1; do # in reality it would be for qtl in 1 1000 

         for ((i=0; i<${#array1[@]}; i++)); do  

        cat QTL.$qtl.txt | awk '$1=='array1[$i]' && $3>='array2[$i]' && 
        $3<='array3[$i]' {print$0}' >  Gene.${array4[$i]}_QTL.$qtl.txt; 
        done; 
    done

在 awk 中,$1 是染色体,$3 是位置,因此基于这些进行过滤。

所以我对基因 A 的 QTL.1.txt 的预期输出是

cat Gene.A_QTL.1.txt
Chr Rs  BP  beta    se  pvalue  
11  rs4616056   134945709   0.129518    0.371824    0.910326   
11  rs11823417  134945710   0.103462    0.41737 0.845826  
11  rs80294507  134945765   0.150336    0.735363    0.961403  
11  rs61907173  134946034   0.104531    0.158224    0.884548  

对于基因 B 的 QTL.1.txt 将是

cat Gene.B_QTL.1.txt
Chr Rs  BP  beta    se  pvalue
11  rs80294507  134945765   0.150336    0.735363    0.961403  
11  rs61907173  134946034   0.104531    0.158224    0.884548  
11  rs147621717 134946277   0.105365    0.196168    0.86476

我最终得到的是空文件,因为我要求根据数组的值过滤这些列的方式不起作用。

非常感谢任何帮助! 提前谢谢你

【问题讨论】:

  • 欢迎来到 SO,您能否将您的输入和输出发布在 CODE TAGS 中,以便更好地了解我们,以便为您提供指导。
  • 您好,谢谢您的建议。这样更好吗?
  • 请在您的帖子中也添加预期的输出。
  • 应该 array2[$] 是 array2[$i]
  • 谢谢,我编辑了array2并添加了预期的输出。

标签: arrays loops awk conditional


【解决方案1】:

混合使用 bash 和 awk 来解析文件并不总是最好的方法。

这里只有 awk 的解决方案。

假设您已将信息分配给文件中的 bash 数组:

$ cat info
11 134945765 154945765 Gene1
12 134945522 174945522 Gene2

您可以使用以下 awk 脚本对数据文件执行查找:

awk 'NR==FNR{
       for(i=2;i<=NF;i++) 
         a[$1,i]=$i
         next
     } 
     a[$1,2]<=$3 && a[$1,3]>=$3{
        print $0 > a[$1,4]"_QTL"
     }' info QTL.1.txt

这将创建一个包含以下内容的文件:

$ cat Gene1_QTL
11 rs80294507 134945765 0.150336 0.735363 0.961403
11 rs61907173 134946034 0.104531 0.158224 0.884548
11 rs147621717 134946277 0.105365 0.196168 0.86476

也许不完全是您正在查看的内容,但我希望这对您有所帮助...

【讨论】:

  • 非常感谢您!我对 a[$1,4]"_QTL" 有一个小问题——它不喜欢 "_QTL" 但没有它也可以工作。但是,我需要通过基因列表和每个 QTL 运行它。是否可以只使用 awk?
  • 用文件名改变“_QTL”。
  • 嗨 karafka,非常感谢您的建议,但我仍然收到一个错误:awk: syntax error at source line 7 context is print $0 > >>> array[$1,4]FILENAME
【解决方案2】:

如果多个基因位于同一染色体中,您可能需要执行以下操作(使用基因名称而不是 chr 作为 Key):

awk 'NR==FNR{
        chr[$4]=$1;
        start[$4]=$2;
        end[$4]=$3;
     } 
     NR!=FNR{
        for (var in chr){
            name=var"_"FILENAME;
            if(chr[var]==$1 && start[var] <=$3 && end[var]>=$3){
                print $0 > name;
            }
        }
    }'  info QTL

【讨论】:

  • 非常感谢!
猜你喜欢
  • 2017-03-13
  • 1970-01-01
  • 2017-08-03
  • 1970-01-01
  • 2019-02-17
  • 1970-01-01
  • 1970-01-01
  • 2021-01-27
  • 2020-01-12
相关资源
最近更新 更多