计算生物学与生物信息学漫谈-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文件过滤与质控。

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

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/452603.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

现代物流管理:SpringBoot技术突破

3系统分析 3.1可行性分析 通过对本智能物流管理系统实行的目的初步调查和分析&#xff0c;提出可行性方案并对其一一进行论证。我们在这里主要从技术可行性、经济可行性、操作可行性等方面进行分析。 3.1.1技术可行性 本智能物流管理系统采用SSM框架&#xff0c;JAVA作为开发语…

【SQL|大数据|数据清洗|过滤】where条件中 “ != “ 和 “ NOT IN() ” 对NULL的处理

对数据进行清洗过滤的时候&#xff0c;NULL往往是一个很特殊的存在&#xff0c;对NULL值的存在通常有以下三种方式 1、保留NULL 2、过滤掉NULL 3、将NULL替换为其他符合业务需求的默认常量 下面是一些常用处理NULL的方式&#xff1a; 如下图所示数据源&#xff1a; car_vin&…

嵌入式入门学习——6Protues点亮数码管,认识位码和段码,分辨共阴还是共阳(数字时钟第一步)

0 系列文章入口 嵌入式入门学习——0快速入门&#xff0c;Let‘s Do It&#xff01; 首先新建基于Arduino UNO的protues工程&#xff0c;见本系列第3篇文章 1 点“P”按钮找器件 2 输入“seg”或“digit”查找数码管器件 3 找到我们想要的6位7段数码管 4如图A、B…DP都是段码…

ESP32-C3 入门笔记04:gpio_key 按键 (ESP-IDF + VSCode)

1.GPIO简介 ESP32-C3是QFN32封装&#xff0c;GPIO引脚一共有22个&#xff0c;从GPIO0到GPIO21。 理论上&#xff0c;所有的IO都可以复用为任何外设功能&#xff0c;但有些引脚用作连接芯片内部FLASH或者外部FLASH功能时&#xff0c;官方不建议用作其它用途。 通过开发板的原…

自由学习记录(11)

Surface Effector 2D Platform Effector 2D 向上跳跃穿过天花板的功能 平台效应器不用变Trigger&#xff0c;因为本来就是要有碰撞的 use one way grouping是让所有相关联的碰撞器都可以单面跳墙 默认不勾选&#xff0c;左右两边没有摩擦力和弹力&#xff0c;要自己先设置sid…

三菱PLC伺服-停止位置不正确故障排查

停止位置不正确时&#xff0c;请确认以下项目。 1)请确认伺服放大器(驱动单元)的电子齿轮的设定是否正确。 2&#xff09;请确认原点位置是否偏移。 1、设计近点信号(DOG)时&#xff0c;请考虑有足够为0N的时间能充分减速到爬行速度。该指令在DOG的前端开始减速到爬行速度&…

计算机毕业设计 基于Web的景区管理系统的设计与实现 Java+SpringBoot+Vue 前后端分离 文档报告 代码讲解 安装调试

&#x1f34a;作者&#xff1a;计算机编程-吉哥 &#x1f34a;简介&#xff1a;专业从事JavaWeb程序开发&#xff0c;微信小程序开发&#xff0c;定制化项目、 源码、代码讲解、文档撰写、ppt制作。做自己喜欢的事&#xff0c;生活就是快乐的。 &#x1f34a;心愿&#xff1a;点…

市场上几个跨平台开发框架?

跨平台桌面应用开发框架是一种工具或框架&#xff0c;它允许开发者使用一种统一的代码库或语言来创建能够在多个操作系统上运行的桌面应用程序。传统上&#xff0c;开发者需要为每个操作系统编写不同的代码&#xff0c;使用不同的开发工具和语言。而跨平台桌面应用开发框架通过…

【网络】IP协议的地址管理

【网络】IP协议的地址管理 一. IP协议格式二. 地址管理1.动态分配IP地址2.NAT机制2.1 NAT机制下网络的请求/响应 3. 网段划分3.1 特殊的IP地址 4.路由选择5.DNS域名解析系统 一. IP协议格式 4位版本号(version): 指定IP协议的版本&#xff08;IPv4/IPv6&#xff09;, 对于IPv4来…

一起搭WPF架构之livechart的MVVM使用介绍

