计算生物学与生物信息学漫谈-2-测序深度/读长质量和Fasta处理

上一篇文章中我们介绍了测序技术的由来与发展，那么在介绍第三代测序的时候，我们提到了关于测序深度和读长的问题，那么本篇文章就详解介绍一下。

计算生物学与生物信息学漫谈-1-测序一路走来-CSDN博客

---

目录

1.测序深度SEQUENCING DEPTH

（1）reads

（2）Coverage

2.Base Call Quality

（1）什么是call

（2）Phred quality score

3.FASTQ 文件

---

1.测序深度SEQUENCING DEPTH

（1）reads

受测序水平限制，测序时需先将基因组打断成DNA片段，然后再建库测序。reads（读长）指的是测序仪单次测序所得到的碱基序列，也就是一连串的ATCGGGTA之类的，这些序列并不是基因组的完整组成部分，而是通过高通量测序技术从基因组中获取的短序列片段。不同测序仪器产生的 reads 长度可能有所不同。不同的测序仪器，reads长度不一样。对整个基因组进行测序，就会产生成百上千万的reads。

（2）Coverage

不同的测序应用中，生物结果和测序数据的解读在很大程度上受到覆盖基因组区域的测序读数数量的影响。通常，多个序列会在基因组的某些区域重叠。测序深度衡量的是平均读数丰度，计算方法是将所有与基因组匹配的测序短读数的碱基数除以该基因组的长度（如果已知基因组大小）。

如果reads读数长度相等，则测序覆盖度（测序深度）计算公式为：

如果reads读书长度不相等，则测序覆盖度计算为：

n是reads数。

测序覆盖率表示为基因组被测序的次数（例如，1X、2X、20X等）。 测序深度影响基因组组装的完整性、从头组装和参考引导组装的准确性、检测到的基因数量、RNA-Seq中的基因表达水平、变异调用、全基因组测序中的基因分型、宏基因组学中的微生物鉴定和多样性分析，以及表观遗传学中蛋白质-DNA相互作用的识别。因此，在进行序列分析之前，研究测序深度非常重要。碱基被测序的次数越多，数据的质量就越好。

2.Base Call Quality

（1）什么是call

在基因测序和生物信息学领域，“Base Call Quality”中的“call”指的是碱基判读，即从测序过程中获得的原始数据中识别出每个碱基（A、T、C、G）的过程。

（2）Phred quality score

得到碱基在序列中位置的过程称为base calling，我们之前讲过的各种测序，无非就是想得到序列的碱基顺序，但是目前有的方法及测序的仪器和样本，都会造成最后的结果存在错误，所以许多base calling的软件中，会计算Phred Quality Score来量化发生错误的可能性，且这个指标把难以量化的可能性转变为了数字参数：

其中 p 是碱基调用出错的概率。

Phred质量分数使用ASCII单个字符进行编码。所有ASCII字符都有一个与之关联的十进制数字。然而，由于前32个ASCII字符是非打印字符，而整数33是惊叹号“!”的十进制数字，因此Q=0就是惊叹号，并且以“!”开始的编码称为Phred+33编码。

Illumina 1.8及更高版本使用这种Phred+33编码（Q33）来在FASTQ文件中编码base calling的质量。较旧的Illumina版本（例如Solexa）使用Phred+64编码，在这种编码中，字符“@”（其十进制数字为64）对应于Q=0。

表格中显示了Phred质量分数（Q）、相应的概率（P）以及十进制数字和ASCII代码。例如，当调用碱基的概率为0.1时，Phred分数将是10（Q=10），但不是给出数字10，而是将该质量分数编码为加号“+”。

较高的Q分数表示错误概率较小，而较低的Q分数表示碱基调用的质量较低，更有可能碱基被错误地调用。例如，质量分数为20表示在100次中有1次出错的概率（1%的错误率），相当于99%的调用准确性。一般来说，Q分数为30被认为是高通量测序（HTS）中良好质量的基准。第二章表显示了一些Q分数及其对应的错误概率、碱基调用准确性和解释。

3.FASTQ 文件

