匹配所有字母前的数字运算--大工程~

huang6894

LZ生物信息的，想处理这样一个bam文本：

1. FCD2FEHACXX:8:2301:1875:46038#79_118    99      chr1    112365  0       16M1D70M2S      =       112469  192     CTTCTCCAACCAAATTAGGAACAATTATATGGCCACATAGTATCGAATCAAGTTTATAATTTTAAAATAATTGGGAGATTTTGTTGTT        ?A@@@A@@@A@@AA@ABA@AAAA@A?@@@A$AA@AA@@?A@?@A>@A@AAAA=AA?A@AAAAAAAAA@?A@ABBB@C=@AAAC@AB?A
   
2. FCD2FEHACXX:8:1104:6915:94857#79_118    147     chr1    73823   37      4S84M   =       73710   -201    AGCTGTAGAGAGAGGCACCCTTTTTTTTTTTTTTAATTATACTTTAAGTTTTAGGGTACATGTGCACCTTGTGCAGGTTAGTTACATA        ??A?>=<879<>=6/2=AA@AAAAAAAAAAAA?A@A?@@?BAA?BAAAAA?A@@A?@A@@B@AA?ABAAA@BB@@AA@AAA@?@??>@
   
3. FCD2FEHACXX:8:2307:11230:38678#79_118   14      chr1    37911759        29      44M4I13M4D8M1I13M5S     =       37911601        -240    CCCTCTGCCTCTCCTTCTCTTTACCTCCCTCCCTCTTTCCCTCTTTCCCTCTCTCTCTCCCTCCCTCTCTTCCTCTCTCTCTCCCTCT        ABA??<@AA?C@A=??CABA@<1AA=AAA>BABAA@?@BBCAC@@=A@A@?=?@A@B=AAA=@AB@C@AA@AB@A?A@A@A@A@@>AA
   

复制代码

每一行都有上面11列，列与列之间以tab分隔，现在我想以第六列（16M1D70M2S\88M\44M4I13M4D8M1I13M5S)这种数字计算“S”之前所有数字，其中M无意义，D表示忽略，I表示加，也就是说像16M1D70M2S这样的就是16 + 70 =86；4S84M表示0；44M4I13M4D8M1I13M5S表示：44 + 4 + 13 + 8  + 1 + 13 =83；得到这个数字（以A表示）之后输出第10列和11列从第A个开始S个字符，比如像16M1D70M2S这样的就是输出第87个开始2个字符TT和?A   ；4S84M表示0开始4个字符AGCT 和??A?；44M4I13M4D8M1I13M5S输出84开始5个字符CCTCT和@@>AA；（注意，匹配的字符必须大于1，小于或者等于1的行忽略） ;判断第9列如果是正数的话，N=1;第9列如果是负数的话，N=2;
最后输出的第一个文件结果是：
@第一列_N\n第10列匹配的字符\n+\n第11列匹配的字符\n\n:

1. @FCD2FEHACXX:8:2301:1875:46038#79_118_1   
   
2. TT
   
3. +
   
4. ?A
   
5. 
   
6. FCD2FEHACXX:8:1104:6915:94857#79_118_2
   
7. AGCT
   
8. +
   
9. ??A?
   
10. 
    
11. FCD2FEHACXX:8:2307:11230:38678#79_118_2
    
12. CCTCT
    
13. +
    
14. @@>AA
    
15.   

复制代码

第二个结果文件还是以第六列（16M1D70M2S\88M\44M4I13M4D8M1I13M5S)这种数字计算，这一次“S”之前所有数字忽略I，D表示加，最后是像16M1D70M2S这样的就是16 +1 + 70 =87；4S84M表示0；44M4I13M4D8M1I13M5S表示：44 + 13 +4 + 8 + 13 =82；得到这个数字（以A表示）之后输出第四列的数字加上这个数字之后的结果：
16M1D70M2S：112365+87=112452；
4S84M：73823+0=73823；
44M4I13M4D8M1I13M5S：37911759+82=37911841；
最后输出：

1. @FCD2FEHACXX:8:2301:1875:46038#79_118_1   16M1D70M2S  112452
   
2. FCD2FEHACXX:8:1104:6915:94857#79_118_2  4S84M  73823；
   
3. FCD2FEHACXX:8:2307:11230:38678#79_118_2  44M4I13M4D8M1I13M5S 37911841  

复制代码

---------------------------------------------------------
补充说明：
如果遇到30S23M35S这种多个S的情况，需要计算多次，一次是30S58M的，一次是53M35S的~
——————————————————————————————————

这种大工程，用perl做了一次，600多行啊！！！600多行啊啊啊啊~
有没有高手尝试一下用shell做一下，教教我吧~谢谢咯~

yestreenstars

试一下：

1. awk '{s=0;t=$6;N=$9>0?1:2;gsub(/[0-9]*[DIM]/,"",t);l=split(t,a,"S");if(l>2){if(a[1]>10)print $1"_"N"\n"substr($10,1,a[1])"\n+\n"substr($11,1,a[1])"\n";if(a[2]>10)print $1"_"N"\n"substr($10,length($10)-a[2]+1)"\n+\n"substr($11,length($11)-a[2]+1)"\n"}else{if(a[1]>10){gsub(/[0-9]*D|[0-9]*S.*/,"",$6);split($6,b,"[MI]");for(i in b)s+=b[i];print $1"_"N"\n"substr($10,s+1,a[1])"\n+\n"substr($11,s+1,a[1])"\n"}}}'

复制代码

代码有待优化~

yestreenstars

信息量好大~

lkk_super

这种大工程，用perl做了一次，600多行啊！！！600多行啊啊啊啊~
然后有啥问题? 既然 perl 都搞定了为啥还要搞成 shell ，淡淡忧伤啊

关阴月飞

表示理解不了.....

huang6894

回复 3# lkk_super


    为了学习。。。

huang6894

回复 4# 关阴月飞


    关大大~是你出手的时候了~

huang6894

回复 2# yestreenstars
快助我一臂之力！！！


   

blackold

还是小工程来钱快。

huang6894

回复 8# blackold


    米有办法~最高只能悬赏20~其实就为了引起大家注意~