使用 Bioperl 改变 fasta 文件中特定位置的核苷酸?
Using Bioperl to alter nucleotides at specific positions in fasta file?
我正在尝试修改 Bioperl 脚本以更改 fasta 文件中特定位置的核苷酸并输出具有更改序列的新文件。
fasta 输入示例:
>seq1
AAATAAA
更改文件的核苷酸位置示例:
##fileformat=VCFv4.1
##samtoolsVersion=0.1.18 (r982:295)
#CHROM POS REF ALT
seq_1 4 G A
我的脚本输出应该是:
seq_1 AAAGAAA
这是我当前的脚本:
#!/usr/bin/env perl
use strict;
use warnings;
use Bio::SeqIO;
use Bio::Tools::CodonTable;
use Bio::Seq;
my $original = shift @ARGV;
my $vcf = shift @ARGV;
my $outname = shift @ARGV;
# read in fasta file with gene sequences
my $in = Bio::SeqIO->new(-file => "$original" , '-format' => 'Fasta');
my $out = Bio::SeqIO->new('-format' => 'Fasta');
open (my $fh2, $vcf) or die "Error, cannot open file $vcf";
my @vcf= <$fh2>;
close ($fh2);
my $pos2;
while ( my $seq = $in->next_seq() ) {
my $id = $seq->id;
my $sequence = $seq->seq(); # get the sequence from the fasta file
# Search sequence in the vcf file and get the position of the SNP
foreach my $vcfline(@vcf){
if($vcfline =~ /$id/){
if($vcfline !~ /^#/){
$vcfline=~ s/\R//g;
my @vcfline= split(' ', $vcfline);
my $comp= $vcfline[0];
my $pos= $vcfline[1];
my $REF= $vcfline[2];
my $pos2=$pos-1; # correct position
# mutate the sequence
my $seq3=substr($sequence,$pos2,1,$REF);
open(OUT, ">> $outname");
print OUT
"$id\t$seq3\n";
close OUT;
}}}}
这目前只打印出一个带有序列 ID 和新核苷酸的文件(取自核苷酸变化文件的第 4 列),但我想要包含核苷酸变化的新序列。
很抱歉,我对 Perl 知之甚少,而且才刚刚开始使用 Bioperl,因此非常感谢您提供有关如何更改此脚本的指导。如果输出可以是 fasta 格式就更好了?我只是在改编别人的剧本时才设法做到这一点!谢谢
你得到这个结果是因为 substr returns 只是被替换的值,而不是它进行替换的整个字符串。很简单,您不需要将 substr 的 return 值存储在 $seq3 中,因为(如您所见)它只是复制 $REF 中的内容:只需打印 $sequence 即可。
print OUT "$id\t$sequence\n";
我正在尝试修改 Bioperl 脚本以更改 fasta 文件中特定位置的核苷酸并输出具有更改序列的新文件。
fasta 输入示例:
>seq1
AAATAAA
更改文件的核苷酸位置示例:
##fileformat=VCFv4.1
##samtoolsVersion=0.1.18 (r982:295)
#CHROM POS REF ALT
seq_1 4 G A
我的脚本输出应该是:
seq_1 AAAGAAA
这是我当前的脚本:
#!/usr/bin/env perl
use strict;
use warnings;
use Bio::SeqIO;
use Bio::Tools::CodonTable;
use Bio::Seq;
my $original = shift @ARGV;
my $vcf = shift @ARGV;
my $outname = shift @ARGV;
# read in fasta file with gene sequences
my $in = Bio::SeqIO->new(-file => "$original" , '-format' => 'Fasta');
my $out = Bio::SeqIO->new('-format' => 'Fasta');
open (my $fh2, $vcf) or die "Error, cannot open file $vcf";
my @vcf= <$fh2>;
close ($fh2);
my $pos2;
while ( my $seq = $in->next_seq() ) {
my $id = $seq->id;
my $sequence = $seq->seq(); # get the sequence from the fasta file
# Search sequence in the vcf file and get the position of the SNP
foreach my $vcfline(@vcf){
if($vcfline =~ /$id/){
if($vcfline !~ /^#/){
$vcfline=~ s/\R//g;
my @vcfline= split(' ', $vcfline);
my $comp= $vcfline[0];
my $pos= $vcfline[1];
my $REF= $vcfline[2];
my $pos2=$pos-1; # correct position
# mutate the sequence
my $seq3=substr($sequence,$pos2,1,$REF);
open(OUT, ">> $outname");
print OUT
"$id\t$seq3\n";
close OUT;
}}}}
这目前只打印出一个带有序列 ID 和新核苷酸的文件(取自核苷酸变化文件的第 4 列),但我想要包含核苷酸变化的新序列。
很抱歉,我对 Perl 知之甚少,而且才刚刚开始使用 Bioperl,因此非常感谢您提供有关如何更改此脚本的指导。如果输出可以是 fasta 格式就更好了?我只是在改编别人的剧本时才设法做到这一点!谢谢
你得到这个结果是因为 substr returns 只是被替换的值,而不是它进行替换的整个字符串。很简单,您不需要将 substr 的 return 值存储在 $seq3 中,因为(如您所见)它只是复制 $REF 中的内容:只需打印 $sequence 即可。
print OUT "$id\t$sequence\n";