求教 范围匹配

catloverat99

file1最后一列 rel_pos的值是否落在另一个文件file2中的区间(例如第一行 267为起点 346为钟点）。
ENSRNOG000000xx1是不同的ID，是两个文件都含有的. 注：一定要ID先匹配再考虑范围
最后输出 2种
1. 以file2为主，如果file1的rel_pos的值落在另一个文件file2中的区间，在file2添加一列，填yes，否则no
2.以file1 为主，如果file1的rel_pos的值落在另一个文件file2中的区间，在file1添加一列，填yes，否则no。

感谢之前大神的代码，$ awk 'FNR==NR{a[$NF]=$0;next}{YN="No";for(n in a)if($2<=+n&&+n<=$3){YN="Yes";break}print $0,YN}' file1 file2
结果似乎不太对，请问是不是这里忽略了id的匹配? 如何把id给匹配上

file1
value.V1        value.V2        value.V3        value.V4        value.V5        value.V6        value.V7        value.V8        value.V9        chromosome        start        end        dmr_length        rel_pos
chr1        17618115        T        C        CC        CC        TT        TT        ENSRNOG00000013436_451        chr1        17617848        17618216        369        268
chr1        17618162        A        C        CC        CC        AA        AA        ENSRNOG00000013436_451        chr1        17617848        17618216        369        315
chr1        24197584        A        G        GG        GG        AG        AG        ENSRNOG00000016371_102        chr1        24197021        24198790        1770        564
chr1        24197597        C        T        TT        TT        CT        CT        ENSRNOG00000016371_102        chr1        24197021        24198790        1770        577
chr1        24197733        G        C        CC        CC        GC        GG        ENSRNOG00000016371_102        chr1        24197021        24198790        1770        713
chr1        24197750        G        A        AA        AA        GG        GG        ENSRNOG00000016371_102        chr1        24197021        24198790        1770        730

file2
V$GLI_Q3 267 346 g2_ENSRNOG00000013436_451
V$GKLF_Q4 345 351 g2_ENSRNOG00000013436_451
V$IK_Q5_01 345 351 g2_ENSRNOG00000013436_451
V$TTF1_Q5_01 347 353 g2_ENSRNOG00000013436_451
V$COE1_Q6 347 360 g2_ENSRNOG00000013436_451
V$TBX5_01 13 24 g2_ENSRNOG00000016371_102
V$RBPJK_01 27 37 g2_ENSRNOG00000016371_102
V$DR4_Q2 34 50 g2_ENSRNOG00000016371_102
V$BBX_04 117 133 g2_ENSRNOG00000016371_102
V$BBX_04 118 134 g2_ENSRNOG00000016371_102

catloverat99

没有大神来帮忙看看么。。。

tc1989tc

1. cat awk.awk
   
2. #!/usr/bin/awk -f
   
3. 
   
4. # Recoder per line base on file2's id
   
5. NR==FNR {
   
6.     file1num ++
   
7.     file1[FNR]=$0
   
8.     if ($1 ~ "value")
   
9.         next
   
10. 
    
11.     if (id1lines[$9]) {
    
12.         id1lines[$9]=id1lines[id]""FNR"-"
    
13.     } else {
    
14.         id1lines[$9]=FNR"-"
    
15.     }
    
16.     next   
    
17. }
    
18. {
    
19.     file2[FNR]=$0
    
20.     file2num ++
    
21.     split($4,a,"_")
    
22.     id=a[2]"_"a[3]
    
23.     if (id2lines[id]) {
    
24.         id2lines[id]=id2lines[id]""FNR"-"
    
25.     } else {
    
26.         id2lines[id]=FNR"-"
    
27.     }
    
28. }
    
29. END {
    
30.     addyesnotofile1(file1,file2)
    
31.     addyesnotofile2(file1,file2)   
    
32. }
    
33. function addyesnotofile1(file1,file2) {
    
34.     len = file1num
    
35.     for (i=1; i <= len; i++) {
    
36.         if (file1[i] ~ "value") {
    
37.             print file1[i] >> "file1.tmp"
    
38.             continue
    
39.         }
    
40.         # found all lines in file2 base on keyid
    
41.         split(file1[i],a)
    
42.         keyid = a[9]
    
43.         repos = a[14]
    
44. 
    
45.         # check if keyid exist in file2
    
46.         if (! id2lines[keyid]) {
    
47.             print file1[i]," NO" >> "file1.tmp"
    
48.             continue
    
49.         }
    
50.         found=0
    
51.         # check if range fit
    
52.         split(id2lines[keyid],a,"-")
    
53.         for (j in a) {
    
54.             split(file2[j],b)
    
55.             min = b[2]
    
56.             max = b[3]
    
57.             if ((min <= repos) && (repos <= max)) {
    
58.                 print file1[i]," YES" >> "file1.tmp"
    
59.                 found = 1
    
60.                 break
    
61.             }
    
62.         }
    
63. 
    
64.         if (found) {
    
65.             continue   
    
66.         }
    
67.         # not fit ranges
    
68.         print file1[i]," NO" >> "file1.tmp"
    
69.     }
    
70.         
    
71. }
    
72. 
    
73. function addyesnotofile2(file1,file2) {
    
74.     len = file2num
    
75.     for (i=1;i<=len;i++) {
    
76.         # give file2 keyid
    
77.         split(file2[i""],a)
    
78.         min = a[2]
    
79.         max = a[3]
    
80.         split(a[4],b,"_")
    
81.         keyid = b[2]"_"b[3]
    
82.         if (! id1lines[keyid]) {
    
83.             print file2[i]," NO" >> "file2.tmp"
    
84.             continue
    
85.         }
    
86.         found = 0
    
87.         # Travel all lines base on keyid
    
88.         split(id1lines[keyid],b,"-")
    
89.         for (j in b) {
    
90.             # give repos value
    
91.             split(file1[j],a)
    
92.             repos = a[14]
    
93.             if ((min <= repos) && (repos <= max)) {
    
94.                 found = 1
    
95.                 print file2[i]," YES" >> "file2.tmp"
    
96.                 break
    
97.             }
    
98.         }
    
99.         if (found) {
    
100.             continue
     
101.         }
     
102.         # not fit ranges
     
103.         print file2[i]," NO" >> "file2.tmp"
     
104.     }   
     
105. }

复制代码

catloverat99

谢谢！太感谢了，我先试试：）

haooooaaa

按行比较？

1. awk 'function get(h,n,w){
   
2.     split(h, u, n)
   
3.     return u[w]
   
4. }
   
