一步一步教你学习Perl

Posted 2021-04-17 生信人

perl在处理字符串方面优势明显，尽管目前好多人推荐写Python，但从目前来说好多流程都是用perl来写的，你周围也有好多人用perl，所以今天小编讲一下，怎么写perl。

今天主题为处理FASTA序列相关脚本。

1、根据序列ID提取相应的序列

#! /usr/bin/perl 

use Bio::SeqIO;###载入模块，use XXX

my $genome=$ARGV[0];

my $seqID=$ARGV[1];

my $seq=$ARGV[2];

my %hash;

my @ID;

###ARGV数组为系统默认数目，用于从屏幕依次读入数据，分别存储到第一个元素，第二个元素，第三个元素；然后本处分别赋给了三个变量。

open (OUT1, ">$seq")|| die "cannot open $seq:$!"；

###打开输出文件写法。注意双引号中“>”，与下面的打开输入文件（打开文件为“<”）区别。

open (IN2, "<$seqID")|| die "cannot open $seqID:$!";

###打开输入文件写法。

my $ina = Bio::SeqIO->new(-file =>$genome, -format => 'fasta');

while(my $obj = $ina->next_seq()){

        my $id = $obj->id;

        my $seq = $obj->seq;

      $hash{$id}=$seq;

}

###上面的为Bio::SeqIO这一模块的用法；$id和$seq存储了序列的ID和序列；以ID为键，序列为值，存储到哈希中；所谓哈希可以简单认为是一个对应关系，如身份证号->小编。注意一个身份证号只能对应小编我一个人（换句话说你是不可能同时拥有两个身份证号的），这是哈希的基本准则。存储到哈希好处是我可以随便根据ID调出对应的序列。

while(<IN2>){###读取要提取的ID的文件了，<>表示按行读取

   next if(/^$/);###跳过空行。^行的首端；$行的末端

   my @tmp=split;###将行进行分割，分隔符号为空格或者tab，然后存储临时数组中

     push @ID, @tmp;###将数组@tmp进一步移入到@ID数组中

}

for($i=0;$i<=$#ID;$i++){###for循环了

print OUT1 ">$ID[$i]\n$hash{$ID[$i]}\n";###哈希用上；

}

close OUT1;###关闭文件

close IN1;###关闭文件

close IN2;###关闭文件

2、计算序列长度脚本

#! /usr/bin/perl 

use Bio::SeqIO;

my $seq=$ARGV[0];

my $Len=$ARGV[1];

open (OUT, ">$Len")|| die "cannot open $Len:$!";

my $ina = Bio::SeqIO->new(-file => $seq, -format => 'fasta');

while(my $obj = $ina->next_seq()){

        my $id = $obj->id;

        my $seq = $obj->seq;

       my seq_len=length($seq);

       print OUT ">$id\t$seq_len\n";

}

close OUT；

####本例与上例高度相似，只是使用了之前讲述的函数length；

3、提取指定位置序列脚本

#! /usr/bin/perl 

use Bio::SeqIO;

my $POS=$ARGV[0];

####输入要求:

###ID start stop strand

###seq1 100 300 +

my $SEQ=$ARGV[1];

my $Subseq=$ARGV[2];

my %hash;

my @position;

open (IN1, "<$POS")|| die "cannot open $POS:$!";

open (OUT, ">$Subseq")|| die "cannot open $Subseq:$!";

while (<IN1>){

my @tmp=split;

        push @position, [@tmp];

}

my $ina = Bio::SeqIO->new(-file => $SEQ, -format => 'fasta');

while(my $obj = $ina->next_seq()){

          my $id = $obj->id;

          my $seq = $obj->seq;

  $hash{$id}=$seq;

}

for($i=0;$i<=$#position;$i++){

  my  $ID=">$position[$i][0]_$position[$i][1]";

my $len=$position[$i][2]-$position[$i][1]+1;

  my $extract_seq=substr($hash{$position[$i][0]},$position[$i][1]-1,$len);

###substr字符串操作函数

if($position[$i][3] eq "-"){

$extract_seq=reverse_complementary($extract_seq);

}

print OUT "$ID\n$extract_seq\n";

}

close IN1;

close OUT;

####子函数写法

sub reverse_complementary{

my ($inputseq)=@_;

my $sequence = reverse(substr($inputseq,0));

$sequence =~ tr/ACGTUNacgtun/TGCAANtgcaan/d;

$inputseq=substr($sequence, 0);

return $inputseq;

}#####反向互补序列的输出！

你 GET到了吗？不懂请留言！