求一个切割序列的脚本

skyxbxb

序列：几万条，fasta格式，如下：
>gene=XLOC_000001
CTCCGTCCTGGATGCCGTCCTCCCTGGTGTCAGCCTCTTGGTCAAAAATTTTGGGGTCCTTTGTCTCCTTCGGTGGTGTCCATTGGTGCGGGTGCATGGATTCTGTGCGCGGGTGCATGGTGGCGGTCAAAACAGGTAGTATGGTGGCGGTCTCTGGATATGCGGGCTGCTCTTCCTGCGGTGATGGATCGTCGGCGGGGAATATAAGGCTCAAGAGGGGTTGATGTTTGCGGCTGGCTTTAATTGCATGTTTGGCGGCGGTTTTTTGTGGTAGGCTGTCTTGTATTGGCTTCGGTTCTTCTAATACTATTTTCGGTAAATCATTTCTAAAGGTGGTGTATGATGTTGAAAtatcttgattttctttgtttgtttcaTTGGAGCTTGTCGGTGTtatatcgggggtatcatatagtcctggtttattagttgtagttttctgtatatctggtggcgggtcgttgggtgctgggttccgggttttttgttgattcaatcc
>gene=XLOC_000002
TAGGGATTTGCTTGTGATCATGTTTCTATAAAGTCGtcattttatcaattcgaaGATCAGCCTTCAAAAATTTAGGTCAACTTATACATTTCAACTCACGTTGCTACACGTTTCAACTTAGCAAAAGGACAGAGggaattagtattattattctctagagcCCTGGAAATATCCCTGTTATTGTTCTTGGTACGGTACCTTATACAAAAGGGACtccataatacaaaaacaaatcagtCCTTGCCGGACAAagctagagaaagaagaagaagaagaagacgacgacaggaggagaagaagaacaagcgagaagaagaagataatagaaaaagg

求一个脚本可以把每条序列按顺序从头到尾切割成一些2000bp的短序列。

skyxbxb

比如说一条长序列 CTCCGTCCTGGATGCCGTCCTCCCTGGTGTCAGCCTCTTGGTCAAAAATTTTGGGGTC
我希望切割成
CTCCG
TCCTG
GATGC
CGTCC
TCCCT
GGTGT
CAGCC
.......
固定长度的短序列，例子中长度为5，我希望切割成长度为2000以内如1900的长度。
各位大神有何妙招

Okelani

perl cut.pl  your_data_file

1. #!/usr/bin/perl
   
2. use warnings;
   
3. use 5.010;
   
4. my $Bp = 5;    # 切割成短序列
   
5. 
   
6. while (<>) {
   
7.     print; $_ = <>; chomp;
   
8.     my $pos = 0;
   
9.     while ( $pos < length ) {
   
10.         say substr $_, $pos, $Bp;
    
11.         $pos += $Bp;
    
12.     }
    
13. }
    
14. 
    

复制代码

skyxbxb

回复 3# Okelani

多谢大神！不好意思我之前没说清楚，能不能切割序列并且给它加上编号呢，比如：
原文件：
>gene=XLOC_000001
CTCCGTCCTGGATGCCGTCCTCCCTGGTGTCAGCCTCTTGGTCAAAAATTTTGGGGTC
结果文件：
>gene=XLOC_000001.1
CTCCG
>gene=XLOC_000001.2
TCCTG
>gene=XLOC_000001.3
GATGC
>gene=XLOC_000001.4
CGTCC
.......
长序列切割成固定长度的短序列，并且编号（编号形式任意，包含原来的编号就好，比如>gene=XLOC_000001_1也行）。

多谢大神！

Okelani

回复 4# skyxbxb

perl cut.pl your_data
perl cut.pl your_data 100

1. #!/usr/bin/perl
   
2. use warnings;
   
3. use 5.010;
   
4. my $Bp = $ARGV[1] ? pop : 80;    # 切割成短序列, default = 80
   
5. 
   
6. while (<>) {
   
7.     # print;    # 原文件  >gene=XLOC_000001
   
8.     chomp;
   
9.     $_ .= '.';
   
10.     chomp( my $seq = <> );
    
11.     my ( $pos, $i ) = ( 0, 1 ); # $i = 编号
    
12.     while ( $pos < length($seq) ) {
    
13.         say $_, $i++;
    
14.         say substr $seq, $pos, $Bp;
    
15.         $pos += $Bp;
    
16.     }
    
17. }
    

复制代码

skyxbxb

回复 5# Okelani

非常感谢大神！！！

凋谢1991

回复 3# Okelani
好心的大神能不能帮我一下，帖子如下：
文件中定位查找特殊字段并输出到表格中
http://bbs.chinaunix.net/forum.p ... mp;fromuid=31517907
(出处: http://bbs.chinaunix.net/)

Okelani

回复 7# 凋谢1991

1. input data
2. procedure
3. output data

b114213903

1. #!/usr/bin/perl
   
2. use strict;
   
3. use Bio::SeqIO;
   
4. 
   
5. my ($file,$length)=@ARGV[0,1];  # perl Cut_fasta.pl  fasta_file  50
   
6. 
   
7. my $SeqIO=Bio::SeqIO->new(-file=>"$file",-format=>'fasta');
   
8. while (my $Seq=$SeqIO->next_seq){
   
9.         my $id=$Seq->id;
   
10.         my $len=$Seq->length;
    
11.         for (my $i=1;$i<=($len/$length+1);$i++){
    
12.                 my $start=$length*($i-1)+1;
    
13.                 my $end=($len<$length*$i)?$len:$length*$i;
    
14.                 my $seq=$Seq->subseq($start,$end);
    
15.                 print ">$id.$i\n$seq\n"
    
16.         }
    
17. }
    
18. $SeqIO->close();

复制代码