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 に答える