-1

私は 1 つの fasta ファイルと 1 つのテキスト ファイルを持っています。 fasta ファイルには fasta 形式のシーケンスが含まれ、テキスト ファイルには遺伝子の名前が含まれています。'>' 記号の後の fasta ファイルのシーケンスの名前をテキスト ファイルの遺伝子名に置き換えたいと思います。私はスクリプトを書いていますが、なぜそれが機能しないのかわかりません。誰かが私を助けてくれますか?

print"Enter annotated file...";
$f1=<STDIN>;
print"Enter sequence file...";
$f2=<STDIN>;
open(FILE1,$f1) || die"Can't open $f1";
@annotfile=<FILE1>;
open(FILE2,$f2) || die"Can't open $f2";
@seqfile=<FILE2>;
@d=split('\t',@annotfile[0]);

for($i=0;$i<scalar(@annotfile);$i++)
{
@curr_all=split('\t',@annotfile[$i]);
@curr_id[$i]=@curr_all[0];
@gene_nm[$i]=@curr_all[1];
}
for($j=0;$j<scalar(@seqfile);$j++)
{   
$id=@curr_id[$j];
$gene=@gene_nm[$j];


@seqfile[$j]=~s/$id[$j]/$gene[$j]/g;
print @seqfile[$j];
}

私のファイルは次のようになります。

annot.txt

pool75_contig_389 ユビキチン リガーゼ
e3a
_

goat300.fasta

goat300.fasta


>pool75_contig_704
CCCTTTCTCCCTTCCCAACATTCAGAGATACTGAATCGAAACTCTTACTGTCTGTTAGAT
GACAAAGAGTTATCCATCCTACATACTCCAATTTCCTTCCGCAACTTGTGATTTCGCCGC
TTGAATCTTGACGCCGTGCGTCCACAGTTTGTTGTGTTTTATCAATCAAGGTCATTATCA
ACCGAAGACGCTATCTATTTTCTTGGCGAAGCTCTCGGAAAGGAGCCATCGAAATGGAAG
TATTTCTCAAGAAAGTCCGCGAGTTATCCCGGAAGCAGTTC
>pool75_contig_389
GACCTATACCGGACCGTCACTGAAAGNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
ACGATCCAGGCATGGAGTTGTGGTGACGAGTAGGAGGGTCACCGTGGTGAGCGGGAAGCC
TCGGGCGTGAGCCTGGGTGGAGCCGCCACGGGTGCAGATCTTGGTGGTAGTAGCAAATAT
TCAAGTGAGAACCTTGAAGGCCGAGGTGGAGAAGGNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNTCATTTGTAT
CGCCCGGAAAACGTCACAAGAACGGGAGTTGCGTACAGAA
>pool75_contig_1977
AAGGGACACCGTTGGGTGAGGCGAGCTGCGTTCCTCGAACCATGGCTTCAAAAAGCGACT
TAGACCGTCAGATTGAACAGCTCAGGGCCTGCAAGCTCATTACAGAGGATGAGGTTAAGG
CACTCTGCGCTAAGGCGCGTGAGATTTTAATTGAAGAGAGTAATGTCCAGTGCGTGGACT
CACCTGTCACGGTTTGTGGCGATATCCACGGCCAGTTTTACGACTTGATTGAACTGTTTA
AAGTGGGCGGAGATGTTC
>pool75_contig_3064
TTACTATTTCTGGGCCTTAAGACTGGCTTAGTCGCTTACGACCCTTATAACAATGTAGAT
GTATATTATAAGGATCTTCCTGATGGTGCTAACGCTATGTTAATTTATTCAAACTCACCG
ACAAAGGAACAGAATATGCTTTGGCAGGTGGAAACTGTTCGATAATTGGATTGAACGACG
GCGGATGCGAGGTATTTTGGACAGTCACTGGCGACTCCGTTTGCTCTCTTTGCTCGATTA
AATCCGACAGCGATAAGTCAAGAGATTTTGTGGTTGGCTCTGAAGATTTTGACATCCGAA
TCTTCCATGGGGATGCCATAATATATGAAATCACGGAGTCTGATG
>pool75_contig_2499
AAGAGAAGAGGTGAGTTTGAGTATTGTTTGTGTGTGTGTGGTTGGGTGAGTGTGTGGTAT
GTGGTGTATGTGTGTGATGAATGTATGTGAAAGAGAGTGATGAATCTCATGGATATGTTC
GAGTTCGTGGTTTCCATTGATCGGTTATAGCCGAGATGATGGATGTGTTCCATGTGTCTG
ATTTCAGTTTAGGATTGTGTTGATGATGTTGATGATGAAAATTGTTGATGGTGATGACGA
TAGTGATGATGATGACGATGTTTCGGATAATGGTGATGATGATGATGGTTCCGACGATGA
TGTTTCGCTTGATGATGGTGATAATGATGACTCCGAAAATAACGTTGACTCGGATGAG
4

