要求:
FASTQ文件长度过滤:编写一个程序,读取FASTQ文件,过滤掉碱基序列长度在60-80之外的序列,将长度在60-80之内的序列输出到结果文件中。
FASTQ格式文件如下:每四行表示一个测序序列,第二行是碱基序列。
@M2,HWI-7001455:326:HGTFVADXX:1:1101:2666:105721:N:0:GCGCTA
TTTAGTTTTGTAGTAATTGTTTGTAGTAAAATTTGTATTAGTTTTTTTATTTGTA
+
CCCFFFDDHHFHHIJJJJJHHIJJJJIIGEIJIJJIIJJJJHIJJJJJIJIJIHI
编程要求如下:
1) 程序命名为fq_filter.pl
2) 程序采用Getopt从命令行输入,参数共有3个,-infile,-outfile,-h;
3) 其中-infile用于接收输入的FASTQ文件名;
4) -outfile用于给出符合条件的结果文件名,输出结果格式要与输入的FASTQ格式相同。
5) -h用于给出程序的使用说明;
fq_filter.pl:
#!/usr/bin/perl -w
#use strict;
use Getopt::Long;
my ($chr,$help);
my $infile="";
my $outfile="";
GetOptions(
"help|h" => \&USAGE,
"infile:s" => \$infile,
"outfile:s" => \$outfile,
);
open IN,"<$infile"||die;
open OUT,">$outfile";
my $i=0;
my $seqID='';
my $sequence='';
my $len=0;
my $mark='';
my $qual='';
while(<IN>){
chomp;
if($i%4 == 0){
$seqID=$_;}
if($i%4 == 1){
$sequence=$_;
$len=length($sequence);}
if($i%4 == 2){
$mark=$_;}
if($i%4 == 3){
if($len>=60 && $len<=80){
$qual=$_;
print OUT "$seqID\n$sequence\n$mark\n$qual\n";}
}
$i=$i+1
}
sub USAGE{
my $usage=<<"USAGE";
------------------------------------------------------------------------------------
Program: thir_test.pl
Date:2017-05-19
Usage:
-infile <fasta format> infile
-outfile <fasta format> outfile
-h help
Example1: perl thir_test.pl -infile data.fa -outfile out.fa
Example2: perl thir_test.pl -h
------------------------------------------------------------------------------------
USAGE
print $usage;
exit;
}
运行命令:
perl fq_filter.pl -infile test.fq -outfile out.fa