有时候我们需要得到序列的反向互补序列,可以用下面的脚本来批量处理序列。
用法:
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;
}
如果觉得我的文章对您有用,请随意打赏。你的支持将鼓励我继续创作!