测序数据处理 —— 介绍

最新推荐文章于 2024-09-15 01:19:00 发布

名本无名

最新推荐文章于 2024-09-15 01:19:00 发布

阅读量1k

点赞数 30

分类专栏：测序数据处理文章标签：大数据

本文链接：https://blog.csdn.net/dxs18459111694/article/details/138787625

版权

测序数据处理专栏收录该内容

7 篇文章 0 订阅

订阅专栏

测序数据处理 —— 介绍

前言

有了前面的基础知识铺垫，现在可以正式进入生物信息的领域。首先，我们先介绍目前应用最广泛的基因测序技术。所谓“测序”，“测”字的重点在于怎么测，用什么手段和方法来测，由此发展出了不同的测序方法，常见的基因测序方法包括:

Sanger 测序：也称为合成终止测序，是第一代测序技术，通过不同的 ddNTP（二硫代氨基嘧啶核苷酸）在 DNA 合成时停止链式反应，以测序目标 DNA 序列。但是它对于长片段的测序不太适用。
二代测序（NGS）：这种方法引入了可逆终止末端，从而实现边合成边测序，可以同时测序数百万个 DNA 片段，但所测的片段一般都比较短，是目前应用最广泛的测序技术；
三代测序：可以直接对单个 DNA 分子进行测序，且测序过程无需进行 PCR 扩增，有效避免了因扩增而引入的偏差和错误，提供了更准确和更长的读长，主要包括 PacBio SMRT、Oxford Nanopore、 MGI DNBSEQ T7 等技术；
单细胞测序：通过从单个细胞中分离出 DNA、RNA 等分子，进行高通量测序，以获取单个细胞的基因组、转录组和表观组等信息。
空间转录组测序：空间转录组技术通常使用微流控技术将组织样本分割成小块，并在每个小块上进行转录组测序。这样可以得到每个小块的基因表达数据，并通过图像处理技术将基因表达数据与组织结构关联起来

而“序”字则表明测的对象是什么，是一段段序列，包括 DNA 、RNA 或蛋白质。DNA 测序可用于研究遗传病、肿瘤、微生物学和人类进化等领域；RNA 测序可用于研究基因表达、转录组学和疾病诊断等领域；蛋白质测序技术可用于研究蛋白质的组成、结构和功能，是生物医学研究和临床诊断中重要的工具。

当“测”和“序”作为两个分离的概念时，“测”属于实验层面，如何获取生物数据，属于数据的上游，而“序”则是分析层面，如何分析序列数据并从中挖掘出有价值的信息，属于下游分析。但两者之间却是不可分割的，当我们评估数据质量时，不得不返回去考虑不同实验阶段所可能产生的误差，从数据的产生阶段进行分析，这对于我们开发算法，提高分析结果的质量至关重要。

受限于技术原因，大部分测序方法都是将序列打断成许多长度较短的片段再进行测序，最新的三代测序可以检测更长的序列，实现对单个分子的测序，而不需要打断，但存在通量低、错误率高等问题，还需要改进。抛开一切测序技术不谈，测序的结果都是获取到序列，不管是 DNA、RNA 还是蛋白质序列。所以，我们做下游的数据处理和分析时，首先要面对的便是序列数据，那这些序列数据一般都是什么样子的呢？

序列数据

最常见的测序数据存储格式有 GenBank、FASTA 和 FASTQ，前两种一般代表已知的序列信息，比如人类的参考基因组序列，已知的某个基因或某条染色体序列。而 FASTQ 格式通常用于存储从测序仪下机后的数据。

GenBank

GenBank 是一个开放式的生物序列数据库，收录了来自世界各地的生物序列数据，包括基因组、转录组、蛋白质等，现由美国国家生物技术信息中心 (NCBI) 维护。GenBank 也是一种比较古老的格式了，当初为了人类阅读以及计算机处理的方便，将其设计为固定宽度的格式，前面是个字符串作为一列，用于存储数据的标识字段，后面的所有字符作为一列，用于详细描述字段的信息，其格式大致如下

LOCUS       SCU49845     5028 bp    DNA             PLN       21-JUN-1999
DEFINITION  Saccharomyces cerevisiae TCP1-beta gene, partial cds, and Axl2p
            (AXL2) and Rev7p (REV7) genes, complete cds.
