判断fastq文件是Phred33还是Phred64编码-脚本02

很多从NCBI上下载的公共数据，都面临数据质量控制以及过滤等问题。在此面临一个问题，这个fastq文件质量编码是什么，如何根据质量编码过滤数据？

01 FASTQ文件简介

FASTQ文件是一种文本格式，常用于存储高通量测序（NGS）数据。每个读取序列（read）包括四行信息，具体格式如下：

1. 第一行：以`@`开头，后面跟着读取序列的ID以及一些可选的描述信息。
2. 第二行：读取的碱基序列（A、T、C、G）。
3. 第三行：以`+`开头，后面可以是空的或附带注释。
4. 第四行：质量得分的编码，表示碱基读取的准确性。

seqkit安装与使用 v2.5.1（生物信息学工具-003）

@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGT
+
!''*((((***+))%%%++)(%%%%).1***-+*''

在FASTQ文件中，第四行的符号表示每个碱基的质量得分，这些符号通过ASCII码进行编码，质量得分的计算方式基于Phred质量评分。

02 Phred质量评分（Phred Quality Score）

Phred质量评分（`Q`值）用于表示碱基读取的准确性，具体公式如下：

Q = -10 * log10(P)

其中，`P`是错误读取碱基的概率。较高的`Q`值代表更低的错误概率。

| Phred值 | 错误概率 | 正确率 |
|---------|----------|--------|
| Q10     | 1/10     | 90%    |
| Q20     | 1/100    | 99%    |
| Q30     | 1/1000   | 99.9%  |
| Q40     | 1/10000  | 99.99% |

Namely，质量值Q是测序错误率的对数*-10。假如错误率是0.01，则Q值为20。可见，错误率越低，其Q值越高。即Q值越高越可靠。

Q10——0.1
Q20——0.01
Q30——0.001
Q40——0.0001

03 Phred编码：Phred33与Phred64

Phred编码将质量得分转化为ASCII字符进行存储，主要有两种常见的编码方式：Phred33和Phred64。

3.1. Phred33编码

范围：ASCII字符33到126（可打印字符的范围）。
字符映射：`Q`值加33后对应ASCII字符。即`Q=0`对应的ASCII字符是`!`，`Q=40`对应的ASCII字符是`I`。
使用：Phred33编码是目前常用的编码方式，Illumina平台在2011年后的数据基本都使用这种编码。 

3.2. Phred64编码

范围：ASCII字符64到126。
字符映射：`Q`值加64后对应ASCII字符。即`Q=0`对应的ASCII字符是`@`，`Q=40`对应的ASCII字符是`h`。
使用：Phred64编码较少使用，主要在早期的Illumina数据中使用，后来被Phred33编码所取代。 

04 如何区分Phred33和Phred64

可以通过查看质量得分的最小和最大ASCII字符来判断是使用Phred33还是Phred64编码。常见的步骤如下：

1. 提取质量得分：从FASTQ文件中提取质量得分的第四行。
2. 转换为ASCII值：将质量得分的每个字符转换为相应的ASCII码值。
3. 判断范围：
   - 若ASCII值范围在33到73之间，则为Phred33编码。
   - 若ASCII值范围在64到104之间，则为Phred64编码。

05 常见质量控制工具

FastQC：一个常用工具，用于检查FASTQ文件的质量、检测Phred编码方式、评估序列的碱基质量分布等。
Seqtk：可以用于转换和修复FASTQ文件，包括编码格式的转换。

seqtk安装与使用-seqtk-1.4（bioinfomatics tools-012）

5.1 Phred编码转换示例

假设你有一个Phred64编码的文件，并且需要将其转换为Phred33编码，可以使用以下命令（以`seqtk`工具为例）：

seqtk seq -Q64 -V file_phred64.fastq > file_phred33.fastq

06 总结

FASTQ文件是存储NGS数据的标准格式，每个读取序列包含四行：序列ID、碱基序列、分隔符和质量得分。
Phred质量评分用于表示碱基读取的准确性，分数越高，表示读取的准确性越高。
Phred33 和 Phred64 是两种不同的编码方式，现代测序数据大多使用Phred33编码。

07 一个shell脚本快速识别

vi script_qc.sh

#!/bin/bash

# Check if a file is passed as an argument
if [ -z "$1" ]; then
    echo "Usage: $0 <filename>"
    exit 1
fi

# Ensure the file exists
if [ ! -f "$1" ]; then
    echo "File not found!"
    exit 1
fi

# Processing the file to determine the quality score encoding
less "$1" | head -n 1000 | awk '{if(NR % 4 == 0) printf("%s", $0);}' \
| od -A n -t u1 -v \
| awk '
BEGIN { min=100; max=0; }
{
    for(i=1; i<=NF; i++) {
        if($i > max) max=$i;
        if($i < min) min=$i;
    }
}
END {
    if(max <= 126 && min < 59) 
        print "Phred33";
    else if(max > 73 && min >= 64) 
        print "Phred64";
    else if(min >= 59 && min < 64 && max > 73) 
        print "Solexa64";
    else 
        print "Unknown score encoding";
}'

运行

bash ./script_qc.sh filename.fastq

即可判断是Phred33还是Phred64编码，一般是Phred33！