一起搭WPF架构之livechart使用介绍 前言ModelViewModelView界面设计界面后端 效果总结 前言 简单的架构搭建已经快接近尾声了&#xff0c;考虑设计使用图表的形式将SQLite数据库中的数据展示出来。前期已经介绍了livechart的安装&#xff0c;今天就详细介绍一下livechart的使用…

第6篇:无线与移动网络

目录 引言 6.1 无线网络的基础概念 6.2 无线局域网&#xff08;WLAN&#xff09;与IEEE 802.11 6.3 蓝牙与无线个域网&#xff08;WPAN&#xff09; 6.4 无线城域网&#xff08;WMAN&#xff09;与WiMax 6.5 ZigBee与智能家居 6.6 移动蜂窝网络&#xff08;3G/4G/5G&…

SpringCloudStream使用StreamBridge实现延时队列

利用RabbitMQ实现消息的延迟队列 一、安装RabbitMQ 1、安装rabbitmq 安装可以看https://blog.csdn.net/qq_38618691/article/details/118223851,进行安装。 2、安装插件 安装完毕后,exchange是不支持延迟类型的,需要手动安装插件,需要和安装的rabbitmq版本一致 https:…

深入探讨C++多线程性能优化

深入探讨C多线程性能优化 在现代软件开发中&#xff0c;多线程编程已成为提升应用程序性能和响应速度的关键技术之一。尤其在C领域&#xff0c;多线程编程不仅能充分利用多核处理器的优势&#xff0c;还能显著提高计算密集型任务的效率。然而&#xff0c;多线程编程也带来了诸…

62天框架安全(学习)

发现学了之后没有去复习&#xff0c;每天都要问自己学了什么&#xff0c;复习了吗&#xff0c;下次还能记住吗 一下内容来自【小迪安全2023】第62天:服务攻防-框架安全&CVE复现&Spring&Struts&Laravel&ThinkPHP_小迪安全文档2023-CSDN博客 一个网站的源码…

排序算法 —— 直接插入排序

目录 1.直接插入排序的思想 2.直接插入排序的实现 实现分析 实现代码 3.直接插入排序的分析 时间复杂度分析 空间复杂度分析 稳定性 1.直接插入排序的思想 直接插入排序的思想就是把待排序的元素按其关键码值的大小依次插入到一个已经排好序的有序序列中&#xff0c…

一种基于OCR图像识别技术的发票采集管理系统及方法

一种基于OCR图像识别技术的发票采集管理系统及方法 摘要 本发明涉及了一种基于OCR图像识别技术的发票采集管理系统及方法&#xff0c;该系统的发票信息采集单元采集发票图片信息数据&#xff0c;OCR图像识别单元基于OCR图像识别技术并结合人工智能深度学习算法对发票图片信息数…

vscode默认添加python项目的源目录路径到执行环境(解决ModuleNotFoundError: No module named问题)

0. 问题描述 vscode中编写python脚本&#xff0c;导入工程目录下的其他模块&#xff0c;出现ModuleNotFoundError: No module named 错误 在test2的ccc.py文件中执行print(sys.path) 查看路径 返回结果发现并无’/home/xxx/first_demo’的路径&#xff0c;所以test2下面的文…

Vue-router 路由守卫执行流程图

vue-router 路由守卫执行的流程图&#xff08;个人理解&#xff09; 图1 - 图2

【MR开发】在Pico设备上接入MRTK3(一)——在Unity工程中导入MRTK3依赖

写在前面的话 在Pico上接入MRTK3&#xff0c;目前已有大佬开源。 https://github.com/Phantomxm2021/PicoMRTK3 也有值得推荐的文章。 MRTK3在PICO4上的使用小结 但由于在MacOS上使用MRTK3&#xff0c;无法通过Mixed Reality Feature Tool工具管理MRTK3安装包。 故记录一下…

jmeter使用文档

文章目录 一、安装使用1、下载2、bin/jmeter.properties介绍 二、windows使用1、微调&#xff08;1&#xff09;界面样式&#xff08;2&#xff09;修改语言 2、简单使用3、各组件详解&#xff08;1&#xff09;CSV 数据文件配置&#xff08;2&#xff09;BeanShell取样器 三、…