Chinaunix

标题: 匹配所有字母前的数字运算--大工程~ [打印本页]
作者: huang6894    时间: 2013-11-15 11:57
标题: 匹配所有字母前的数字运算--大工程~
本帖最后由 huang6894 于 2013-11-15 13:48 编辑

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

  6. FCD2FEHACXX:8:1104:6915:94857#79_118_2
  7. AGCT
  8. +
  9. ??A?

  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    时间: 2013-11-15 11:57
试一下：
  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    时间: 2013-11-15 12:04
信息量好大~
作者: lkk_super    时间: 2013-11-15 12:23
这种大工程，用perl做了一次，600多行啊！！！600多行啊啊啊啊~
然后有啥问题? 既然 perl 都搞定了为啥还要搞成 shell ，淡淡忧伤啊
作者: 关阴月飞    时间: 2013-11-15 12:39
表示理解不了.....
作者: huang6894    时间: 2013-11-15 12:58
回复 3# lkk_super


    为了学习。。。
作者: huang6894    时间: 2013-11-15 12:58
回复 4# 关阴月飞


    关大大~是你出手的时候了~
作者: huang6894    时间: 2013-11-15 12:59
回复 2# yestreenstars
快助我一臂之力！！！


   
作者: blackold    时间: 2013-11-15 13:11
还是小工程来钱快。
作者: huang6894    时间: 2013-11-15 13:19
回复 8# blackold


    米有办法~最高只能悬赏20~其实就为了引起大家注意~
作者: klainogn    时间: 2013-11-15 13:27
本帖最后由 klainogn 于 2013-11-15 15:58 编辑
  1. awk '{match($6,/[0-9]+S/);len=substr($6,RSTART,RLENGTH)+0;str=$1"_"gsub("-","",$9)+1}NR==FNR{gsub(/[0-9]*D|[0-9]*S86M|[0-9]*S/,"",$6);split($6,a,/[A-Z]/);for(i in a) start+=a[i];print str"\n"substr($10,start,len)"\n+\n"substr($11,start,len)"\n";start=0;next}{str2=$6;gsub(/[0-9]*I|[0-9]*S86M|[0-9]*S/,"",$6);split($6,a,/[A-Z]/);for(i in a) start+=a[i];print str,str2,$4+start;start=0}' file file

  3. FCD2FEHACXX:8:2301:1875:46038#79_118_1
  4. GT
  5. +
  6. B?

  8. FCD2FEHACXX:8:1104:6915:94857#79_118_2
  9. AG
  10. +
  11. ??

  13. FCD2FEHACXX:8:2307:11230:38678#79_118_2
  14. CCCTC
  15. +
  16. A@@>A

  18. FCD2FEHACXX:8:2301:1875:46038#79_118_1 16M1D70M2S 112452
  19. FCD2FEHACXX:8:1104:6915:94857#79_118_2 2S86M 73823
  20. FCD2FEHACXX:8:2307:11230:38678#79_118_2 44M4I13M4D8M1I13M5S 37911841
复制代码
作者: yestreenstars    时间: 2013-11-15 13:39
楼主你给的数据第10列和第11列的长度不同的？
  1. [root@localhost ~]# awk '{print length($10),length($11)}' i
  2. 88 87
  3. 88 87
  4. 88 88
  5. [root@localhost ~]#
复制代码
作者: huang6894    时间: 2013-11-15 13:40
回复 10# klainogn


    大大，gsub(/[0-9]*D|[0-9]*S86M|[0-9]*S/,"",$6);这一句话。。。好像不具有普适性呀，好像是不是只适合当前测试数据呀？
作者: huang6894    时间: 2013-11-15 13:41
回复 11# yestreenstars


    sorry~可能漏写了~反正是很怪异的数字~两列是一样长的~
作者: yestreenstars    时间: 2013-11-15 13:42
第二行，明明是2S，为什么你的结果里确实4S的结果？
作者: yestreenstars    时间: 2013-11-15 13:43
回复 13# huang6894
那你补上吧，难怪我的测试结果不对~

   
作者: klainogn    时间: 2013-11-15 13:50
回复 12# huang6894
补充说明：
如果遇到30S23M35S这种多个S的情况，需要计算多次，一次是30S58M的，一次是53M35S的~
这个点还需要 星辰 哥发挥一下


   
作者: huang6894    时间: 2013-11-15 13:50
本帖最后由 huang6894 于 2013-11-15 14:09 编辑

