求序列重复的问题？

leigh111

大家好：
  我有一些字母组成的序列。我想让它们每60个连续字母为一组，而且上一组的最后25个连续字母 和 下一组的前面25个连续的字母应该是重复的。以此类推。 这程序现在没一点思路。是不是应该建一个子程序？
    望大家给予帮助。
例子是这样的：
qwertyyuuiiooasdfghjklzxcvbnmqweeaqwdrfgthyjuikmnbvcxzsaaqwertgfdsssyuimnbvcxzaswqerttyyuuuijmnbvcx

我要的结果是这样的：
qwertyyuuiiooasdfghjklzxcvbnmqweeaqwdrfgthyjuikmnbvcxzsaaq  #比如说这是60个，那么第二组是：
aqwdrfgthyjuikmnbvcxzsaaqwertgfdsssyuimnbvcxzaswqerttyyuuu #加粗的连续字母就是重复的。这样以此类推。

longbow0

是不是想把一个长的序列，如基因组，分成 60 bp 的片段，相邻 2 段的首尾 25 字符 一样?

1. 
   
2. 
   
3. my $seq = 'qwertyyuuiioo...';
   
4. my $offset = 0；
   
5. my $len = 60;
   
6. 
   
7. my $i = 0;
   
8. 
   
9. while ( $i < ( length($seq)/35) ) {
   
10.     my $subseq = substr($seq, $offset, $len);
    
11. 
    
12.     (do sth.)
    
13.    
    
14.     $offset +=35;
    
15.     $i++;
    
16. }
    

复制代码

leigh111

对的。您很厉害。我刚接触生物这块。呵呵。这里要用到关于很多perl 的东西。
这是我刚才写的。写到这里真是没思路了。呵呵

use strict;
use warnings;
open(FH,'D:\sequences.fasta') || die "$!";
my $n = 0;
my $length = 60;
my $i = 0;

while (<FH>) {
    my $value = substr ($_,$n,60);
    last if (length $value < 35);
    print "$value\n";
}

您的方法还是好。我仔细研究一下。谢谢您了

兰花仙子

回复 2# longbow0


    的确不错。不过这个length($seq)/35 + 1这个+1似乎是多余的。

1. use strict;
   
2. use Data::Dumper;
   
3. 
   
4. my $seq = "qwertyyuuiiooasdfghjklzxcvbnmqweeaqwdrfgthyjuikmnbvcxzsaaqwertgfdsssyuimnbvcxzaswqerttyyuuuijmnbvcx";
   
5. 
   
6. my $offset = 0;
   
7. my $len = 60;
   
8. 
   
9. my $i = 0;
   
10. my @result;
    
11. 
    
12. while ( $i < ( length($seq)/35) ) {
    
13.     my $subseq = substr($seq, $offset, $len);
    
14.     push @result,$subseq;
    
15. 
    
16.     $offset +=35;
    
17.     $i++;
    
18. }
    
19. 
    
20. print Dumper \@result;

复制代码

output:

$VAR1 = [
          'qwertyyuuiiooasdfghjklzxcvbnmqweeaqwdrfgthyjuikmnbvcxzsaaqwe',
          'wdrfgthyjuikmnbvcxzsaaqwertgfdsssyuimnbvcxzaswqerttyyuuuijmn',
          'imnbvcxzaswqerttyyuuuijmnbvcx'
        ];

longbow0

回复  longbow0


    的确不错。不过这个length($seq)/35 + 1这个+1似乎是多余的。




output:
...
兰花仙子 发表于 2010-03-13 17:08

   谢谢 仙子 提醒，确实是考虑不周。

leigh111

还得求大家一个忙，就是读入一个文件的话。
比如说文件里是：
fgheiojwdcasdncjnksacnsadahhawkdnjkaskjdhukwhuiadckjsncmnawkaw
twgdhjdvsbcmabwhdjbcmnsb hfgvjagwgdrfhjfmasvdghafywgdvsvcdhvhj
wadjwbjdbahvwdyawvbhjdvhwjagdyagdhjvqawyfdyqvvghsfdyasfyudgw
awdawdcadawdcdffehfuiehwhbefehfuewj

因为while是一行一行读取的。当open这个文件时，我想让它们组成一个字符串，就是第一行和第二行没有换行符的。依此类推。
像这样：
my $seq = "fgheiojwdcasdncjnksacnsadahhawkdnjkaskjdhukwhuiadckjsncmnawkawfgheiojwdcasdncjnksacnsadahhawkdnjkaskjdhukwhuiadckjsncmnawkawwadjwbjdbahvwdyawvbhjdvhwjagdyagdhjvqawyfdyqvvghsfdyasfyudgwawdawdcadawdcdffehfuiehwhbefehfuewj";
就是去掉每一行最后的换行符。该怎么实现呢？
我用 chop函数不行。用join 函数也没实现成功。

longbow0

处理序列文件的时候，像 GenBank、Fasta 等，最好是用 BioPerl 读写，可以避免在输入输出格式上的麻烦。

1. 
   
2. #
   
3. # 从 Fasta 文件中读取序列
   
4. #
   
5. use Bio::SeqIO;
   
6. 
   
7. my $o_seqi = Bio::SeqIO->new(
   
8.     -file => 'seq.fasta',
   
9.     -format => 'fasta',
   
10. );
    
11. 
    
12. while (my $o_seq = $o_seqi->next_seq) {
    
13.     my $seq = $o_seq->seq;    # 这里 $seq 就是你所需要的除去了换行、空白等等的序列字符串
    
14.     ...
    
15. }
    

复制代码

BioPerl - http://www.bioperl.org/

CPAN
Bio::SeqIO - http://search.cpan.org/~cjfields/BioPerl-1.6.1/Bio/SeqIO.pm
Bio::Seq - http://search.cpan.org/~cjfields/BioPerl-1.6.1/Bio/Seq.pm

yybmsrs

还得求大家一个忙，就是读入一个文件的话。
比如说文件里是：
fgheiojwdcasdncjnksacnsadahhawkdnjkaskjd ...
leigh111 发表于 2010-03-14 03:54

    用chomp()可以去除换行，然后join

leigh111

对的。搞个半天。其实是自己 忽略一个小点
open(FH,'D:\sequences.fasta') || die "$!";
my $string;
while (<FH>) {
    chomp;
  $string .= $_;   #这一步就搞定了。呵呵。我说怎么chomp不行呢。

谢谢大家的帮助。