像Illumina这样的测序技术提供了实时分析（RTA）软件，该软件将单个碱基调用数据存储在称为BCL文件的中间文件中。当测序运行完成后，这些BCL文件会被过滤，如果样品是多重化的还会进行解复用，然后转换成名为FASTQ的序列文件格式。

对于单端运行的每个样本，将有一个FASTQ文件；而对于双端运行的每个样本，则会有两个FASTQ文件（R1和R2）：R1文件用于正向读取，R2文件用于反向读取。 FASTQ文件通常会被压缩，它们可能具有“.fastq.gz”的文件扩展名。FASTQ是一种可读的文件格式，已成为大多数高通量测序（HTS）技术输出存储的 facto标准。

一个FASTQ文件由多个记录组成，每个记录包含四行数据，如图所示。

FASTQ文件中每个记录的第一行以“@”符号开始，这一行被称为读取标识符，因为它标识了序列（读取）。

一个典型的由Illumina仪器生成的读取的FASTQ标识符行如下所示：

字符参数	描述
@	读标识符行的开始
<instrument>	设备ID或序列ID
<run num>	仪器运行的次数
<flowcell ID>	流动池ID
<lane>	读取序列所在的泳道编号
<tile>	读取序列所在的tile编号
<x>	DNA簇的X坐标
<y>	DNA簇的Y坐标
<UMI>	仅当使用唯一的分子标识符（UMI）时
<read>	读取编号（单端读取为1或双端读取为2）
<filtered>	如果读取通过了过滤则为Y，如果没有通过则为N
<control num>	0（没有任何控制位被激活）或偶数

表中描述了Illumina FASTQ标识符行的元素，而上图显示了一个包含三个读取记录的示例FASTQ文件。在索引序列中观察到的序列会被写入FASTQ标题以代替样本编号。这些信息对于故障排除和解复用非常有用。然而，这些元数据元素可能会被其他元素修改或替换，特别是在提交到数据库或被用户修改时。 FASTQ文件的第二行包含了测序仪推断出的碱基。这些碱基包括腺嘌呤、胞嘧啶、鸟嘌呤和胸腺嘧啶，分别用A、C、G和T表示。如果某个位置的碱基不明确（由于测序错误而未确定），则可能会包含字符N。 第三行以加号“+”开始，它可能包含其他附加的元数据或相同的标识符行元素。

FASTQ文件的第四行包含ASCII编码的字符串，代表每个碱基的Phred质量分数。每个ASCII字符的数值对应于序列行中碱基的质量分数。 研究人员通常从测序仪器获取原始测序数据用于自己的研究。原始测序数据也可以从数据库下载，科学家和研究机构会将自己的原始数据存档并公开提供。无论哪种情况，原始测序数据通常都是以FASTQ文件的形式获得的。

NCBI SRA数据库是数百种物种原始数据的最大数据库之一。FASTQ文件以序列读取档案（SRA）格式存储，可以使用SRA-toolkit下载和提取，这是由NCBI开发的一系列程序集合，可以从“https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi”提供的说明下载和安装。

主页 - SRA - NCBI --- Home - SRA - NCBI (nih.gov)

为了演示的目的，我们将从NCBI SRA数据库下载原始数据。我们将使用一个单端FASTQ文件，其运行ID为“SRR030834”

该FASTQ文件包含从格陵兰岛Qeqertasussuk出土的4000年前已灭绝的Saqqaq古爱斯基摩人头发束中测序得到的reads。为了保持文件组织性，您可以创建目录“fastqs”，然后使用“fasterq-dump”命令下载FASTQ文件（确保您已在计算机上安装了SRA-toolkit，并且它在路径上）：

Linux服务器安装SRAToolkit教程-CSDN博客

 mkdir fastqscd fastqsmkdir singlecd singlefasterq-dump --verbose SRR030834

 FASTQ文件可能包含多达数百万条目，其大小可能是几兆字节或吉字节，这通常使它们太大而无法在普通文本编辑器中打开。一般来说，除非有必要进行故障排除或出于好奇，否则无需打开FASTQ文件。要显示大型FASTQ文件，我们可以使用一些Unix或Linux命令，如“less”或“more”来逐页显示非常大的文本文件，或使用“cat”来显示文件的内容。