5. function pp(q,r,v){
   
6.     for(i=0;++i<=length(q);){
   
7.         if(q[i]==r[i]){
   
8.             rel=get(F1A[i]," ",14)
   
9.             little=get(F2A[i]," ",2)
   
10.             large=get(F2A[i]," ",3)
    
11.             print (rel>little&&rel<large)? v[i]" YES":v[i]" NO"
    
12.         }else{
    
13.             print v[i]" NO"
    
14.         }
    
15.     }
    
16. }
    
17. NR==1&&/^ch/{f=NF>5?1:0}
    
18. NF>5&&/^ch/{F1[FNR]=get($(NF-5),"_",1);F1A[FNR]=$0}NF<5{F2[FNR]=get($NF,"_",2);F2A[FNR]=$0}
    
19. END{
    
20.     if(f){
    
21.         pp(F1,F2,F1A)
    
22.     }else{
    
23.         pp(F2,F1,F2A)
    
24.     }
    
25. }' a b

复制代码

catloverat99

回复 3# tc1989tc


刚才试着运行了一下,结果发现如果不同的 ENS ID 如果范围相同，结果也变成yes. 我需要先保证id相同的情况下，看file 1是否在file 2的范围。麻烦大神，再给看看。下面是运行结果。

file1的结果
value.V1        value.V2        value.V3        value.V4        value.V5        value.V6        value.V7        value.V8        value.V9        chromosome        start        end        dmr_length        rel_pos
chr1        17618115        T        C        CC        CC        TT        TT        ENSRNOG00000013436_451        chr1        17617848        17618216        369        268  YES
chr1        17618162        A        C        CC        CC        AA        AA        ENSRNOG00000013436_451        chr1        17617848        17618216        369        315  YES
chr1        24197584        A        G        GG        GG        AG        AG        ENSRNOG00000016371_102        chr1        24197021        24198790        1770        564  NO
chr1        24197597        C        T        TT        TT        CT        CT        ENSRNOG00000016371_102        chr1        24197021        24198790        1770        577  NO
chr1        24197733        G        C        CC        CC        GC        GG        ENSRNOG00000016371_102        chr1        24197021        24198790        1770        713  NO
chr1        24197750        G        A        AA        AA        GG        GG        ENSRNOG00000016371_102        chr1        24197021        24198790        1770        730  NO

file2的结果,倒数第二行也是yes. 这应该是匹配了范围，没有考虑id的结果
V$GLI_Q3 267 346 g2_ENSRNOG00000013436_451  YES
V$GKLF_Q4 345 351 g2_ENSRNOG00000013436_451  NO
V$IK_Q5_01 345 351 g2_ENSRNOG00000013436_451  NO
V$TTF1_Q5_01 347 353 g2_ENSRNOG00000013436_451  NO
V$COE1_Q6 347 360 g2_ENSRNOG00000013436_451  NO
V$TBX5_01 13 24 g2_ENSRNOG00000016371_102  NO
V$RBPJK_01 27 37 g2_ENSRNOG00000016371_102  NO
V$DR4_Q2 34 50 g2_ENSRNOG00000016371_102  NO
V$BBX_04 267 346 g2_ENSRNOG00000016371_102  YES
V$BBX_04 118 134 g2_ENSRNOG00000016371_102  NO

catloverat99

回复 5# haooooaaa


  结尾的a 和 b是要换成我的文件名file1 file2？抱歉我的水平太有限。

haooooaaa

回复 7# catloverat99


    是

catloverat99

回复 8# haooooaaa

结果怎么都是no?


V$GLI_Q3 267 346 g2_ENSRNOG00000013436_451 NO
V$GKLF_Q4 345 351 g2_ENSRNOG00000013436_451 NO
V$IK_Q5_01 345 351 g2_ENSRNOG00000013436_451 NO
V$TTF1_Q5_01 347 353 g2_ENSRNOG00000013436_451 NO
V$COE1_Q6 347 360 g2_ENSRNOG00000013436_451 NO
V$TBX5_01 13 24 g2_ENSRNOG00000016371_102 NO
V$RBPJK_01 27 37 g2_ENSRNOG00000016371_102 NO
V$DR4_Q2 34 50 g2_ENSRNOG00000016371_102 NO
V$BBX_04 117 133 g2_ENSRNOG00000016371_102 NO
V$BBX_04 267 346 g2_ENSRNOG00000016371_102 NO

haooooaaa

回复 9# catloverat99


    我测试，file2 第一个是YES，