1. 入门课程
image.png
2. 测序原理
（待填坑）
3. 分析流程
image.png
第1题，与FASTQ与FASTA格式有关
1.0 掌握fasta格式
-
概述一下，fasta格式是一种非常简单的储存序列的格式，可以储存核酸序列（DNA/RNA）也可以储存蛋白质的氨基酸序列（Amino Acid sequence，简称AA序列），主要分成2个部分。
-
举个例子
1. >sp|P69905|HBA_HUMAN Hemoglobin subunit alpha OS=Homo sapiens GN=HBA1
2. MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSHGSAQVKGHG
3. KKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLLSHCLLVTLAAHLPAEFTP
4. AVHASLDKFLASVSTVLTSKYR
-- 第1行，以“>”开始，存储序列的描述信息
-- 第2行，序列部分，中间前后都可以有空格。序列部分一般在70-80之间，小于120.
-- 其实实际操作中，程序处理的时候都是自动去掉空格和换行符，把序列读成1行再处理，所以，我也干过把整条人类染色体都放到一行的233举动，这么算下来，一行可以有240*10E6这么长
- 再举个例子-蛋白序列
1. >sp|P69905|HBA_HUMAN Hemoglobin subunit alpha OS=Homo sapiens GN=HBA1
2. MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSHGSAQVKGHG
3. KKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLLSHCLLVTLAAHLPAEFTP
4. AVHASLDKFLASVSTVLTSKYR
这个例子是从UniRef数据库中下载的人类血红蛋白α亚基的序列。
1. >sp|P69905|HBA_HUMAN Hemoglobin subunit alpha OS=Homo sapiens GN=HBA1
信息解读：
-- 第一行，
- sp代表该蛋白是在Swiss-Prot数据库（该数据库中的蛋白经过文献+实验核实过，可信度高）
- P69905是该序列在uniref数据库中的编号
- HBA_HUMAN 是这个序列的简称
- Hemoglobin subunit alpha是全称
- OS=Homo sapiens是物种
- GN=HBA1是指基因的名字是HBA1
后面的内容就是这个HBA1基因对应的蛋白序列。
字母代表氨基酸种类。
A alanine P proline
B aspartate/asparagine Q glutamine
C cystine R arginine
D aspartate S serine
E glutamate T threonine
F phenylalanine U selenocysteine
G glycine V valine
H histidine W tryptophan
I isoleucine Y tyrosine
K lysine Z glutamate/glutamine
L leucine X any
M methionine * translation stop
N asparagine - gap of indeterminate length
- 再举个例子-核酸序列
为了方便大家理解，我们继续使用人类血红蛋白a亚基对应的mRNA序列，这个序列是从NCBI RefSeq数据库中下载的。
>gi|13650073|gb|AF349571.1| Homo sapiens hemoglobin alpha-1 globin chain (HBA1) mRNA, complete cds
CCCACAGACTCAGAGAGAACCCACCATGGTGCTGTCTCCTGACGACAAGACCAACGTCAAGGCCGCCTGG
GGTAAGGTCGGCGCGCACGCTGGCGAGTATGGTGCGGAGGCCCTGGAGAGGATGTTCCTGTCCTTCCCCA
CCACCAAGACCTACTTCCCGCACTTCGACCTGAGCCACGGCTCTGCCCAGGTTAAGGGCCACGGCAAGAA
GGTGGCCGACGCGCTGACCAACGCCGTGGCGCACGTGGACGACATGCCCAACGCGCTGTCCGCCCTGAGC
GACCTGCACGCGCACAAGCTTCGGGTGGACCCGGTCAACTTCAAGCTCCTAAGCCACTGCCTGCTGGTGA
CCCTGGCCGCCCACCTCCCCGCCGAGTTCACCCCTGCGGTGCACGCCTCCCTGGACAAGTTCCTGGCTTC
TGTGAGCACCGTGCTGACCTCCAAATACCGTTAAGCTGGAGCCTCGGTGGCCATGCTTCTTGCCCCTTTG
G
对应第一行：
- gi|13650073：所有来自于NCBI的序列都有一个gi号，就是数据库中的流水号，具有唯一性。
- gb|AF349571.1：是genebank编号的信息。
- Homo sapiens hemoglobin alpha-1 globin chain (HBA1) mRNA, complete cds：后面就是对序列的一个通俗的描述。这里使用的是mRNA，包含完整的CDS（Coding sequence）区域。
疑问：为什么mRNA数据还是用T而不用U表示？
解释：其实这是为了保证数据的统一性，因为U只是在RNA中替换了原来的T，所以为了下游的方便分析处理，无论RNA序列还是DNA序列都是使用T而不是U。
** 对于核苷酸序列的信息，我们一般使用下面的对应表。**
A adenosine C cytidine G guanine
T thymidine N A/G/C/T (any) U uridine
K G/T (keto) S G/C (strong) Y T/C (pyrimidine)
M A/C (amino) W A/T (weak) R G/A (purine)
B G/T/C D G/A/T H A/C/T
V G/C/A - gap of indeterminate length
基本上FASTA就介绍到这里，这个格式主要是把序列储蓄到数据库中的一种格式，但是它不适合储存我们刚刚测到的测序数据。一个很重要的原因就是，它没有序列的质量信息。
那一般带有测序质量信息的FASTQ格式就成了储存测序数据的常用格式啦！
1.1 掌握FASTQ格式
1)格式有什么特点？
下面是一个Illumina平台测序的真实数据，其中包含了1条reads的信息。
1. @ST-E00126:128:HJFLHCCXX:2:1101:7405:1133
2. TTGCAAAAAATTTCTCTCATTCTGTAGGTTGCCTGTTCACTCTGATGATAGTTTGTTTTGG
3. +
4. FFKKKFKKFKF<KK<F,AFKKKKK7FFK77<FKK,<F7K,,7AF<FF7FKK7AA,7<FA,,
fastq内容格式有4行：
-- 第1行：主要储存序列测序时的坐标等信息；
举个例子：
@ST-E00126:128:HJFLHCCXX:2:1101:7405:1133
1. @，开始的标记符号;
2. ST-E00126:128:HJFLHCCXX，测序仪唯一的设备名称;
3. 2，lane的编号，说明位于flowcell的第2个lane；
4. 1101，tile的坐标，代表了第2个lane中的第1101个tile
5. 7405:1133，该read在第2个lane中第1101个tile中的什么位置，该tile中的X坐标是7405，该tile中的Y坐标是1133
-- 第2行：就是测序得到的序列信息，一般用ATCGN来表示，其中N用于荧光信号干扰无法判断到底是哪个碱基时的代表符号；
-- 第3行：以“+”开始，可以储存一些附加信息，但目前的测序fastq文件这一行一般是空的，但是“+”不可省略；
-- 第4行：储存的是质量信息，与第2行的碱基序列是一一对应的，其中的每一个符号对应的ASCII值是经过换算的phred值，可以简单理解为对应位置碱基的测序质量值，越大说明测序的质量越好。不同的版本对应的phred值范围不同。
详细谈谈FASTQ质量值的计算方法
在测序仪进行测序的时候，会自动根据荧光信号的强弱给出一个参考的测序错误概率（error probility，P）根据定义来说，P值肯定是越小越好。我们怎么储存他们呢？直接储存成小数点？比如1%储存程0.01？这肯定是不高效的，因为一个碱基就得占用4个字符。
解决办法是什么呢？
科学家想的办法是：
1.将P取log10之后再乘以-10，得到的结果为Q。
1. 比如，P=1%，那么对应的Q=-10*log10（0.01）=20
2.把这个Q加上33或者64转换成一个新的数值，称为Phred，最后把Phred对应的ASCII字符对应到这个碱基。
Phred值简单来说就是：reads质量得分
Fastq格式里的reads质量得分编码方式有好几种，现在Illumina用的一般是Phred33，但偶尔还会遇到Phred64（旧版本）的。
1. 如Q=20，Phred = 20 + 33 = 53，对应的符号是”5”
这样就可以用1个符号与1个碱基一一对应，是不是很聪明？
各版本不同Phred对应的ASCII值
-
在计算Q值和加上33/64的时候，不同测序仪，产生的数据不同，大概如下所示：
-- Solexa标准
image.png
-- Illumina标准
image.png
不同测序仪的不同Phred值对应的ASCII表
image.png
说明：
1. 对于每个碱基的质量编码标示，不同的软件采用不同的方案，目前有5种方案：
2. Sanger，Phred quality score，值的范围从0到92，对应的ASCII码从33到126，但是对于测序数据（raw read data）质量得分通常小于60，序列拼接或者mapping可能用到更大的分数。
3. Solexa/Illumina 1.0, Solexa/Illumina quality score，值的范围从-5到63，对应的ASCII码从59到126，对于测序数据，得分一般在-5到40之间；
4. Illumina 1.3+，Phred quality score，值的范围从0到62对应的ASCII码从64到126，低于测序数据，得分在0到40之间；
5. Illumina 1.5+，Phred quality score，但是0到2作为另外的标示，详见http://solexaqa.sourceforge.net/questions.htm#illumina
6. Illumina 1.8+
知道了通过测序之后我们拿到的序列信息是通过FASTQ格式储存的。那么我们怎么衡量高通量测序中，1次测序中测序质量到底是怎么样的？我们应该对数据进行哪些修正？这个时候，我们就需要介绍下一期的神器FastQC工具！
2) 什么是phred值，怎么计算？
是评估这个bp测序质量的值，测序仪通过判断荧光信号的颜色来判断碱基的种类，ATCG分别对应红黄蓝绿，信号强弱不同，在这种情况下对每个结果的判断的正确性都存在一个概率值，这个值被储存为ASCII码形式，转化方式如下：
- 将该碱基判断错误概率值P取log10之后再乘以-10，得到的结果为Q。
比如，P=1%，那么对应的Q=-10*log10（0.01）=20（这个计算公式illumina平台使用，Solexa系列测序仪使用不同的公示来计算质量值：Q=-10log(P/1-P)）
- 把这个Q加上33或者64转成一个新的数值，称为Phred，最后把Phred对应的ASCII字符对应到这个碱基。
如Q=20，Phred = 20 + 33 = 53，53在ASCII码表里对应的ASCII符号是”5”
3) phred33 与 phred64是什么意思？
P ---> Q + 33/64 = Phred33/64 ---> ASCII
P--在测序仪进行测序的时候，会自动根据荧光信号的强弱给出一个参考的**测序错误概率**（error probility，P）;
Q--将该碱基判断错误概率值P取log10之后再乘以-10，得到的结果为Q；比如，P=1%，那么对应的Q=-10*log10（0.01）=20（这个计算公式illumina平台使用，Solexa系列测序仪使用不同的公示来计算质量值）
Phred33/64--把这个Q加上33或者64转成一个新的数值，称为Phred，最后把Phred对应的ASCII字符对应到这个碱基的测序错误概率。
质量字符的ASCII值和质量得分的关系有如下两种：可以粗略分为 Phred+33和Phred+64，这里的33和64就是指ASCII值转换为Q该减去的数值。
在处理测序数据时，因为一些软件会根据碱基质量得分的不同做不同的处理，常要指定正确的编码方式，有必要对质量字符与质量得分的关系（Phred+33或Phred+64）作出正确的判断。当然，如果处理的是最近两年产生的测序数据，基本上都是Phred+33的，但从NCBI SRA数据库下载的较早的数据可能不同，需要注意。
1.2 什么序列适合用FASTA保存，什么序列适合FASTQ保存？
单纯的蛋白或者核酸的序列信息一般用FASTA格式保存，而测序文件一般用包含仪器信息和测序质量的FASTQ格式保存。
网友评论