2 に答える 2

0

あなたのコードには多くの警告があり、アプローチは非効率的でした。最初に、実際に動作する Perl プログラムをお見せしましょう。後で説明します。

#!/usr/bin/perl
use strict;
use warnings;

# Read the annotations file
print"Enter annotated file...\n";
# my $f1 = <STDIN>;
my $f1 = 'annot.txt';
open(my $fh_annotations, '<', $f1) or die "Can't open $f1";
my @annotfile = <$fh_annotations>;
close $fh_annotations;

# Read the sequence file
print"Enter sequence file...\n";
# my $f2 = <STDIN>;
my $f2 = 'goat300.fasta';
open(my $fh_genes, '<', $f2) or die "Can't open $f2";
my @seqfile = <$fh_genes>;
close $fh_genes;

# Process the annotations data
my %names; # this hash is going to hold the names
foreach my $line (@annotfile) {
  chomp $line;                      # remove newline
  my @fields = split /\t/, $line;   # split into array
  $names{$fields[0]} = $fields[1];  # save in the hash as key->value pair
}

# Process the sequence data
foreach my $line (@seqfile) {
  # Look at each line
  if ($line =~ m/>(.+)$/) {
    # If there is a heading there, remember it...
    if (exists $names{$1}) {
      # ... check if we know a name for it and replace it in the line
      $line =~ s/($1)/$names{$1}/;
    }
  }
  # output the line (this would be done to another filehandle)
  print $line;
}

これにより、両方のファイルが読み取られ、メモリに保存されます。しかし、名前に対して 2 つの配列を構築しようとする代わりに、キーと値のペアである hashを使用しました。数値の代わりに名前があり、特定の並べ替えがない配列のようなものと考えてください。

これらの名前を設定したら、シーケンス ファイルを処理できます。各行を見て、>記号を探して見出しがあるかどうかを確認するだけです。そこにある場合 ($1かっこのために入ります)、ハッシュにハッシュ エントリ (with exists)があるかどうかを調べます%names。その場合、見出しを適切な名前に置き換えることができます。

その後、それを新しいファイルに書き出すことができます。印刷してるだけです。

他にもいくつかのテクニックを使用しました。残念ながら、人々が BioPerl のコンテキストで入手する文献はかなり時代遅れです。このアドバイスを参考にしてください。あなたの人生が楽になります。

  • 常に と を使用strictしwarningsます。彼らはあなたのコードの問題について教えてくれます。
  • 変数は必ず で宣言してくださいmy。これは、問題の先頭に変数を設定する必要がある他の言語とは異なります。必要な場所で宣言できます。変数は特定のスコープ内にのみ存在します。これは、最も近い囲み{と}括弧の間、またはブロックを意味します。
  • セキュリティのために、3 つの引数のオープンおよびレキシカル ファイル ハンドルを使用します。詳細はこちらをご覧ください。
  • Perl は、Cループforeachの代替として提供しています。forこの場合、それは物事をずっと簡単にしました。

このプログラムについてもう 1 つ: このサンプル データはかなり短いものでしたが、実際のデータはもっと大きいと思います。メモリ不足にならないように、読み取り中にシーケンス ファイルを処理することを検討してください。他の操作を行う場合を除き、すべての行を保存する必要はありません。

open my $fh_out, '>', $filename_out or die $!;
open my $fh_in, '<', $filename_in or die $!;
while (my $line = <$fh_in>) {
  # do stuff with the line, like your regex
  print $fh_out $line;
}
close $fh_in;
close $fh_out;
于 2013-03-23T10:38:40.457 に答える
0

自分で行うのではなく、 Bio::SeqIOを使用してFastaデータセットを解析することを検討してください。Bio::SeqIO はこのタスクのために生きており、十分に開発されています。さらに、バイオインフォマティクスに携わっているのであれば、Bio::SeqIO を理解するのに役立ちます。これを踏まえて、次のことを考慮してください。

use strict;
use warnings;
use Bio::SeqIO;

