私は、次世代シーケンシング イルミナのデータを分析するために、さまざまなスクリプトと入力パラメーターを一緒にパイプすることを含むバイオインフォマティクス プロジェクトに取り組んでいます。最初のスクリプトのデバッグに助けが必要です。そのタスクは、qseq ファイルを解析し、「適切な」サンプルを fastq 形式に変換し、出力を一時的な txt ファイル (ディスク) に保存することです。
デバッグの目的で、パイプ スキームは次のようになります。
# the input parameters are "fed" into the script and the output is written
# to tmp.txt
script01.pl [input parameters] > tmp.txt
このコマンドをターミナルに入力し、tmp.txt ファイルを調べて、スクリプトが期待どおりの結果を出力するかどうかを確認します。しかし、プロジェクト全体では、すべてのスクリプトをまとめてパイプ処理するラッパー スクリプトと呼ばれるものがあります。script01 は、スクリプト 02-06 で個別に処理する必要がある end1 および end2 データの出力を tmp ファイルに保存します。
これがコードです。何が起こっているのかを理解できるように、説明的なコメントを追加しました。また、表示する qseq ファイルはありませんが、フィールドを正しく解析していると仮定してください。
#!/usr/bin/perl
use strict; use warnings;
my $end1_temp=shift; # this variable is the location of the temporary file
my $end2_temp=shift; # this variable is the location of the temporary file
my $qseq_file=shift;
my $barcode=shift;
# declare and initialize variables
my $overhang= 'C[AT]GC';
my $end1_trim_offset= length($barcode)+4;
my $end2_trim_offset= 4;
my %to_keep=(); # an empty hash
my @line=(); # an empty array
# open the qseq file, or exit
open (QSEQ_1_FILE, $qseq_file) or die "couldn't open $qseq_file for read\n";
# also, open (for output) the end1 temp file, so we can write to it while
# we process the end1 input file above
(open END_1_FILEOUT,">$end1_temp") or die "couldn't open $end1_temp for write\n";
# reads each line of the qseq data file one at a time.
# assume each sample is kept on a separate line
while(<QSEQ_1_FILE>){
chomp; @line=split; # index and formats the qseq fields for parsing
# skip samples those that didn't pass QC
$line[10]>0 or next;
# process the samples here. look at only the samples that pass the quality
# control and whose sequences contain the barcode+overhang at the beginning
# of the string, otherwise skip to the next sample in the data file
# (i.e start the next iteration of the loop)
# trim the barcode+overhang from seq and qual
$line[8]= substr($line[8], $end1_trim_offset); #sequence
$line[9]= substr($line[9], $end1_trim_offset); #quality
# unique sequence identifier. don't worry about what $line[1-6] represent
# just know that $identifier is unique for each of the 'good' samples
my $identifier = $line[0].'-'.$line[1].'-'.$line[2].'_'.$line[3].'-'.$line[4].'-'.$line[5].'_#'.$line[6];
# store the identifiers of the 'good' samples in a hash.
# the hash should contain the identifier as the key and numbers (1,2,3,etc.)
# as the values. the following increments the hash values for each identifier.
$to_keep{$identifier}++;
# the following below is suppose to write information to the filehandle in fastq format:
# @[the identifier]/1
# sequence [ATGCAGTAAT...]
# +[the identifier]/1
# quality [ASCII characters]
print END_1_FILEOUT '@' . "$identifier/1\n$line[8]\n" . '+' . identifier/1\n$line[9]\n";
}
print "Found " . int(keys %to_keep) . " reads from end1\n";
# close the filehandles
close QSEQ_FILE; close END_1_FILEOUT;
この時点で、「適切な」シーケンスのみの識別子を含むハッシュがあり、fastq データが保存場所に書き込まれています。Script01 は、fastq 出力をディスク上の一時ファイルに保存することを想定しています。
$end1_temp = '~/tmp/sampleD1_end1.fq'とend2_temp = '~/tmp/sampleD1_end2.fq'
質問:print END_1_FILE上記の行は fastq データをファイルハンドルに書き込んでいますか、それとも $end1_temp 変数に書き込んでいますか? $end1_temp 変数と $end2_temp 変数を script02 に渡す必要があるためです。また、デバッグの目的で、script01 からの fastq 出力を表示するにはどうすればよいですか?
ヘルプが必要な残りのコードは次のとおりです。これは同じスクリプト上にあり、上記のコードに直接従います。
# if the sequence is paired (has both end1 and end2 data), then the qseq_2 file exists and the conditions evaluates to true
if ($end2_temp) {
# changes the qseq filename from file name from end1 to end2 data
# don't worry about why it works
$qseq_file=~ s/'_1_'/'_2_'/;
open (QSEQ_2_FILE, $qseq_file) or die "could not open $qseq_file for read\n";
# open (for output) the end2 temp file, so we can write to it while we process
# the end2 input file above
open (END_2_FILEOUT, ">$end2_temp") or die "could not open $end2_temp for write\n";
# reads each line of the end2 file one at a time
while(<QSEQ_2_FILE>){
# skip comments
/^\#/ and next;
chomp; #keep the qseq fields on one line.
my @line= split; #indexes the qseq fields.
# unique sequence identifier that preserves the sequencing information
# in other words, samples from end1 and end2 will have the same unique
# identifier because they contain the exact same fields in columns 0-6
# or $line[0-6]. also end1 and end2 have the same number of samples
my $identifier = $line[0] .'-'.$line[1].'-'.$line[2].'_'.$line[3].'-'.$line[4].'-'.$line[5].'_#'.$line[6];
# recall the %to_keep hash which stored the identifiers of the 'good' samples;
# it was inside the QSEQ_1_FILE loop
# here I only want the end2 samples whose identifiers match the identifiers
# from the end1 samples
# skip sample where end1 didn't pass; does this work??
# the condition is suppose to evaluate to false if the identifiers don't match
$to_keep{$identifier} or next;
# trim the barcode+overhang from seq and qual
$line[8]=substr($line[8], $end2_trim_offset); # sequence
$line[9]=substr($line[9], $end2_trim_offset); # quality
print END_2_FILEOUT '@' . "$identifier/2\n$line[8]\n" . '+' . "$identifier/2\n$line[9]\n";
}
close QSEQ_2_FILE; close END_2_FILEOUT;
}
これで script01 は終了です。この時点で、end1 と end2 の「適切な」サンプルの fastq データを別の保存場所に書き込む必要がありました。Script01 は、fastq 出力をディスク上の 2 つの一時ファイルに保存することを想定しています。私の質問はデバッグ目的だと思いますが、script01 によって作成された tmp ファイルを表示するにはどうすればよいですか?
最後に、コマンドscript01.pl [input parameters] > tmp.txtを Linux ターミナルに入力すると、script01 の出力が tmp.txt に保存されます。"Found X read from end1" は、end1 の読み取りを処理した後にスクリプトが出力するものです。ここで、X は %to_keep ハッシュ内の読み取りの数です。
tmp.txt を見ると、Found 0 reads from end1.0 が表示されるので、ハッシュに何も保存されていないことを意味します。end1 から約 630 万回の読み取りを保存すると想定されています。読み取りがハッシュに保存されない理由を誰かが理解するのを手伝ってくれますか?
問題は、ハッシュに保存する必要があるかどうかを決定するために使用している基準を読み取りが通過していないことだと思います。別の問題は、識別子の保存方法にある可能性があります。
皆さんはそれを見て、私が見逃したかもしれない何かがないか確認できますか?
ありがとう。私の質問に対する提案や回答は大歓迎です。