回复 14# yestreenstars


    对不起对不起，实在太不仔细了~已经改过来了~谢谢大大~~
  1. [root@huang]$ cat c.txt
  2. 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
  3. 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@?@??>@
  4. 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
  5. [root@huang]$ awk '{print length($10),length($11)}' c.txt
  6. 88 88
  7. 88 88
  8. 88 88
复制代码
作者: huang6894    时间: 2013-11-15 13:51
回复 15# yestreenstars


    谢谢大大，改过来了~谢谢谢谢
作者: yestreenstars    时间: 2013-11-15 13:54
回复 18# huang6894
还有我14楼问的问题呢~

   
作者: Herowinter    时间: 2013-11-15 13:55
压力好大，表示看懂需求要一段时间，
今天上班比较忙，只能拜读各位的代码了。
作者: huang6894    时间: 2013-11-15 14:09
回复 19# yestreenstars


    已修正~是4S的！对不起呀，测试数据没弄好~
作者: yestreenstars    时间: 2013-11-15 14:10
先发第一种结果，多个S的处理看得不是很明白：
  1. [root@localhost ~]# cat i
  2. 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
  3. 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@?@??>@
  4. 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
  5. [root@localhost ~]# awk '{s=0;n=gensub(/.*(^|[^0-9])([0-9]*)S.*/,"\\2",$6);gsub(/[0-9]*D|[0-9]*S.*/,"",$6);split($6,a,"[MI]");for(i in a)s+=a[i];N=$9>0?1:2;print $1"_"N"\n"substr($10,s+1,n)"\n+\n"substr($11,s+1,n)"\n"}' i
  6. FCD2FEHACXX:8:2301:1875:46038#79_118_1
  7. TT
  8. +
  9. ?A

  11. FCD2FEHACXX:8:1104:6915:94857#79_118_2
  12. AGCT
  13. +
  14. ??A?

  16. FCD2FEHACXX:8:2307:11230:38678#79_118_2
  17. CCTCT
  18. +
  19. @@>AA

  21. [root@localhost ~]#
复制代码
作者: yestreenstars    时间: 2013-11-15 14:34
关于你说的特殊情况（有多个S的）能否给个案例？
作者: huang6894    时间: 2013-11-15 14:34
回复 22# yestreenstars


    大大，有个问题~你的脚本处理测试数据是通过的~可是不知道为什么，和测试数据一样的实际文本却不能处理呢？

方便的话，帮我看看我的数据可以么?谢谢你了！！！
use.zip (148.74 KB, 下载次数: 14)
  1. awk '{s=0;n=gensub(/.*(^|[^0-9])([0-9]*)S.*/,"\\2",$6);gsub(/[0-9]*D|[0-9]*S.*/,"",$6);split($6,a,"[MI]");for(i in a)s+=a[i];N=$9>0?1:2;print $1"_"N"\n"substr($10,s+1,n)"\n+\n"substr($11,s+1,n)"\n"}' c.txt
  2. FCD2FEHACXX:8:2301:1875:46038#79_118_1
  3. TT
  4. +
  5. ?A

  7. FCD2FEHACXX:8:1104:6915:94857#79_118_2
  8. AG
  9. +
  10. ??

  12. FCD2FEHACXX:8:2307:11230:38678#79_118_2
  13. CCTCT
  14. +
  15. @@>AA

  17. awk '{s=0;n=gensub(/.*(^|[^0-9])([0-9]*)S.*/,"\\2",$6);gsub(/[0-9]*D|[0-9]*S.*/,"",$6);split($6,a,"[MI]");for(i in a)s+=a[i];N=$9>0?1:2;print $1"_"N"\n"substr($10,s+1,n)"\n+\n"substr($11,s+1,n)"\n"}' use.txt |le
  18. FCD2FEHACXX:8:1203:10339:145301#79_118_2

  20. +


  23. FCD2FEHACXX:8:1202:16979:45106#79_118_2

  25. +


  28. FCD2FEHACXX:8:1205:12796:100011#79_118_1

  30. +


  33. FCD2FEHACXX:8:1201:19531:81831#79_118_2

  35. +



复制代码
作者: huang6894    时间: 2013-11-15 14:43
本帖最后由 huang6894 于 2013-11-15 14:52 编辑

