删除匹配的字符序列和其对应的质量序列

iamline

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

把你的问题抽象成一个 Perl 语言的算法问题。而不是讲解生物学。

iamline

104359176 发表于 2013-09-25 09:29
把你的问题抽象成一个 Perl 语言的算法问题。而不是讲解生物学。

其实我是perl新手，也不知道要用什么算法来做，整体的要求是每4行为一个单元，删除每个单元里第二行的接头序列及其对应的质量值序列（在第4行），其余不变，求大侠指点

104359176

回复 3# iamline


    不好意思，我是生物学的门外汉，你说的那些术语，我根本不懂。

xiumu2280

其实楼主就是想要个去除测序接头的程序···

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

复制代码

iamline

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

非常感谢，帮助很大！！

pitonas

5 楼
很牛啊