短读质量控制 Read Filtration by Quality
由于各种因素,测序过程中一定会存在着错误,要么读错,要么漏读、要读多了。解决这种错误可以从源头,也就是测序仪改进,也可以通过生物信息学的手段,将可能错误的序列过滤掉。
测序仪的下机数据一般都是FASTQ,第二列存放序列,第四列存放对应碱基的质量。由于空间有限,所以无法直接以0.01%这类形式存放概率,必须要做一些转换,从P值先还算成q值
import math
p = 0.001
q = - 10math.log10(p)
# 30
Solexa在1.3版本之前的换算方式是-10log10(p/(1-p)),操作非常的溜
从p值换算成q值,依旧需要两个位置进行存放质量,于是需要进一步编码。这一步采用了ASCII码中字符的位置信息来对应q值。但是ASCII码前面32位是不可见的控制字符,肯定是不能用的,于是就需要往后挪挪,那么挪多少呢?不同测序公司又开始搞自己的一套了。
虽然最后的故事是illumina代表了测序届的半壁江山,格式最后都是Illumina 1.8+,采用Phred+33的形式,但是如果用公共数据的时候一定要小心
# Phred+33
chr(30+33)
#'?'
这一步是根据测序质量对低质量的read进行过滤,Rosalind推荐FASTX-Toolkit,这也是我最早使用的质控工具,但是在使用过程前,我们需要简单的判断下这个测序格式是Phred+33还是Phred+64。
grep 2 rosalind_filt_1_dataset.txt #有结果
grep X rosalind_filt_1_dataset.txt # 无结果
# 基本上断定这个是Phred33
但是编译好的v0.0.13的FASTX-Toolkit的fastq_quality_filter
默认是处理Phred+64,毕竟0.013版本是2012年开发出来的,那个时候主流就是Phred+64。所以只能去自己编译v0.014
题目给我的是p=78,q=24,所以程序按照如下方式运行
tail -n +2 rosalind_filt.txt > rosalind_filt.fq
~/opt/biosoft/fastx_toolkit-0.0.14/bin/fastq_quality_filter -q 24 -p 78 -i rosalind_filt.fq | grep -c '^@Ro'
过滤低质量碱基
如果使用Fastqc发现序列前后几个碱基质量不太好时,我们可以使用trimmomatic
过滤掉按照一定的阈值对read前后进行过滤
问题: 给定一个phred33编码的FASTQ文件,和碱基质量阈值q,给出read前后过滤的文件
解决方法: 用trimmomaitc就行了,例如java -classpath trimmomatic-0.22.jar org.usadellab.trimmomatic.TrimmomaticSE -phred33 data/s1.fq data/tmp.fq TRAILING:30 MINLEN:50
就是过滤前后低于30的碱基,然后删掉不足50的read。
tail -n +2 rosalind_bphr.txt > rosalind_bphr.fq
trimmomatic SE -phred33 rosalind_bfil.fq tt.fq LEADING:22 TRAILING:22