Chinaunix

标题: 删除匹配的字符序列和其对应的质量序列 [打印本页]

作者: iamline    时间: 2013-09-24 17:56
标题: 删除匹配的字符序列和其对应的质量序列
本帖最后由 iamline 于 2013-09-24 17:59 编辑

dna测序的序列文件如下所示:

more seq.file
  1. @HWI-ST531R:212:C22J6ACXX:8:2314:4159:37409
  2. AAAAAAAAAATACTGTTTATGATGATATCATTGAAGCTGCGGGTGAATTAGGTGATAAAACGCTACTTCCTGTTTTAGATACTATGTTGTACAAGTTTGAT
  3. +
  4. C@CFFFFFHHHHHIIHGIIJGIJIBGEIIGJCGHIJIJGIJJJ8>EDDFCED6@ACCCDDD@@BDDDDDDCDDCDCCCCADDDDECDDDCDDDCCDCDCAC
  5. @HWI-ST531R:212:C22J6ACXX:8:2102:3421:60710
  6. AAAAAAAAACCATCCAAATCTGGATGGCTTTTCATAATTCTGAGAAATTAGCTGCGCTGGCGCACCGCTTCAAATAAGCAAATTCCGGTCGCAACCGAAAC
  7. +
  8. @@CFFFFFHHHHHJJJJJJJJJJJJJJJJJJJJIJIIJJIJJIJIJJJJIIIJJHHFFDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD
  9. @HWI-ST531R:212:C22J6ACXX:8:1303:19050:54218
  10. AAAAAAAAACCGCACCAGAATTCCCATCGATCCAAGGCCGTTGGTGTTCCTTTCCCCACAAACCGCCCCAACGTCTCAATCGAGACTCCCCCTCACAGCAG
  11. +
  12. ??<BDDDDFFFF08@FFBG>D8?FFFF9;;==5@@=(66:55?B(5;?BA;AA>AA<<;<?B9?5@08-<@BB<2<9(4+:28+25@:>8>B@@B<A>?##
  13. @HWI-ST531R:212:C22J6ACXX:8:2108:4215:95315 1:N:0:GAGTGG
  14. AAAAAAAAACCGGGCAATTGCCCGGTGCTGGAAGGGATCTCGTTATTCGGGGGAATAAGGTAGATGTGGATAGTCCAGCCAGTGCGTCAAATAGTGAGAGA
  15. +
  16. <<@DDDDDAH@DHG:FHFGIIGCGI@FHHIHIIFHE6?1;?@?==ACDCA@BB5<8>C@C33:@@A@>>:>:?::>4>9<<A::>9(55>A:>A4(+4@32
  17. @HWI-ST531R:212:C22J6ACXX:8:2208:5557:5216 1:N:0:GAGTGG
  18. AAAAAAAAACCTGAAAAAAACGGCCTGACGTGAATCAAGCAACTTTTTTCAGGTTTTGCCCGCTTAGTGCGGTAACAATCCTTTACTCAGTAATAATATTT
  19. +
  20. @@CFFFFFFDFFFEHGIGIHGHHGIIGEDEHFEHFBC@?DBA(;@CCDDC@CDCCCDDCDDD;DB7C>C@BB05?A>@CDDDDCCDDDDCCCCDDDDECDD
  21. @HWI-ST531R:212:C22J6ACXX:8:2214:1365:3705
  22. AAAAAAAAACCTGAAAAAAACGGCCTGACGTGAATCAAGCAATTTTTTTCAGGTTTTGCCCGCTTAGTGCGGTAACAATCCTTTACTCAGTAATAATATTT
  23. +
  24. @@@FDDFFF:DHFHGIGHEHEFGEAFHGIIHGCEDFFDE@A@>;ACDDD:>@C(>CD?@@CD@B@BDCA>;;0;B>CC@C<CCDCDAA>CACCCDDEC>CB
  25. @HWI-ST531R:212:C22J6ACXX:8:1107:2469:40337
  26. AAAAAAAAACCTGAAAAAAACGGCCTGACGTGAATCAAGCAATTTTTTTCAGGTTTTGCCCGCTTAGTGCGGTAACAATCCTTTACTCAGTAATAATATTT
  27. +
  28. CCCFFFFFHHHHHJJJJJJIJIJJJJJJJJHGHGFFEFEEECEEEEDDDCDDDACDDDDDDDDDDDDDDDDDBDDDDDDDDDDDDDDCDDDEDEEDEEEED
  29. @HWI-ST531R:212:C22J6ACXX:8:2211:2015:43926
  30. AAAAAAAAACCTGAAAAAAACGGCCTGACGTGAATCAAGCAATTTTTTTCAGGTTTTGCCCGCTTAGTGCGGTAACAATCCTTTACTCAGTAATAATATTT
  31. +
  32. CCCFFFFFFHHAFGGHDHIIIIIGHGGHGGHHAHEDFFFCAECCCEDDD>:CCCCDDDDDDD7@BDDCCDBB;B@CCDDCCDDDDDDCCCADDDEDEDEEE
复制代码
…………

即所谓的fastq格式,每4 行信息代表一条测序的序列,第一行和第三行表示序列ID(其中,第一行以“@” 开头而第三行以“+” 开头;第三行中 ID可以省略,但“+” 不能省略),第二行为碱基序列,第四行是第二行中序列的每个碱基所对应的测序ASCII质量值。

现在有一段很短的接头序列ATTC,如果匹配上seq.file中的序列,即匹配到每4行所代表的序列中的第2行碱基序列,则删除seq.file中的存在的接头序列ATTC和这部分接头序列对应的质量值,使得第2行的碱基序列和第4行的质量序列最终还是一一对应

这个目标有点像公共软件fastx_clipper的功能,只是这里想用perl来实现,该如何才能实现呢?希望得到大家的指点!!
作者: 104359176    时间: 2013-09-25 09:29
把你的问题抽象成一个 Perl 语言的算法问题。而不是讲解生物学。
作者: iamline    时间: 2013-09-25 09:42
104359176 发表于 2013-09-25 09:29
把你的问题抽象成一个 Perl 语言的算法问题。而不是讲解生物学。


其实我是perl新手,也不知道要用什么算法来做,整体的要求是每4行为一个单元,删除每个单元里第二行的接头序列及其对应的质量值序列(在第4行),其余不变,求大侠指点
作者: 104359176    时间: 2013-09-25 10:30
回复 3# iamline


    不好意思,我是生物学的门外汉,你说的那些术语,我根本不懂。
作者: xiumu2280    时间: 2013-09-25 10:40
本帖最后由 xiumu2280 于 2013-09-25 10:45 编辑

其实楼主就是想要个去除测序接头的程序···
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;

  4. while (<DATA>) {
  5.         chomp;
  6.         my $line = $_;
  7.         my $line1 = <DATA>;
  8.         my $line2 = <DATA>;
  9.         my $line3 = <DATA>;
  10.         if ($line1=~/^ATCA/) {
  11.                 $line1=substr ($line1,4);
  12.                 $line3=substr ($line3,4);
  13.         }
  14.         print "$line\n$line1$line2$line3";
  15. }
  16. __DATA__
  17. @HWI-ST531R:212:C22J6ACXX:8:2314:4159:37409
  18. AAAAAAAAAATACTGTTTATGATGATATCATTGAAGCTGCGGGTGAATTAGGTGATAAAACGCTACTTCCTGTTTTAGATACTATGTTGTACAAGTTTGAT
  19. +
  20. C@CFFFFFHHHHHIIHGIIJGIJIBGEIIGJCGHIJIJGIJJJ8>EDDFCED6@ACCCDDD@@BDDDDDDCDDCDCCCCADDDDECDDDCDDDCCDCDCAC
  21. @HWI-ST531R:212:C22J6ACXX:8:2102:3421:60710
  22. AAAAAAAAACCATCCAAATCTGGATGGCTTTTCATAATTCTGAGAAATTAGCTGCGCTGGCGCACCGCTTCAAATAAGCAAATTCCGGTCGCAACCGAAAC
  23. +
  24. @@CFFFFFHHHHHJJJJJJJJJJJJJJJJJJJJIJIIJJIJJIJIJJJJIIIJJHHFFDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD
  25. @HWI-ST531R:212:C22J6ACXX:8:1303:19050:54218
  26. AAAAAAAAACCGCACCAGAATTCCCATCGATCCAAGGCCGTTGGTGTTCCTTTCCCCACAAACCGCCCCAACGTCTCAATCGAGACTCCCCCTCACAGCAG
  27. +
  28. ??<BDDDDFFFF08@FFBG>D8?FFFF9;;==5@@=(66:55?B(5;?BA;AA>AA<<;<?B9?5@08-<@BB<2<9(4+:28+25@:>8>B@@B<A>?##
  29. @HWI-ST531R:212:C22J6ACXX:8:2108:4215:95315 1:N:0:GAGTGG
  30. AAAAAAAAACCGGGCAATTGCCCGGTGCTGGAAGGGATCTCGTTATTCGGGGGAATAAGGTAGATGTGGATAGTCCAGCCAGTGCGTCAAATAGTGAGAGA
  31. +
  32. <<@DDDDDAH@DHG:FHFGIIGCGI@FHHIHIIFHE6?1;?@?==ACDCA@BB5<8>C@C33:@@A@>>:>:?::>4>9<<A::>9(55>A:>A4(+4@32
  33. @HWI-ST531R:212:C22J6ACXX:8:2208:5557:5216 1:N:0:GAGTGG
  34. AAAAAAAAACCTGAAAAAAACGGCCTGACGTGAATCAAGCAACTTTTTTCAGGTTTTGCCCGCTTAGTGCGGTAACAATCCTTTACTCAGTAATAATATTT
  35. +
  36. @@CFFFFFFDFFFEHGIGIHGHHGIIGEDEHFEHFBC@?DBA(;@CCDDC@CDCCCDDCDDD;DB7C>C@BB05?A>@CDDDDCCDDDDCCCCDDDDECDD
  37. @HWI-ST531R:212:C22J6ACXX:8:2214:1365:3705
  38. AAAAAAAAACCTGAAAAAAACGGCCTGACGTGAATCAAGCAATTTTTTTCAGGTTTTGCCCGCTTAGTGCGGTAACAATCCTTTACTCAGTAATAATATTT
  39. +
  40. @@@FDDFFF:DHFHGIGHEHEFGEAFHGIIHGCEDFFDE@A@>;ACDDD:>@C(>CD?@@CD@B@BDCA>;;0;B>CC@C<CCDCDAA>CACCCDDEC>CB
  41. @HWI-ST531R:212:C22J6ACXX:8:1107:2469:40337
  42. ATCAAAAAAACCTGAAAAAAACGGCCTGACGTGAATCAAGCAATTTTTTTCAGGTTTTGCCCGCTTAGTGCGGTAACAATCCTTTACTCAGTAATAATATTT
  43. +
  44. CCCFFFFFHHHHHJJJJJJIJIJJJJJJJJHGHGFFEFEEECEEEEDDDCDDDACDDDDDDDDDDDDDDDDDBDDDDDDDDDDDDDDCDDDEDEEDEEEED
  45. @HWI-ST531R:212:C22J6ACXX:8:2211:2015:43926
  46. AAAAAAAAACCTGAAAAAAACGGCCTGACGTGAATCAAGCAATTTTTTTCAGGTTTTGCCCGCTTAGTGCGGTAACAATCCTTTACTCAGTAATAATATTT
  47. +
  48. CCCFFFFFFHHAFGGHDHIIIIIGHGGHGGHHAHEDFFFCAECCCEDDD>:CCCCDDDDDDD7@BDDCCDBB;B@CCDDCCDDDDDDCCCADDDEDEDEEE
复制代码

作者: iamline    时间: 2013-09-25 14:03
xiumu2280 发表于 2013-09-25 10:40
其实楼主就是想要个去除测序接头的程序···


非常感谢,帮助很大!!
作者: pitonas    时间: 2013-09-25 16:01
5 楼
很牛啊




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