短读质量控制 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