搜索
查看: 3177|回复: 2

fastq数据perl脚本过滤

[复制链接]

13

主题

66

帖子

276

积分

中级会员

Rank: 3Rank: 3

积分
276
QQ
发表于 2017-1-15 22:20:27 | 显示全部楼层 |阅读模式
     这里是使用perl脚本实现对fastq文件的过滤功能,具体思路很简单,针对PE测序,同时打开双端reads的fastq文件,然后按顺序对应依次读入4行,进行双端过滤!因为考虑到PE reads在是使用bwa map时针对单端read会报错,所以在过滤过程的中,一端read质量不合格,那么不管另一段质量如何,都应过滤掉!
    脚本如下,中间涉及到我根据测序文库构建试剂盒的分析说明而用到的过滤参数:
    [mw_shl_code=perl,true]
1 #!/home/huizhen/bin/perl/bin
  2
  3 use strict;
  4 use warnings;
  5 use Getopt:ong;
  6
  7 my $q ||= 20;
  8 my $p ||= 80;
  9 my $len ||= 45;
10 my $N ||= 10;
11 my($fq1,$fq2,$out1,$out2,$help);
12
13 GetOptions(
14     "q=i" => \$q,
15     "p=i" => \$p,
16     "len=i" => \$len,
17     "N=i" => \$N,
18     "fq1=s" => \$fq1,
19     "fq2=s" => \$fq2,
20     "out1=s" => \$out1,
21     "out2=s" => \$out2,
22     "help|?" => \$help,
23     );
24 my $usage=<<END;
25 Usages:
26     perl $0 [options]
27 Options:
28     -q <the min quality to keep>
29     -p <the min percent base to achive the min quality>
30     -len <the min length of base to keep>
31     -N <max N base to keep>
32     -fq1 <fastq_r1>
33     -fq2 <fastq_r2>
34     -out1 <filtered_fastq_r1>
35     -out2 <filtered_fastq_r2>
36     -help <this info>
37 END
38 die $usage if (! $fq1||! $fq2||! $out1||! $out2||$help);
39
40 open FQ1,"<",$fq1 or die "cant open!";
41 open FQ2,"<",$fq2 or die "cant open!";
42 open OUT1,">",$out1 or die "cant open!";
43 open OUT2,">",$out2 or die "cant open!";
44
45 my($total_before,$total_after) = (0,0);
46 while(1){
47     my $name_1 = <FQ1>;
48     my $seq_1 = <FQ1>;
49     my $plus_1 = <FQ1>;
50     my $qual_1 = <FQ1>;
51     my $name_2 = <FQ2>;
52     my $seq_2 = <FQ2>;
53     my $plus_2 = <FQ2>;
54     my $qual_2 = <FQ2>;
55     last unless $name_1;
56     $total_before ++;
57     #print"$name_1$seq_1$plus_1$qual_1";
58     chomp $qual_1;
59     chomp $qual_2;
60     my $len_1 = length($qual_1);
61     my $len_2 = length($qual_2);
62     next if ($len_1 < $len || $len_2 < $len);
63     my $n_1 = $seq_1 =~ tr/N/N/;
64     next if ($n_1/$len_1 >= $N/100);
65     my $n_2 = $seq_2 =~ tr/N/N/;
66     next if ($n_2/$len_2 >= $N/100);
67     my @f = split //,$qual_1;
68     my @r = split //,$qual_2;
69     my ($quality_1,$quality_2) = (0,0);
70     my ($count_1,$count_2) = (0,0);
71     my $sign = 0;
72     foreach(@f){
73         chomp;
74         $quality_1 = ord($_);
75         $quality_1 -= 33;
76         $count_1 ++ if $quality_1 >= $q;
77     }
78     $sign ++ if($count_1/$len_1 >= $p/100);
79     foreach(@r){
80         chomp;
81         $quality_2 = ord($_);
82         $quality_2 -= 33;
83         $count_2 ++ if $quality_2 >= $q;
84     }
85     $sign ++ if($count_2/$len_2 >= $p/100);
86     if($sign == 2){
87         print OUT1 "$name_1$seq_1$plus_1$qual_1\n";
88         print OUT2 "$name_2$seq_2$plus_2$qual_2\n";
89         $total_after ++;
90     }
91 }
92 close FQ1;
93 close FQ2;
94 close OUT1;
95 close OUT2;
96 my $percent =sprintf"%0.2f", $total_after/$total_before;
97 print"The reads inputtotal_before\nThe reads passedtotal_after($percent)\n";[/mw_shl_code]



上一篇:实现模糊匹配-高级难度
下一篇:【R处理芯片数据】--总结
苛求远离完美
回复

使用道具 举报

13

主题

66

帖子

276

积分

中级会员

Rank: 3Rank: 3

积分
276
QQ
 楼主| 发表于 2017-1-15 22:21:14 | 显示全部楼层
其实还可以添加更多的过滤分析,都是这个套路!
苛求远离完美
回复 支持 1 反对 0

使用道具 举报

633

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2017-1-16 11:02:52 | 显示全部楼层
赞,你这个流程,可以当做一个生信人必练数据处理题目了!
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2020-2-18 01:07 , Processed in 0.027024 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.