使用Bioperl改变fasta文件中特定位置的核苷酸？

Question

我正在尝试调整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格式,那会更好吗？我只是设法做到这一点,因为我正在调整其他人的脚本!谢谢.

Answer 1

你得到这个结果是因为substr只返回被替换的值,而不是它替换它的整个字符串.很简单,你不需要在$ seq3中存储substr的返回值,因为(正如你所知)它只是复制$ REF中的内容:只需打印$ sequence.

print OUT "$id\t$sequence\n";