仿写fastqc软件的部分功能（上）

前面我们介绍了fastqc这个软件的使用方法 http://www.bio-info-trainee.com/?p=95 ，这是一个java软件，但是有些人服务器没有配置好这个java环境，导致无法使用，这里我贴出几个perl代码，也能实现fastqc的部分功能

统一测试文件是illumina的phred33格式的fastq文件，共100000/4=25000条reads，读长都是101个碱基

程序名-fastq2quality.pl

使用命令：perl fastq2quality.pl SRR504517_1.fastq >quality.txt

功能： 把fastq格式的每条原始reads的第四行ascii码质量值，转换为Q值并输出一个矩阵，有多少条reads就有多少行，每条reads的碱基数就是列数。

[perl]
while (<>){

    next unless $.%4==0;

    chomp;

    s/\r//g;

    @F=split//;

    foreach (@F){

        $num=ord($_);

        $num-=33;

        print "$num\t";

    }

    print "\n";

}
[/perl]

 

统计结果如下

程序名-fastq2meanQ.pl

使用命令：perl fastq2meanQ.pl SRR504517_1.fastq

功能： 把fastq格式的原始reads统计每条reads的平均Q值，并画出Q值1到50各有多少条reads的分布图

[perl]
while (<>){

    next unless $.%4==0;

    chomp;

    s/\r//g;

    @F=split//;

    $mean=0;

    $sum=0;

    foreach (@F){

        $num=ord($_);

        $num-=33;

        $sum+=$num;

    }

    $mean=int($sum/@F);

    $hash{$mean}++;

}

print "$_ \t$hash{$_}\n" foreach sort {$a<=>$b}  keys %hash;
[/perl]

 

统计结果如下

 

程序名-fastq2fivenum.pl

使用命令：perl fastq2fivenum.pl  SRR504517_1.fastq

功能： 把fastq格式的每条原始reads的第四行ascii码质量值，转换为Q值，并对每一个位点统计所以reads的四分位数，加上平均数。

[perl]
use List::Util qw/max min sum maxstr minstr shuffle/;

while (<>){

    next unless $.%4==0;

    chomp;

    s/\r//g;

    @F=split//;

    foreach (0..@F-1){

        $num=ord($F[$_]);

        $num-=33;

        $tmp[$_]->{$num}++;

    }

}

print "num\tmean\tmin\tq25\tq50\tq75\tmax\n";

$i=0;

foreach $hash (@tmp){

    $sum_reads=sum values %{$hash};

    $num_q25=int($sum_reads/4);

    $num_q50=int($sum_reads/2);

    $num_q75=int(3*$sum_reads/4);

    $sum_Q=0;

    $sum_value=0;

    foreach (sort {$a<=>$b} keys %{$hash}){

          #print "$_ \t$hash->{$_}------"   

            $sum_Q+=$_ * $hash->{$_};

            $q25_before=($sum_value<$num_q25);

            $q50_before=($sum_value<$num_q50);

            $q75_before=($sum_value<$num_q75);

            $sum_value+=$hash->{$_};

            $q25_last=($sum_value>$num_q25);

            $q50_last=($sum_value>$num_q50);

            $q75_last=($sum_value>$num_q75);

            $q25=$_ if $q25_before && $q25_last;

            $q50=$_ if $q50_before && $q50_last;

            $q75=$_ if $q75_before && $q75_last;            

        }

    $mean=$sum_Q/$sum_reads;

    $min=min keys %{$hash};

    $max=max keys %{$hash};   

    $i++;

    print "$i\t$mean\t$min\t$q25\t$q50\t$q75\t$max\n";

}
[/perl]

 

统计结果文件如下

最后一个，统计GC含量

这个我将会在下一篇讲诉如何用R画图