Bio-Linux-shell详解-2-基本Shell命令快速掌握-CSDN博客

如果FASTQ文件名以“.gz”扩展名结尾，这意味着该文件已使用“gzip”程序压缩。在这种情况下，应分别使用“zless”、“zmore”和“zcat”命令，而不是“less”、“more”和“cat”命令，且无需解压缩文件。

我们还可以使用“head”和“tail”分别显示文件的前几行和最后几行。以下命令将显示文件的前15行：

head -15 SRR030834.fastq

 如果FASTQ文件很大，我们可以使用“gzip”程序将其压缩，以使其大小减少三倍以上。使用gzip压缩“SRR030834.fastq”文件将使其大小减少到不到一G字节。

gzip SRR030834.fastq

使用“gzip -d”可以解压一个已压缩的文件。这个命令会移除原始的压缩文件（例如 .gz 文件），并生成一个解压后的文件。如果你想要保留原始的压缩文件，可以使用“gunzip -c”或者“zcat”命令，并将输出重定向到一个新的文件中。

gzip -d SRR030834.fastq.gz

如果你需要知道FASTQ文件中的记录数，可以使用“cat”或“zcat”与“wc -l”的组合，后者用于计算文本文件中的行数。请记住，FASTQ文件中的一条记录包含4行。

我们可以使用Unix管道符号“|”将“cat”命令的输出传递给“wc -l”命令。以下命令行将计算存储在FASTQ文件中的记录数：

cat SRR030834.fastq | echo $((`wc -l`/4))

如果需要显示目录中多个以“.fastq”为文件扩展名的文件的文件名和读取计数，我们可以使用以下脚本：

 for filename in *.fastq; 
doecho -e “$filename\t `cat $filename | wc -l | awk ‘{print $1 / 
4}’`”done

要以表格格式显示FASTQ文件，您可以使用“cat”命令，然后使用Unix管道将输出传递给“paste”命令，该命令将FASTQ记录的四行转换为表格格式。

cat SRR030834.fastq | paste - - - - > SRR030834_tab.txt

 less -S SRR030834_tab.txt

从FASTQ文件创建表格文件将帮助我们执行多种操作，例如排序条目、过滤掉重复的读取、提取读取ID、序列或质量分数，以及创建FASTA文件。我们期望FASTQ文件的标识符行的格式是一致的。如果您显示“SRR030834tab.txt”，您会注意到标识符行的某些字段是由空格分隔的，如果我们认为空格是列分隔符，那么ID将在第一列，序列将在第四列。然而，在从其他FASTQ文件提取的表格文件中，这个列顺序可能不同。假设我们只想从“SRR030834tab.txt”中提取ID和序列到一个单独的文本文件中，那么我们可以使用“awk”命令如下：

 awk ‘{print $1 “\t” $4}’ SRR030834_tab.txt > SRR030834_seq.txt

“awk”命令从“SRR030834tab.txt”中提取第一列和第四列，并打印出这两列，它们之间用制表符分隔。输出被定向到一个新的文本文件“SRR030834seq.txt”。

Linux命令允许我们执行多步操作。假设我们想从FASTQ文件创建一个FASTA文件；我们可以通过多个步骤来实现。首先，我们需要像上面那样提取ID和序列到一个文件中，然后我们可以移除“@”符号，只留下ID，接着我们需要在每一行的开头添加“>”，且“>”和ID之间没有空格，最后，我们将两列分开，形成FASTA的定义行（defline）和序列，将它们存储在一个文件中，并删除临时文件。

cat SRR030834.fastq | paste - - - - \> SRR030834_tab.tmpawk ‘{print $1 “\t” $4}’ SRR030834_tab.tmp \| sed ‘s/@//g’ > SRR030834_seq.tmpsed -i ‘s/^/>/’ SRR030834_seq.tmpawk ‘{print $1, “\n” $2}’ SRR030834_seq.tmp \> SRR030834.fastarm *.tmp

---

以上就是这次的全部内容，下一次将介绍使用FastaQC工具进行Fasta文件过滤与质控。

有任何问题欢迎与我联系。