给一个fastq ID列表,输出指定ID对应的双端fastq文件:

给一个fastqID列表,输出指定ID对应的双端fastq文件:

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();
  • 发表于 2020-05-25 11:32
  • 阅读 ( 3276 )
  • 分类:perl

0 条评论

请先 登录 后评论
omicsgene
omicsgene

生物信息

698 篇文章

作家榜 »

  1. omicsgene 698 文章
  2. 安生水 347 文章
  3. Daitoue 167 文章
  4. 生物女学霸 120 文章
  5. xun 82 文章
  6. 红橙子 78 文章
  7. rzx 74 文章
  8. CORNERSTONE 72 文章