perl get_fq_by_id.pl <id><fq1><fq2><OUT1><OUT2>
给一个fastqID列表,输出指定ID对应的双端fastq文件:
die "perl $0 <id><fq1><fq2><OUT1><OUT2>" unless(@ARGV==5);
use Bio::SeqIO;
use Bio::Seq;
open my $FQ1 ,"zcat $ARGV[1] |" or die "$!";
open my $FQ2 ,"zcat $ARGV[2]|" or die "$!";
my$fq1=Bio::SeqIO->new(-fh=>$FQ1,-format=>'fastq');
my$fq2=Bio::SeqIO->new(-fh=>$FQ2,-format=>'fastq');
open my $GZ1 ,"| gzip >$ARGV[3]" or die $!;
open my $GZ2 ,"| gzip >$ARGV[4]" or die $!;
my$fqo1=Bio::SeqIO->new(-fh=>$GZ1,-format=>'fastq');
my$fqo2=Bio::SeqIO->new(-fh=>$GZ2,-format=>'fastq');
my%keep=(
);
open IN ,"$ARGV[0]" or die "$!";
while(<IN>){
chomp;
next if /^#/;
my@tmp=split(/\s+/);
$keep{$tmp[0]}=1;
}
close(IN);
my$i=0;
while ( my $obj1=$fq1->next_seq() and my $obj2=$fq2->next_seq() ) {
my$id1=$obj1->id;
my$id2=$obj2->id;
#print "$id1,$id2\n";
#die;
if( exists $keep{$id1} or exists $keep{$id2}){
$fqo1->write_seq($obj1);
$fqo2->write_seq($obj2);
}
}
$fq1->close();
$fq2->close();
$GZ1->close();
$GZ2->close();
如果觉得我的文章对您有用,请随意打赏。你的支持将鼓励我继续创作!