多个文本的信息匹配问题求助

huang6894

各位大大:
我想咨询一下,如果我现在有以下三个文件
__________________________________________________________________________
all.txt:
refseq1 refGene mRNA    1   2   .       +      .       NR_024540       WASH5P
refseq1 refGene exon    2   10   .       +       .       NR_024540       WASH5P  E1     
refseq1 refGene intr    4   6   .       +       .       NR_024540       WASH5P  I1
refseq1 refGene exon    9   11   .       +      .       NR_024540       WASH5P  E2     
refseq2 refGene exon    1  6  .       +       .       NM_001005221    OR4F29  
refseq2 refGene exon    3  6  .       +       .       NM_001005221    OR4F29  E1     
refseq2 refGene exon    7  9  .       +       .       NM_001005221    OR4F29  E1      
refseq2 refGene mRNA    1  3  .       +       .       NM_001005224    OR4F3
refseq2 refGene exon    6  8  .       +       .       NM_001005224    OR4F3   E1     

nm.txt
NR_024540
NM_001005221

seq.txt
>refseq1 chr1 + 1000000 [1,1000000] ""
AGCTCGGTCCCCCCCCCCCCCCCTTTTTTT

>refseq2 chr2 + 1000000 [1,1000000] ""
AGCCCCCCCTCGGTCCCCCCCCCTTTTTTT
_________________________________________________________________________

我想得到的是nm.txt的字符串对应all.txt的第9列且all.txt的第十一列以E开始/^E/(不确定是否每一行都有第11列)的对应行,匹配后根据第一列的信息匹配seq.txt上的序列信息,获取第四列到第五列数值间的碱基,如果同一NM号且碱基序列间有重叠,合并输出.

也就是说最后得到的结果应该是:
_______________________________________________________________

>NR_024540_E1   WASH5P  refseq1 +       2  10  9 11  
GCTCGGTCCC

> NM_001005221_E1  OR4F29   refseq2 +       3  6
CCCC

> NM_001005221_E2  OR4F29   refseq2 +       7  9
CCC

可是我花费很多时间都得不出结果,怎么办呀??求助求助!

liion631818

回复 1# huang6894

1. 
   
2. declare -A seqArr
   
3. declare -A all
   
4. 
   
5. while read line; do
   
6.     if [[ "$line" == "" ]]; then
   
7.                 continue        
   
8.         fi
   
9.     if [[ "$line" =~ ">ref" ]]; then
   
10.         ref=${line%% *}
    
11.         ref=${ref:1}
    
12.         continue
    
13.     fi
    
14.     seqArr["$ref"]="$line"
    
15. done <seq.txt
    
16. 
    
17. while read line; do
    
18.         arr=($line)
    
19.         if [[ ${arr[10]} =~ "E" ]]; then
    
20.                 ref=${arr[0]}
    
21.                 nm=${arr[8]}
    
22.                 if `grep $ref seq.txt &>/dev/null`; then
    
23.                         if `grep $nm nm.txt &>/dev/null`; then
    
24.                
    
25.                                 echo \>${nm}_${arr[10]} ${arr[9]} $ref "+" ${arr[3]} ${arr[4]}
    
26.                                 dna=${seqArr["$ref"]}
    
27.                                 echo ${dna:(( ${arr[3]}-1)):(( ${arr[4]}-${arr[3]}+1 ))}
    
28.                         fi       
    
29.                 fi
    
30.         fi
    
31. done <all.txt
    

复制代码

huang6894

回复 2# liion631818


    大大，不知道为什么：
                                dna=${seqArr["$ref"]}
                                echo ${dna ${arr[3]}-1)) ${arr[4]}-${arr[3]}+1 ))}
这段好像echo不出来

rdcwayx

回复 3# huang6894


    在bash里运行。

1. bash script.sh

复制代码

rdcwayx

回复 2# liion631818


    没有实现这条: 如果同一NM号且碱基序列间有重叠,合并输出.

huang6894

回复 5# rdcwayx


    是呢，怎么办呢？

liion631818

回复 6# huang6894

1. 
   
2. #!/bin/bash
   
3. 
   
4. declare -A seqArr
   
5. declare -A all
   
6. 
   
7. while read line; do
   
8.     if [[ "$line" == "" ]]; then
   
9.                 continue         
   
10.         fi
    
11.     if [[ "$line" =~ ">ref" ]]; then
    
12.         ref=${line%% *}
    
13.         ref=${ref:1}
    
14.         continue
    
15.     fi
    
16.     seqArr["$ref"]="$line"
    
17. done <seq.txt
    
18. 
    
19. while read line; do
    
20. 
    
21.         arr=($line)
    
22.         if [[ ${arr[10]} =~ "E" ]]; then
    
23.                 ref=${arr[0]}
    
24.                 nm=${arr[8]}
    
25.                 if `grep $ref seq.txt &>/dev/null`; then
    
26.                         if `grep $nm nm.txt &>/dev/null`; then
    
27.                                 all[\>${nm}_${arr[10]} ${arr[9]} $ref]="${all[\>${nm}_${arr[10]} ${arr[9]} $ref]} ${arr[3]} ${arr[4]}"
    
28.                         fi        
    
29.                 fi
    
30.         fi
    
31. done <all.txt
    
32. 
    
33. for i in "${!all[@]}"; do
    
34.     set -- $i
    
35.     ref=$3
    
36.     dna=${seqArr["$ref"]}
    
37.     arr1=(`echo ${all[$i]} | xargs -n1 | sort -n` )
    
38.     (( last=${#arr1[@]}-1 ))
    
39.     echo $i" +"${all[$i]}
    
40.     echo ${dna:(( ${arr1[0]}-1)):(( ${arr1[$last]}-${arr1[0]}+1 ))}
    
41. done
    

复制代码

huang6894

回复 7# liion631818


  谢谢大大，我还有个问题：
  我在执行这个脚本的时候，提示
——————————————————————————————————————————————————————————————————————

1. test.sh: line 3: declare: -A: invalid option
   
2. declare: usage: declare [-afFirtx] [-p] [name[=value] ...]
   
3. test.sh: line 4: declare: -A: invalid option
   
4. declare: usage: declare [-afFirtx] [-p] [name[=value] ...]
   
5. test.sh: line 26: >NR_024540_E1 WASH5P refseq1: syntax error: operand expected (error token is ">NR_024540_E1 WASH5P refseq1")

复制代码

——————————————————————————————————————————————————————————————————————
  怎么解决呢？

liion631818

回复 8# huang6894


    你的bash版本太旧了，不支持关联数组

1. 
   
2. $ bash -version
   
3. GNU bash, version 4.2.24(1)-release (i686-pc-linux-gnu)
   
4. Copyright (C) 2011 Free Software Foundation, Inc.
   
5. License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
   
6. 
   
7. This is free software; you are free to change and redistribute it.
   
8. There is NO WARRANTY, to the extent permitted by law.
   
9. 
   

复制代码

huang6894

回复 9# liion631818

--------------------------------------------------------------------
$ bash -version
GNU bash, version 3.2.25(1)-release (x86_64-redhat-linux-gnu)
Copyright (C) 2005 Free Software Foundation, Inc.
-----------------------------------------------------------------
declare -a seqArr
declare -a all
是不是也可以定义数组呢？
————————————————————————————————
我不理解的是：
test.sh: line 26: >NR_024540_E1 WASH5P refseq1: syntax error: operand expected (error token is ">NR_024540_E1 WASH5P refseq1")
这条信息