回复 23# yestreenstars


    比如：
  1. FCD2FEHACXX:8:1303:20504:42704#79_118   99      chr1    15317774        29      20S33M35S  =       15317807        119     CCGAGAGGCGGAGGTTGCAGTGAGCCTAGGTCGCACCACTGCACTCCAGCCTGGTGACAAAGAGAGACTCTGTGTCAAAAAAAAGAAG        ?@<AA@BAA>A@AA>@AA@B?A@AA@@ABA?B>C@CA@CACC@BACA@CCAACB<AAC?AABABAA@CACAC?C?@@AAAAAAAC@AC        BD:Z:LLLLMMMMLLLLMMLLLLLMLMLMLLLLLMLLLLLLLMLLMLLLLLMMMLLLMMLMLLLLGLMMMMMLLLLMMMMLLLGGGGGGLMLL   MD:Z:11G17A3    PG:Z:MarkDuplicates.4   RG:Z:FAMILY_all XG:i:0  BI:Z:OOOOOOOOOOOOOONOOOOOOOOOOOOOOONOOOOOOOOOOOOOONOOOOOOOONOOOOONNOOOOOOONOOOOOOOONNNNNNNOON   AM:i:29 NM:i:2  SM:i:29 XM:i:2  XO:i:0  MQ:i:29 XT:A:M
复制代码
这样的
希望可以输出：

  2. FCD2FEHACXX:8:1303:20504:42704#79_118_1
  3. CCGAGAGGCGGAGGTTGCAG
  4. +
  5. ?@<AA@BAA>A@AA>@AA@B

  7. FCD2FEHACXX:8:1303:20504:42704#79_118_1
  8. GTGACAAAGAGAGACTCTGTGTCAAAAAAAAGAAG
  9. +
  10. B<AAC?AABABAA@CACAC?C?@@AAAAAAAC@AC
复制代码
作者: yestreenstars    时间: 2013-11-15 16:17
  1. awk '{s=0;t=$6;N=$9>0?1:2;gsub(/[0-9]*[DIM]/,"",t);l=split(t,a,"S");if(l>2){print $1"_"N"\n"substr($10,1,a[1])"\n+\n"substr($11,1,a[1])"\n\n"$1"_"N"\n"substr($10,length($10)-a[2]+1)"\n+\n"substr($11,length($11)-a[2]+1)}else{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"}}'
复制代码
不知道有没有bug，试一下吧~
作者: huang6894    时间: 2013-11-15 16:30
回复 26# yestreenstars


    星辰大大真是神一样的存在呀！！！！
    大大，我还想问一下，如果我不想处理那些匹配出来的第10、11列字符数小于10个的，如果小于10个不输出这一行，怎么改？？
作者: yestreenstars    时间: 2013-11-15 16:39
你的第二种结果文件如果有2个S的又怎么处理？
作者: damcool    时间: 2013-11-15 16:45
楼主对要求的解释是在看不懂。
作者: huang6894    时间: 2013-11-15 16:48
回复 29# yestreenstars


比如：
  1. FCD2FEHACXX:8:1303:20504:42704#79_118   99      chr1    15317774        29      20S33M35S  =       15317807        119     CCGAGAGGCGGAGGTTGCAGTGAGCCTAGGTCGCACCACTGCACTCCAGCCTGGTGACAAAGAGAGACTCTGTGTCAAAAAAAAGAAG        ?@<AA@BAA>A@AA>@AA@B?A@AA@@ABA?B>C@CA@CACC@BACA@CCAACB<AAC?AABABAA@CACAC?C?@@AAAAAAAC@AC        BD:Z:LLLLMMMMLLLLMMLLLLLMLMLMLLLLLMLLLLLLLMLLMLLLLLMMMLLLMMLMLLLLGLMMMMMLLLLMMMMLLLGGGGGGLMLL   MD:Z:11G17A3    PG:Z:MarkDuplicates.4   RG:Z:FAMILY_all XG:i:0  BI:Z:OOOOOOOOOOOOOONOOOOOOOOOOOOOOONOOOOOOOOOOOOOONOOOOOOOONOOOOONNOOOOOOONOOOOOOOONNNNNNNOON   AM:i:29 NM:i:2  SM:i:29 XM:i:2  XO:i:0  MQ:i:29 XT:A:M