open my $fh, '<', 'annot.txt' or die $!;
my %annot = map { /(\S+)\s+(.+)/; $1 => $2 } <$fh>;
close $fh;

my $in = Bio::SeqIO->new( -file => 'goat300.fasta', -format => 'Fasta' );

while ( my $seq = $in->next_seq() ) {
    my $seqID = $annot{ $seq->id } // $seq->id;
    print "$seqID\n" . $seq->seq . "\n";
}

データセットの出力:

tumor susceptibility
CCCTTTCTCCCTTCCCAACATTCAGAGATACTGAATCGAAACTCTTACTGTCTGTTAGATGACAAAGAGTTATCCATCCTACATACTCCAATTTCCTTCCGCAACTTGTGATTTCGCCGCTTGAATCTTGACGCCGTGCGTCCACAGTTTGTTGTGTTTTATCAATCAAGGTCATTATCAACCGAAGACGCTATCTATTTTCTTGGCGAAGCTCTCGGAAAGGAGCCATCGAAATGGAAGTATTTCTCAAGAAAGTCCGCGAGTTATCCCGGAAGCAGTTC
ubiquitin ligase e3a
GACCTATACCGGACCGTCACTGAAAGNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNACGATCCAGGCATGGAGTTGTGGTGACGAGTAGGAGGGTCACCGTGGTGAGCGGGAAGCCTCGGGCGTGAGCCTGGGTGGAGCCGCCACGGGTGCAGATCTTGGTGGTAGTAGCAAATATTCAAGTGAGAACCTTGAAGGCCGAGGTGGAGAAGGNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNTCATTTGTATCGCCCGGAAAACGTCACAAGAACGGGAGTTGCGTACAGAA
serine threonine-protein phosphatase 4 catalytic subunit
AAGGGACACCGTTGGGTGAGGCGAGCTGCGTTCCTCGAACCATGGCTTCAAAAAGCGACTTAGACCGTCAGATTGAACAGCTCAGGGCCTGCAAGCTCATTACAGAGGATGAGGTTAAGGCACTCTGCGCTAAGGCGCGTGAGATTTTAATTGAAGAGAGTAATGTCCAGTGCGTGGACTCACCTGTCACGGTTTGTGGCGATATCCACGGCCAGTTTTACGACTTGATTGAACTGTTTAAAGTGGGCGGAGATGTTC
bardet-biedl syndrome 2 protein P
TTACTATTTCTGGGCCTTAAGACTGGCTTAGTCGCTTACGACCCTTATAACAATGTAGATGTATATTATAAGGATCTTCCTGATGGTGCTAACGCTATGTTAATTTATTCAAACTCACCGACAAAGGAACAGAATATGCTTTGGCAGGTGGAAACTGTTCGATAATTGGATTGAACGACGGCGGATGCGAGGTATTTTGGACAGTCACTGGCGACTCCGTTTGCTCTCTTTGCTCGATTAAATCCGACAGCGATAAGTCAAGAGATTTTGTGGTTGGCTCTGAAGATTTTGACATCCGAATCTTCCATGGGGATGCCATAATATATGAAATCACGGAGTCTGATG
succinyl- ligase
AAGAGAAGAGGTGAGTTTGAGTATTGTTTGTGTGTGTGTGGTTGGGTGAGTGTGTGGTATGTGGTGTATGTGTGTGATGAATGTATGTGAAAGAGAGTGATGAATCTCATGGATATGTTCGAGTTCGTGGTTTCCATTGATCGGTTATAGCCGAGATGATGGATGTGTTCCATGTGTCTGATTTCAGTTTAGGATTGTGTTGATGATGTTGATGATGAAAATTGTTGATGGTGATGACGATAGTGATGATGATGACGATGTTTCGGATAATGGTGATGATGATGATGGTTCCGACGATGATGTTTCGCTTGATGATGGTGATAATGATGACTCCGAAAATAACGTTGACTCGGATGAG

ハッシュ%annotは、データの内容を読み取ってキャプチャすることによって初期化されますannot.txt。Bio::SeqIO オブジェクトは、goat300.fastaファイル データを使用して作成されます。ループはwhilefasta シーケンスを反復します。変数$seqIDは、ハッシュ内のキーに関連付けられた値を取得するか%annot、現在のシーケンス ID を保持します (この//表記は、定義済みまたはを意味するため、$seqID が確実に定義されます)。最後に、Fasta レコードが印刷されます。

お役に立てれば！

于 2013-03-23T17:32:12.317 に答える