perl脚本批量生成反向互补序列

有时候我们需要得到序列的反向互补序列,可以用下面的脚本来批量处理序列。

有时候我们需要得到序列的反向互补序列,可以用下面的脚本来批量处理序列。

用法:

perl  fasta_Reverse_complementary.pl  -fa input.fa  -out out.fa

input.fa是输入文件,out.fa是反向互补后的序列

输入文件格式入下所示:

>bta-26a-2
GGCUGUGGCUGGAUUCAAGUAAUCCAGGAUAGGCUGUUUCCAUCUGUGAGGCCUAUUCUU
GAUUACUUGUUUCUGGAGGCAGCU
>bta-18b
CUUGUGUUAAGGUGCAUCUAGUGCAGUUAGUGAAGCAGCUCAGAAUCUACUGCCCUAAAU
GCUCCUUCUGGCACA
>bta-29a
AUGACUGAUUUCUUUUGGUGUUCAGAGUCAAUAUAAUUUUCUAGCACCAUCUGAAAUCGG
UUAU
>bta-7f-2
UGUGGGAUGAGGUAGUAGAUUGUAUAGUUUUAGGGUCAUACCCCAUCUUGGAGAUAACUA
UACAGUCUACUGUCUUUCCCACG

代码如下:

#!/usr/bin/perl -w
use strict;
use Getopt::Long;
use Config::General;
use Bio::SeqIO;
use Bio::Seq;
my $version = "1.3";
## prepare parameters #######################################################################
## -------------------------------------------------------------------------------------------
## GetOptions
my %opts;
GetOptions(\%opts,  "fa=s", "out=s","h");
if(!defined($opts{out}) || !defined($opts{fa}) ||defined($opts{h}))
{
print <<"Usage End.";
Usage
Forced parameter:
-out          outfile                        <outfile>        must be given
-fa fasta file <infile> must be given
Other parameter:
-h           Help document
Usage End.
exit;
}
my $in  = Bio::SeqIO->new(-file => "$opts{fa}" , -format => 'Fasta');
open(OUT,">$opts{out}") ||die "open file $opts{out} faild.\n";

while ( my $seq = $in->next_seq() ) {
my($id,$sequence)=($seq->id,$seq->seq);
$sequence = &reverse_complement_IUPAC($sequence);
print OUT ">$id\n$sequence\n";
}
sub reverse_complement_IUPAC {
        my $dna = shift;
                # reverse the DNA sequence
                my $revcomp = reverse($dna);
                # complement the reversed DNA sequence
                $revcomp =~ tr/ABCDGHMNRSTUVWXYabcdghmnrstuvwxy/TVGHCDKNYSAABWXRtvghcdknysaabwxr/;
        return $revcomp;
}
sub reverse_complement {
        my $dna = shift;
                # reverse the DNA sequence
                my $revcomp = reverse($dna);
                # complement the reversed DNA sequence
                $revcomp =~ tr/ACGTacgt/TGCAtgca/;
        return $revcomp;
}



  • 发表于 2018-08-02 17:14
  • 阅读 ( 5510 )
  • 分类:perl

0 条评论

请先 登录 后评论
安生水
安生水

351 篇文章

作家榜 »

  1. omicsgene 702 文章
  2. 安生水 351 文章
  3. Daitoue 167 文章
  4. 生物女学霸 120 文章
  5. xun 82 文章
  6. rzx 78 文章
  7. 红橙子 78 文章
  8. CORNERSTONE 72 文章