复制代码
我希望得到：
20S33M35S：15317774+0=15317774；
20M33M35S：15317774+53=15317827；

  2. @FCD2FEHACXX:8:1303:20504:42704#79_118 _1   20S33M35S  15317774
  3. @FCD2FEHACXX:8:1303:20504:42704#79_118 _1   20S33M35S  15317827
复制代码
作者: yestreenstars    时间: 2013-11-15 17:11
第二题：
  1. awk '{N=$9>0?1:2;s=0;t=$6;gsub(/[0-9]*[^S0-9]/,"",t);l=split(t,a,"S");if(l<3){t=$6;gsub(/[0-9]*(I|S.*)/,"",t);split(t,b,"[DM]");for(i in b)s+=b[i];print $1"_"N,$6,$4+s}else{print $1"_"N,$6,$4;t=$6;gsub(/[0-9]*(I|S$)/,"",t);split(t,b,"[DMS]");for(i in b)s+=b[i];print $1"_"N,$6,$4+s}}'
复制代码
作者: huang6894    时间: 2013-11-15 17:26
回复 32# yestreenstars


    大哥大哥！ 还是那个问题：如果我不想处理那些匹配出来的第10、11列字符数小于10个的，如果小于10个不输出这一行，怎么改呀？？
   
   谢谢大大~
   ====================
   深深的挫败感啊~我600多行的perl~
作者: yestreenstars    时间: 2013-11-15 17:35
回复 33# huang6894
第二题没有用到第10、11列吧？如果你问的是第一题，请看28楼~

   
作者: huang6894    时间: 2013-11-15 17:41
回复 34# yestreenstars


    嗯嗯，大大，是这样的，我就是说因为第一题运算的结果是用于匹配第10、11列嘛~然后我就想这个数如果小于10的话，我就不要这一行了~就是S前面的数字运算结果必须大于10，比较麻烦的是这一次忽略I，D表示加
作者: yestreenstars    时间: 2013-11-15 17:50
回复 35# huang6894
你的意思是要把两种结果合并在一起？如果是请给个案例，看你要什么样的格式~

   
作者: huang6894    时间: 2013-11-15 17:55
回复 36# yestreenstars


    sorry,大大，我没说清楚，结果我是想生成两个文件的，第一个文件你已经帮我完美解决了~
    第二个文件的话是这个意思：
    比如：
  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
复制代码
这个还是以第六列（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；#这个时候S前面的数字运算结果是87>10,所以输出这行
4S84M：73823+0=73823；#这个时候S前面的数字运算结果是0<10,所以不输出这行
44M4I13M4D8M1I13M5S：37911759+82=37911841；#这个时候S前面的数字运算结果是82>10,所以输出这行
最后输出：
@FCD2FEHACXX:8:2301:1875:46038#79_118_1   16M1D70M2S  112452
@FCD2FEHACXX:8:2307:11230:38678#79_118_2  44M4I13M4D8M1I13M5S 37911841  
可以吗？
作者: damcool    时间: 2013-11-15 17:56
本帖最后由 damcool 于 2013-11-15 21:33 编辑

问题一，可处理多S：
  1. awk '{s=$6;gsub(/[0-9]+D|[0-9]+I|[0-9]+M|S$/,"",s);m=$6;gsub(/[0-9]+S|[0-9]+D/,"",m);n=split(m,a,"[MI]");t=0;if (match(m,"I")) for (i=1;i<=n;i++) t+=a[i];n=split(s,a,"S");for (i=1;i<=n;i++) if (length(substr($10,t,a[i]))>1&& length(substr($11,t,a[i]))>1) printf "%s_%d\n%s\n%s\n%s\n\n",$1,($9>0?1:2),substr($10,t,a[i]),"+",substr($11,t,a[i])}' data
复制代码
问题二，可指定字符串长度(通过strlen指定)：
  1. awk -v strlen=10 '{m=$6;gsub(/D/,"I",m);gsub(/[0-9]+S/,"",m);n=split(m,a,"[MI]");t=0;if (match(m,"I")) for (i=1;i<=n;i++) t+=a[i];if (t>strlen) printf "%s_%d\t%s\t%d\n\n",$1,($9>0?1:2),$6,$4+t}' data