ACCESSION   U49845
VERSION     U49845.1  GI:1293613
KEYWORDS    .
SOURCE      Saccharomyces cerevisiae (baker's yeast)
  ORGANISM  Saccharomyces cerevisiae
            Eukaryota; Fungi; Ascomycota; Saccharomycotina; Saccharomycetes;
            Saccharomycetales; Saccharomycetaceae; Saccharomyces.
REFERENCE   1  (bases 1 to 5028)
  AUTHORS   Torpey,L.E., Gibbs,P.E., Nelson,J. and Lawrence,C.W.
  TITLE     Cloning and sequence of REV7, a gene whose function is required for
            DNA damage-induced mutagenesis in Saccharomyces cerevisiae
  JOURNAL   Yeast 10 (11), 1503-1509 (1994)
  PUBMED    7871890

LOCUS：存储序列的一些元数据信息，包括序列名称、序列长度、分子类型、所属门类以及修改日期
DEFINITION：序列的简要描述，包括来源生物体、基因或蛋白质名称或序列功能的一些描述等信息
ACCESSION：序列的唯一标识符
VERSION：序列标识符后面加上版本号，后面的 GI 是序列标识符，每次更新序列时都会提供一个新的标识
KEYWORDS：描述序列的词或短语。如果条目中不包含关键字，则该字段仅包含一个句点。
SOURCE：自由格式信息，包括生物名称的缩写形式，有时会跟上分子类型标识。
REFERENCE：序列提供者信息及其发表文献

除了上面的基本信息之外，还有序列的特征信息，例如基因的编码区，翻译后的氨基酸序列等

FEATURES             Location/Qualifiers
     source          1..5028
                     /organism="Saccharomyces cerevisiae"
                     /db_xref="taxon:4932"
                     /chromosome="IX"
                     /map="9"
     CDS             <1..206
                     /codon_start=3
                     /product="TCP1-beta"
                     /protein_id="AAA98665.1"
                     /db_xref="GI:1293614"
                     /translation="SSIYNGISTSGLDLNNGTIADMRQLGIVESYKLKRAVVSSASEA
                     AEVLLRVDNIIRARPRTANRQHM"
     gene            687..3158
                     /gene="AXL2"
ORIGIN
        1 gatcctccat atacaacggt atctccacct caggtttaga tctcaacaac ggaaccattg
       61 ccgacatgag acagttaggt atcgtcgaga gttacaagct aaaacgagca gtagtcagct
      121 ctgcatctga agccgctgaa gttctactaa gggtggataa catcatccgt gcaagaccaa
      181 gaaccgccaa tagacaacat atgtaacata tttaggatat acctcgaaaa taataaaccg
      241 ccacactgtc attattataa ttagaaacag aacgcaaaaa ttatccacta tataattcaa
      301 agacgcgaaa aaaaaagaac aacgcgtcat agaacttttg gcaattcgcg tcacaaataa
      361 attttggcaa cttatgtttc ctcttcgagc agtactcgag ccctgtctca agaatgtaat
      421 aatacccatc gtaggtatgg ttaaagatag catctccaca acctcaaagc tccttgccga
      481 gagtcgccct cctttgtcga gtaattttca cttttcatat gagaacttat tttcttattc

GenBank 格式的优势在于它是一种通用格式，可以表示各种信息，同时又保持易读性，但这并不是一种利于数据分析的文件格式。

参考链接：http://www.ncbi.nlm.nih.gov/Sitemap/samplerecord.html。

FASTA

FASTA 格式是一种用于存储 DNA、RNA 和蛋白质序列的文本格式，其中每个核酸或氨基酸都用相应的 ASCII 字符表示，同时还允许在序列之前的描述行中定义序列名称并添加注释。

其基本格式为：描述行 + 序列，例如

>ENST00000619216.1|ENSG00000278267.1|-|-|MIR6859-1-201|MIR6859-1|68|miRNA|
TGTGGGAGAGGAACATGGGCTCAGGACAGCGGGTGTCAGCTTGCCTGACCCCCATGTCGC
CTCTGTAG
>ENST00000410691.1|ENSG00000222623.1|-|-|RNU6-1100P-201|RNU6-1100P|104|snRNA|
ATGCTTGCCTCAGTAGCACACATACTTAAGTTGGAACAATAGAGAGATTGGCACGGCCTC
TGTGAAAGAATGACATGCAAATTTGTGAAGCATTCCATATTTTT

描述行以 > 字符开头，其后紧跟序列名称或唯一标识符，后面还可以添加任意数量的描述信息，包括空格，但必须保证在同一行内。不同数据库对描述行中序列名称格式的定义存在差别，例如， NCBI 的序列格式如下所示

>ref|NC_000019.10|:c41357461-41326782 Homo sapiens chromosome 19, GRCh38.p14 Primary Assembly
AGTAGAAAATGTGGGTTCCTGGAGATGAGTTGATTCTGAGAAATTGCCTAGGTAAATGTGCCATATATTT
CACCAATATAAAGAACGTCTACCTGGTACCTGCCCTGTGTCCAGCCTGTGTTGGGTGATGCCAGGGAGAC
ATTGAGGATCAAGGCAGTGTCCCTTAGCCTTGTCCTCACAGGGTTCACCGTCCTATGGGGGGACACACCT

但是对我们的分析来说其实没那么重要，只要能根据名称找到对应的序列即可。

在注释之后，是由一行或多行构成的序列数据，其中每行的长度应短于 80 个字符。序列可以是核酸或蛋白质序列，还可以包含空白占位或比对用字符。核酸和氨基酸的编码字符要遵循 IUB/IUPAC 标准，支持的核酸编码如下

编码	含义	编码	含义	编码	含义
A	腺嘌呤	T	胸腺嘧啶	U	尿嘧啶
C	胞嘧啶	G	鸟嘌呤	R	嘌呤（A、G）
Y	嘧啶（ C、T、U）	K	酮基（ G、T、U）	M	氨基（ A、C）
S	强结合力（ C、G）	W	弱结合力（ A、T、U）	B	非 A
D	非 C	H	非 G	V	非 U 非 T
N	任意核酸	-	不定长空白占位符

氨基酸编码如下

编码	含义	编码	含义	编码	含义
A（Ala）	丙氨酸	B（Asx）	天冬氨酸或天冬酰胺	C（Cys）	半胱氨酸
D（Asp）	天冬氨酸	E（Glu）	谷氨酸	F（Phe）	苯丙氨酸
G（Gly）	甘氨酸	H（His）	组氨酸	I（Ile）	异亮氨酸
K（Lys）	赖氨酸	L（Leu）	亮氨酸	M（Met）	甲硫氨酸
N（Asn）	天冬酰胺	O（Pyl）	吡咯赖氨酸	P（Pro）	脯氨酸
Q（Gln）	谷氨酰胺	R（Arg）	精氨酸	S（Ser）	丝氨酸
T（Thr）	苏氨酸	U（Sec）	硒半胱氨酸	V（Val）	缬氨酸
W（Trp）	色氨酸	Y（Tyr）	酪氨酸	Z（Glx）	谷氨酸或谷氨酰胺
X	任意	*	转录终止	-	不定长空白占位符

相较于 GenBank 格式，FASTA 更加简单明了，降低了序列操作和分析的难度。需要注意的是，通常核苷酸符号大小写均可，而氨基酸常用大写字母。

参考链接：https://emunix.emich.edu/~mevett/BioinformaticsTools/IUB%20Codes.htm 和 https://en.wikipedia.org/wiki/FASTA_format

FASTQ

FASTQ 是一种用于存储高通量测序数据的文件格式，包括原始测序数据所检测到的每个碱基及其对应的质量分数，该格式每四行表示一条序列（也称为 read）：

第一行以 @ 符号开头，后面跟上序列标识符和可选描述，也可能 barcode 或 UMI 序列等；
第二行是实际检测到的核苷酸序列，一般是 ATCGN，N 表示无法判断是哪个碱基；
第三行以 + 符号开头，后面可选择跟第一行相同的序列标识符，一般这行只有一个 + ；
第四行与第二行一一对应，表示序列中每个碱基的质量分数，。

例如

@A00783:1132:HYV5NDSX3:1:1101:1859:1047 1:N:0:GGACTCCT+TCTACTCT
AGTGAGAGTGTAAGAGGAGAGGCACTCTCATACTTTTCTTGTGGGAATATAAATAGTATAATAATCTCCCTTATTATTTCATA
+
FFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFF:F:FFFFFFFF:FF,FFFFFFFFFFFFFFFFFFFFFFFFFFF:FFF:
@A00783:1132:HYV5NDSX3:1:1101:1877:1047 1:N:0:GGACTCCT+TCTACTCT
AGTGAGAGTGTAAGAGGAGAGGCACTCTCATACTTTTCTTGTGGGAATATAAATAGTATAATAATCTCCCTTATTATTTCATA
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:F

如果是 Illumina 测序平台的话，序列第一行信息一般的形式为

@<instrument>:<run number>:<flowcell ID>:<lane>:<tile>:<x-pos>:<y-pos>:<UMI> <read>:<is filtered>:<control number>:<index>

对于 FASTQ 数据，我们还需要理解第四行字符是如何表示碱基质量的。测序仪在进行测序时，会自动根据荧光信号的强度计算每个碱基被测错的概率。概率值是一个小数，但是小数不适合存储，所以先对概率值 P 进行对数转换
$Q=round(-10*log_{10}(P))$
比如，概率为 0.01，那算出的 Q 值就是 20，但一般整数也要占用 4 个字节，如果能转换为字符，只占一个字节，可以大大减少存储的开销。所以将这个 Q 值加上 33 或 64 将其转换为 ASCII 字符，就可以实现
$P h re d = Q + 33$
比如 Q 值为 20 对应的字符是 5。对应关系大致如下所示

!"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHI 
|    |    |    |    |    |    |    |    | 
0....5...10...15...20...25...30...35...40 
|    |    |    |    |    |    |    |    | 
worst <-...........................-> best

值越大表示测序质量越好