3

我正在尝试调整 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 格式会更好吗？我只是在改编别人的剧本时才做到这一点！谢谢。

4

1 回答 1

5

你得到这个结果是因为 substr 只返回被替换的值，而不是它进行替换的整个字符串。很简单，您不需要将 substr 的返回值存储在 $seq3 中，因为（正如您所发现的）它只是复制了 $REF 中的内容：只需打印 $sequence 即可。

print OUT "$id\t$sequence\n";
于 2015-04-09T15:37:05.923 回答