(10)仿写fastqc-生信菜鸟团博客2周年精选文章集

用仿写软件的方法来学习编程

我首先仿写了fastqc软件,学会了很多基础知识:

仿写fastqc软件的一些功能-R代码

仿写fastqc软件的部分功能-perl代码

仿写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含量

程序名-fastq2meanGC.pl

使用命令:perl fastq2meanGC.pl SRR504517_1.fastq

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

[perl]</pre> while (<>){

next unless $.%4==2;

chomp;

s/\r//g;

@F=split//;

$GC=0;

foreach (@F){

$GC++ if /G/;

$GC++ if /C/;

}

#print "$GC\n";

$GC=int(100*$GC/length);

$hash{$GC}++;

}

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

结果如下所示

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

仿写fastqc软件的一些功能(下)

文件来自于上面perl代码的输出文件,好像算法有点问题,26G的文件居然处理近一个小时才出数据!

R语言本身自带的画图工具都很丑,懒得说了,可以用ggplot2来重新画一个,不是项目要求没有报酬我就懒得画了,大家面前看看画图原理即可。

a=read.table(“meanQ.txt”)

看看数据结构如下

> head(a)

V1 V2

1 2 93879

2 3 17800

3 4 25295

4 5 33259

5 6 55685

6 7 84866

plot(a,type=’l’,col=’red’,ylab=’reads number’,xlab=’mean quality’,main=’mean Q distribution’)

可以看出绝大部分的reads的Q值都在30-35直接,也就是说本次测序挺符合要求的,但是还是需要对那些平均Q20以下的reads过滤掉。

a=read.table(‘meanGC.txt’)

看看数据结构如下

> head(a)

V1 V2

1 0 503

2 1 151

3 2 163

4 3 179

5 4 315

6 5 443

plot(a,type=’l’,col=’red’,ylab=’reads number’,xlab=’reads bp’,main=’GC% distribution’)

可以看出GC含量的分布看起来挺符合正态分布的,大部分reads的GC含量都是在40%-60%直接

a=read.table(‘fivenum.txt’,header=T)

看看数据结构如下

boxplot(t(a[,3:7]),xlab=’reads bp’,ylab=’Q value’,main=’mean Q boxplot’)

可以看出测序质量从1-100bp过去质量越来越差,但是大部分都是高于Q30,但是88bp之后的碱基测序质量不咋地,可能需要trim掉

对于这个数据还可以画一个图

plot(a[,1:2],type=’l’,col=’red’,ylab=’Q value’,xlab=’reads bp’,main=’mean Q value distribution’)

可以看到88bp之后的平均Q值小于30,根据我们的阈值可能要把所有的reads的后面约10个bp的碱基要trim掉

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-01-06

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏玉树芝兰

如何用Python提取中文关键词?

本文一步步为你演示,如何用Python从中文文本中提取关键词。如果你需要对长文“观其大略”,不妨尝试一下。

972
来自专栏机器学习和数学

[Tensorflow] Tensorflow中模型保存与回收的简单总结

今天要聊得是怎么利用TensorFlow来保存我们的模型文件,以及模型文件的回收(读取)。刚开始接触TensorFlow的时候,没在意模型文件的使用,只要能...

3458
来自专栏大数据文摘

GitHub排名前20的Pandas, NumPy 和SciPy函数

2287
来自专栏申龙斌的程序人生

解密310 BTC(2)

价值1400万的比特币猜谜游戏刚火了几天,大奖就被一位高手全部取走,310 BTC的破解过程现在还没有公开,但已经有黑客公布了第二关的解法视频,过程相当复杂,我...

791
来自专栏应用案例

如何用Python提取中文关键词?

本文一步步为你演示,如何用Python从中文文本中提取关键词。如果你需要对长文“观其大略”,不妨尝试一下。 ? 需求 好友最近对自然语言处理感兴趣,因为他打算利...

2368
来自专栏ATYUN订阅号

人工智能为什么能做的事情这么多?密码猜测在它面前也能行得通

密码猜测之所以有效是因为… 人类是可预测的 ? 如果你要求小明设置一个密码。他可能只是简单地把密码设置为“xiaoming”。现在,系统告诉他密码必须包含数字...

3256
来自专栏华章科技

独家 | Python数据分析入门指南

有一个朋友最近问到这个问题,我觉得把它公开出来对其他人也会有帮助。这是给完全不了解Python而想找到从零到一的最简单的路径的人的建议:

753
来自专栏技术小黑屋

Gmail垃圾邮件过滤器文件分享

Gmail垃圾邮件过滤器文件。 A filter file for Gmail to auto-delete spams. 工作后,一直使用Gmail邮件托...

832
来自专栏生信技能树

生物信息学技能面试题(第4题)-多个同样的行列式文件合并起来

相信用过htseq-count的朋友都知道,它是分开对每个样本计算所有的基因表达量,所以会生成一个个独立的文件,我用perl脚本模仿它的结果如下: $ head...

3077
来自专栏机器之心

教程 | 用摄像头和Tensorflow.js在浏览器上实现目标检测

1844

扫码关注云+社区