复制代码
作者: yestreenstars    时间: 2013-11-15 18:06
回复 37# huang6894
试一下：
  1. awk '{N=$9>0?1:2;s=0;t=$6;gsub(/[0-9]*[^S0-9]/,"",t);l=split(t,a,"S");if(l<3){t=$6;gsub(/[0-9]*(I|S.*)/,"",t);split(t,b,"[DM]");for(i in b)s+=b[i];if(s>=10)print $1"_"N,$6,$4+s}else{print $1"_"N,$6,$4;t=$6;gsub(/[0-9]*(I|S$)/,"",t);split(t,b,"[DMS]");for(i in b)s+=b[i];if(s>=10)print $1"_"N,$6,$4+s}}'
复制代码
作者: huang6894    时间: 2013-11-15 18:42
回复 38# damcool


    牛人！
作者: huang6894    时间: 2013-11-15 18:42
回复 39# yestreenstars


    谢谢您！
作者: damcool    时间: 2013-11-15 21:42
huang6894 发表于 2013-11-15 18:42
回复 38# damcool

希望不是为了研究转基因东东来害我们。。。。
作者: yestreenstars    时间: 2013-11-16 02:05
用perl琢磨了一下，我对perl不是很熟，只能写成这样了，第一题：
  1. #!/usr/bin/perl
  2. use 5.010;
  3. open(FILE,"use.txt") or die;
  4. $len = 10;
  5. while (<FILE>) {
  6.         split;
  7.         $sum = 0;
  8.         $N = $_[8] > 0 ? 1 : 2;
  9.         while ($_[5] =~ /(\d+)S/g) {
  10.                 push(@numbers,$1);
  11.         }
  12.         $length = @numbers;
  13.         if ($length > 1) {
  14.                 if ($numbers[0] >= $len) {
  15.                         printf "%s_%s\n%s\n+\n%s\n\n",$_[0],$N,substr($_[9],0,$numbers[0]),substr($_[10],0,$numbers[0]);
  16.                 }
  17.                 if ($numbers[1] >= $len) {
  18.                         printf "%s_%s\n%s\n+\n%s\n\n",$_[0],$N,substr($_[9],length($_[9])-$numbers[1]),substr($_[10],length($_[10])-$numbers[1]);
  19.                 }
  20.         } else {
  21.                 if ($numbers[0] >= $len) {
  22.                         if ($_[5] =~ /S$/) {
  23.                                 while ($_[5] =~ /(\d+)[MI]/g) {
  24.                                         $sum += $1;
  25.                                 }
  26.                         }
  27.                         printf "%s_%s\n%s\n+\n%s\n\n",$_[0],$N,substr($_[9],$sum,$numbers[0]),substr($_[10],$sum,$numbers[0]);
  28.                 }
  29.         }
  30.         undef @numbers;
  31. }
复制代码
作者: yestreenstars    时间: 2013-11-16 02:32
优化了一下第一题的awk：
  1. #!/usr/bin/awk
  2. {
  3.         s=0;t=$6;N=$9>0?1:2;
  4.         gsub(/[0-9]*[DMI]/,"",t);
  5.         l=split(t,a,"S");
  6.         if(l>2){
  7.                 if(a[1]>=10)printf "%s_%s\n%s\n+\n%s\n\n",$1,N,substr($10,1,a[1]),substr($11,1,a[1]);
  8.                 if(a[2]>=10)printf "%s_%s\n%s\n+\n%s\n\n",$1,N,substr($10,length($10)-a[2]+1),substr($11,length($11)-a[2]+1)
  9.         } else {
  10.                 if(a[1]>=10){
  11.                         gsub(/[0-9]*(D|S.*)/,"",$6);
  12.                         split($6,b,"[MI]");
  13.                         for(i in b)s+=b[i];
  14.                         printf "%s_%s\n%s\n+\n%s\n\n",$1,N,substr($10,s+1,a[1]),substr($11,s+1,a[1])
  15.                 }
  16.         }
  17. }
复制代码
作者: huang6894    时间: 2013-11-18 09:30
@damcool谢谢您~@yestreenstars由衷感谢大大的帮助~
作者: pitonas    时间: 2013-11-18 14:40
damcool 发表于 2013-11-15 14:42
希望不是为了研究转基因东东来害我们。。。。

希望不是{:2_179:}



欢迎光临 Chinaunix (http://bbs.chinaunix.net/) Powered by Discuz! X3.2