2019-04-15-生物信息学100个基础问题 - 丁立的博客

参考资料：孟浩巍博士的知乎：https://zhuanlan.zhihu.com/ngs-learning

第一题 FASTQ与FASTA

1.1 掌握FASTQ格式

1.1.1 格式有什么特点？

fastq内容格式有4行：

第1行主要储存序列测序时的坐标等信息；

举个例子：

@ST-E00126:128:HJFLHCCXX:2:1101:7405:1133   
@，开始的标记符号;
ST-E00126:128:HJFLHCCXX，测序仪唯一的设备名称;  
2，lane的编号;           
1101，tail的坐标;
7405，在tail中的X坐标;
1133，在tail中的Y坐标

第2行就是测序得到的序列信息，一般用ATCGN表示，其中N用于荧光信号干扰无法判断到底是哪个碱基时的代表符号；
第3行以“+”开始，可以储存一些附加信息，但目前的测序fastq文件这一行一般是空的；
第4行储存的是质量信息，与第二行的碱基序列是一一对应的，其中的每一个符号对应的ASCII值是经过换算的phred值，可以简单理解为对应位置碱基的测序质量值，越大说明测序的质量越好。

1.1.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”

1.1.3 phred33与phred64是什么意思？

质量字符的ASCII值和质量得分的关系有如下两种：可以粗略分为Phred+33和Phred+64，这里的33和64就是指ASCII值转换为Q该减去的数值。

在处理测序数据时，因为一些软件会根据碱基质量得分的不同做不同的处理，常要指定正确的编码方式，有必要对质量字符与质量得分的关系（Phred+33或Phred+64）作出正确的判断。当然，如果处理的是最近两年产生的测序数据，基本上都是Phred+33的，但从NCBI SRA数据库下载的较早的数据可能不同，需要注意。

1.2 FASTA格式的构成是怎样的，有什么样的规律？

fasta格式用于储存序列，可以储存DNA、RNA和蛋白质序列，一般分为两个部分，第1行是以>开头的序列描述信息，包括数据库中的编号，序列名称，序列类型，剩余的为序列信息,以蛋白质和mRNA序列文件为例:

蛋白质fasta文件

以>开头
sp|P69905 数据库编码
HBA_HUMAN Hemoglobin subunit alpha 蛋白质名称
OS=Homo sapiens 所属物种
GN=HBA1 基因名称
>sp|P69905|HBA_HUMAN Hemoglobin subunit alpha OS=Homo sapiens GN=HBA1
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSHGSAQVKGHGKK
VADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLLSHCLLVTLAAHLPAEFTPAVHASL
DKLASVSTVLTSKYR

核酸序列文件（mRNA序列中的U均用T来代替）

以>开头
gi|13650073 基因ID
gb|AF349571.1 genebank编号
Homo sapiens hemoglobin alpha-1 globin chain (HBA1) 基因名称
mRNA, complete cds 序列类型

>gi|13650073|gb|AF349571.1| Homo sapiens hemoglobin alpha-1 globin chain (HBA1) 
mRNA, complete cds 
CCCACAGACTCAGAGAGAACCCACCATGGTGCTGTCTCCTGACGACAAGACCAACGTCAAGGCCGCCTGGGGTAAG
GTCGGCGCGCACGCTGGCGAGTATGGTGCGGAGGCCCTGGAGAGGATGTTCCTGTCCTTCCCCACCACCAAGACCT
ACTTCCCGCACTTCGACCTGAGCCACGGCTCTGCCCAGGTTAAGGGCCACGGCAAGAAGGTGGCCGACGCGCTGAC
CAACGCCGTGGCGCACGTGGACGACATGCCCAACGCGCTGTCCGCCCTGAGCGACCTGCACGCGCACAAGCTTCGG
GTGGACCCGGTCAACTTCAAGCTCCTAAGCCACTGCCTGCTGGTGACCCTGGCCGCCCACCTCCCCGCCGAGTTCA
CCCCTGCGGTGCACGCCTCCCTGGACAAGTTCCTGGCTTCTGTGAGCACCGTGCTGACCTCCAAATACCGTTAAGC
TGGAGCCTCGGTGGCCATGCTTCTTGCCCCTTTG

1.3 什么序列适合用FASTA保存，什么序列适合用FASTQ保存？

单纯的蛋白或者核酸的序列信息一般用FASTA格式保存，
而测序文件一般用包含仪器信息和测序质量的FASTQ格式保存。

第二题测序技术初探

现在实验室或者公司常用第1代测序与第2代测序，那么：

2.1 第1代测序 sanger 测序法的原理是什么？通量比较低的核心原因是什么？

sanger法测序及双脱氧链终止法，它采取DNA复制原理，通过在DNA复制过程中添加双脱氧三磷酸核苷酸（ddNTP）终止DNA链的延伸，在DNA链不同位置的延伸终止判断该位置的碱基类型。但是凝胶电泳的时间较长，导致sanger法测序通量低。

2.2 作为2006年正式发布的illumina测序技术，或者称为第2代测序技术的代表性技术，其最大的特点是什么？

高通量，成本低，但测序长度较短。

2.3 Illumina测序技术的核心是什么？

核心内容有两个：
一个是桥式PCR，主要用于扩大信号；
另一个是4色荧光可逆终止反应，使illumina测序可以实现边合成边测序的技术。

2.4 Illumina测序技术为什么不能像第1代测序技术一样测500bp以上？

主要的原因有两个：

一方面测序时，经过长时间的PCR，会有不同步的情况；

比如一开始1个cluster中是100个完全一样的DNA链，但是经过1轮增加碱基，其中99个都加入了1个碱基，显示了红色，另外1个没有加入碱基，不显示颜色。这时候整体为红色，我们可以顺利得到结果。随后，在第2轮再加入碱基进行合成的时候，之前没有加入的加入了1个碱基显示红色，剩下的99个显示绿色，这个时候就会出现杂信号。当测序长度不断延长，这个杂信号会越来越多，最后很有可能出现50个红，50个绿色，这时信号不足以判断碱基类型；

第二就是测序过程中合成酶的活性越来越不稳定，后面碱基添加出现问题。

第三题 Illumina测序技术细节探究

目前我们最常使用的就是Illumina公司的测序技术，Illumina公司的测序技术最明显的几个特点是：价格低，通量高，测序读长短。那么我们今天的问题，就是围绕Illumina测序技术的细节来提问的。

3.1 什么是Illumina测序adapter？同一批上机的adapter序列一样吗？它的作用是什么？

adapter的中文意思为适配器或者接口，在illumina测序过程中关键一步是将文库片段固定在flowcell上，然后通过桥式PCR将片段扩增，在被打断成300~500bp的长度的片段末端被补平后adaptor将被添加到片段两端，一方面用于将片段固定在flowcell上，同时adaptor中还包含桥式PCR所需要的引物。

3.2 一个完整的Illumina测序过程是那几步？

完整的测序过程仅包含两步：
第一是桥式PCR扩增，
第二是以4色荧光可逆终止反应为核心技术的测序；

3.3 什么是桥式PCR技术？为什么要进行桥式PCR？

加上adaptor之后的DNA样品与flowcell上固定的oligo（寡链核苷酸）匹配后就被固定在flowcell上，通过桥式PCR进行扩增成cluster,便于后面的荧光测序，主要步骤为：

进行第一轮扩增，将序列补成双链。
加入NaOH强碱性溶液破坏DNA的双链，并洗脱。由于最开始的序列是使用化学键连接的，所以不会被洗。
加入缓冲溶液，这时候序列自由端的部分就会和旁边的oligo进行匹配。
进行一轮PCR，在PCR的过程中，序列是弯成桥状，所以叫桥式PCR，一轮桥式PCR可以使得序列扩增1倍
如此循环下去，就会得到一个具有完全相同序列的cluster

3.4 我们都说，测序结果会包含index，那么index是什么？有什么作用？

一条lane能测得的数据量在30G左右，而一个样品的测序量一般不会这么大，所以在建库的时候对每一种样品的接头加上不同的标签序列，这个标签就叫做Index，有了index就可以同时在一个lane中测多种数据了，后期可以根据index将数据分开；

3.5 我们所说的flowcell，lane，tile都是什么意思？

flowcell 是指Illumina测序时，测序反应发生的位置，1个flowcell含有8条lane
lane 每一个flowcell上都有8条泳道，用于测序反应，可以添加试剂，洗脱等等
tile 每一次测序荧光扫描的最小单位

3.6 Illumina测序结果质量表示方法采用的是Phred33还是Phred64？

最新的测序质量结果一般都为Phred33，但是早期的测序数据可能出现Phred64。

第四题 Illumina测序技术细节探究 II

4.1 Illumina目前主流的测序仪都有哪几种型号？各自大概的通量是多少？（也就是1个run能跑出多少数据）

目前主流的测序仪及其通量主要是：

Hiseq2500（50-1000Gb）
Hiseq3000（125-750Gb）
Hiseq4000（125-1500Gb）
Hiseq X Five（900-1800Gb）
Hiseq X Ten（900-1800Gb）。

4.2 Illumina目前的测序技术，最核心的就是边合成边测序，即我们常说的 Sequencing by synthesis （SBS），那么为什么能够实现SBS？

经过桥式PCR之后同一段序列已经成簇，下一段就是开始进行测序，这一步比较简单，就是加入primer，然后添加经过特殊处理的ATCG四种碱基。

特殊的地方有两点：

一个是碱基部分加入了荧光基团，可以激发出不同的颜色，
另一个是脱氧核糖3号位加入了叠氮基团而不是常规的羟基，这个叠氮集团保证了每次只能够在序列上添加1个碱基。

这样每1轮测序，保证只有1个碱基加入的当前测序链。这时候测序仪会发出激发光，并扫描荧光。因为一个cluster中所有的序列是一样的，所以理论上，这时候cluster中发出的荧光应该颜色一致。随后加入试剂，将脱氧核糖3号位的—N2改变成—OH，然后切掉部分荧光基团，使其在下一轮反应中，不再发出荧光。如此往复，就可以测出序列的内容。

4.3 “Illumina测序技术为什么不能像第1代测序技术一样测500bp以上？”，这里面主要涉及到两种错误，一种叫phasing，一种叫pre-phasing，分别是什么意思？

通俗来讲，phasing表示本来同步添加的碱基有一些没加上，而pre-phasing则是加多了，都会导致当前bp的荧光检测出现噪音，造成phasing的主要原因是合成酶的活性降低，而pre-phasing则可能是叠氮基团性质不稳定，转化为羟基在一步检测中添加了不止一个碱基。

第五题测序建库的adapter

5.1 adapter是什么意思？adapter与primer有什么区别？

adapter在中文是适配器或者接口的意思，在前面的内容中已经提到将测序序列打碎成片断后要将末端补平然后添加adapter，用于与flowcell上的oligo匹配固定并为后续桥式PCR做准备，而前面提到的Index与adapter之间的位置关系一般为adapter1-Index-fragment-adapter2，adapter2通过与oligo互补连接在flowcell上，在进行完桥式PCR之后进行测序时，添加primer，这一段primer的序列是与Index互补的而非adapter1，所以最终拿到的测序结果应该是Index+fragment+adapter2或者Index+部分fragment。

5.2 比如最终的测序结果是 AATTCCGGATCGATCG…，那么adapter的序列可能出现在哪一端，还是两端都有可能出现？为什么？

一般出现在3’端，在上面第1题中已经说到，最终的测序结果应该是Index+fragment+adapter2或者Index+部分fragment，也就是说测序的方向是从5’到3’，adapter只可能出现在3’端。

第六题读懂FastQC报告 Part I

6.1 图中的横坐标表示什么意思？

横轴是测序序列第1个碱基到第150个碱基。

6.2 图中的纵坐标表示什么意思？

纵坐标表示每一bp所对应的测序质量值。
测序质量值（Q）：Q=-10 * log10（error P），
即碱基判断错误概率值P取log10之后再乘以-10，
得到的结果再加上phred值对应ASCII表所得到的值就是该碱基测序的质量值。
即Q20表示1%的错误率，Q30表示0.1%的错误率。

6.3 图中的蓝色线是什么意思？

蓝色的细线是各个位置的质量值的平均值的连线。

6.4 图中的box下面的bar，上面的bar，箱体的下沿，箱体的上沿，箱体内部的横线分别代表什么意思？

每1个boxplot，都是该位置的所有序列的测序质量的一个统计；
上面的bar是90%分位数；
下面的bar是10%分位数；
箱子中间的横线是50%分位数；
箱体上沿是75%分位数；
箱体下沿是25%分位数。
分位数的一个简单解释是：如果一组数的25%分位数是a，意味着a超过了这组数中25%数字的大小。

6.5 图1与图2最主要的区别在哪里？结合我们之前的问题，为什么会出现这种情况？

相比于reads1的测序结果，图2中reads2测序质量均匀性差，准确率低。
主要原因是：
read2的测序是在reads150bp测序完成后，
forward strands再通过1次桥式PCR合成reverse strands；在这之后再进行荧光测序；
测序质量差的主要原因是因为长时间测序结束后，合成酶的活性降低，
导致合成时加不上一些碱基，最终同步性变差，主要属于phasing错误。

第七题读懂FastQC报告 Part II

boxplot一般是认为FastQC几张必看的质控图之一。一般情况下FastQC的结果会包含下面几个图，而我们主要会看下图圈出来的几个。

上一题讨论了“Per base sequence quality”，本题先来讨论 “Per base sequence content

7.1 图1与图2中横坐标是什么意思？纵坐标是什么意思？

横轴代表1到150bp；纵轴代表ATCG在该bp的百分比。

7.2 图1是1个正常的DNA 全基因组测序结果，为什么前面的几bp线是波动的？后面的线是平衡的？

根据wason-crick配对原则，A和r应该相等，G和C应该相等；
但是一般测序的时候，刚开始测序仪状态不稳定，很可能出现不平衡的情况。
像这种情况，如果测序的得分很高，可以不进行trim开始部分的序列信息；
如果测序得分很低，需要进行trim开始部分的序列信息。

7.3 图2是1个特殊RNA建库的测序结果，4条线出现波动更可能是什么原因造成的？

RNA是单链核苷酸，GC或AT的含量并没有直接的关系，
4条线的波动是由于模板中ATCG的量不同导致的（就是所谓Gc不平衡）。

7.4 在图1中你能不能看出一个恒定的量？（提示，同一物种间相同，不同物种间一般不同）如果能看出来，这个量是什么？数值大约是多少？

GC含量在同一物种中是一个恒定值。图中Gc总体比例大约在42%（目测）。

第八题读懂FastQC报告 Part III

本题来研读2张图：

第1张图是：Per sequence GC content

图1-1 human 全基因组测序FastQC 结果图

图1-2 human 全基因组测序FastQC 结果图

第2张图是：Sequence Length Distribution

8.1 图1-1 与图1-2 中的横坐标是什么意思？纵坐标是什么意思？

横轴是0-100%；纵轴是拥有相对Gc含量的序列所对应的数量。

8.2 图1-1中是human全基因组测序，结合昨天的问题，那么peak的中间大约应该在多少？

上一次的问题中GC总体比例大约在42左右，
所以如果这是human全基因组测序，那么peak在横坐标对应42的位置比较好。

8.3 图1-2与图1-1有哪些显著的不同？造成这些不同的原因有可能是什么？遇到这个问题，我们通常应该做些什么？

图1-1有一个peak，并且与理论值（蓝线）基本重合；
而图1-2有两个，并且其中一个peak与蓝线相差很多，

当红色的线出现双峰，基本是混入了其他物种的DNA序列

遇到这个问题，首先进行mapping统计有多大比例reads map到了目标参考基因组上，
如果比例非常低说明污染严重，数据不可用；
如果大部分reads都map成功，剩余一部分可以通过blast检查是混入了哪些污染物，
过滤掉这些reads就可以，不影响后续的分析。

8.4 图2-1的横坐标是什么意思？纵坐标是什么意思？

横坐标代表序列长度，纵坐标代表长度为某一bp的序列所对应的数量。

8.5 图2-1是刚下机的fastq数据进行FastQC 结果图，有什么特点？为什么会出现这样的结果？如果对刚下机的fastq数据进行cutadapter，图2-1还会是这样的结果吗？为什么？

测序仪成功下机的数据都是整齐的一定长度的序列，比如最常用的i1lumina x Ten是双端150bp；
测序过程当中产生的不足150bp的序列在下机时已经被过滤掉了；
如果进行cut adapter，序列的长度将不一致，因为reada中包含信息的insert的长度并不完全一致，
150bp的测序长度是否已经包含了adapter的序列是未知的，因此cutadapter之后的reads长度不同。

8.6 能力扩展题：

请想办法，计算Human genome 19（hg19）每一条染色体的GC含量。

思路1
-UCSC或Ensembl下载Human genome 19（hg19）染色体序列；
-直接把下载的序列使用Pastqc检测序列“质量”，
-report中sequence content across all bases会显示GC结果。

思路2
-Ucsc或Emble下载Human genome 19（hg19）染色体序列；
-写python程序依次读取序列内容；
-统计每一条染色体的GC含量；

第九题读懂FastQC报告中的duplicate问题

本题来详细聊聊duplicate问题。

duplicate的产生主要是因为Illumina建库的过程中，一般会需要使用PCR来帮助扩增插入序列的浓度。在扩增的过程中，如果PCR扩增轮数过大，就会出现duplicate的问题，即产生一模一样的若干条序列。

FastQC中“Sequence Duplication Levels”图是用来刻画duplicate情况的。

9.1 图1中的横坐标是什么意思，纵坐标是什么意思？

横坐标代表序列重复水平；纵坐标代表重复水平序列占所有序列的百分比。

9.2 图1中的红线和蓝线分别代表什么意思？

红线代表去duplicate之后序列理论重复性分布（服从possion distribution或者
binomial distribution）情况，
蓝线代表全部的序列重复性分布情况。

9.3 图1中的duplicate是全部序列的duplicate的情况吗？还是随机筛选了一部分？为什么要这样做？

是选择的每一个文件里前100，000条序列作为样本进行的计算，因为样本本身很大，
前100000已经能够代表样本的重复性。

9.4 如果让你写程序，判断1个fastq文件中duplicate的比例，你的大概思路是什么？

#Python风格的伪代码：

#第1步对序列进行排序
sort the FAsTQ file by the sequence，and names as sorted_file；

#第2步对排序的序列统计是否为duplicate
total_num=1 
duplication_num=0 

reads_1=sorted_file.readline（）
for reads 2 in sorted_file 
  if reads 1 ==reads 2 
    duplication_num=duplication_num+1 
  else 
    reads1=reads_2 
  total_num=total_num+1 
print total_num 
print duplication_num

9.5 既然谈到了duplicate的问题，那就存在remove duplicate的问题，什么情况下应该去duplicate，什么情况下不去除？

DNA-Seq中序列如果是随机打断需要考虑deduplicaion；葡切的样本一般不需要考虑这个问题；
RNA-Seq一般不考虑remove duplication（有paper专门讨论过这个问题）；
单细胞测序需要建库过程中需要添加random barcode，且必须考虑duplication。

第十题读懂FastQC报告之adapter与kmer

在第5题的时候讨论过adapter的问题，我们知道adapter的最主要的作用是为了能够与flowcell连接，方便进行桥式PCR。那么我们的fastq文件中到底含不含adapter呢？FastQC报告就能告诉我们。同时呢，本题还会讨论kmer的问题，相关的报告FastQC也会输出出来。

Part I adapter部分 图 1-1 1个正常的adapter报告

图 1-2 1个RNA-Seq的adapter报告

10.1 关于adapter的问题：

10.1.1 Illumina的通用adapter序列是什么？图1-1与图1-2中的各种不同颜色的图例是什么意思？

Illumina Paired End Adapters（cannot be used for multiplexing）

Top adapter 
5' ACACTCTTTCCCTACACGACGCTCTTCCGATC*T 3'

Bottom adapter 
5' P-GATCGGAAGAGCGGTTCAGCAGGAATGCCGAG 3'

**来源**

http：//bioinformatics.cvr.ac.uk/blog/i1lumina-adapter-and-primer-sequences/

不同颜色的图例代表不同的测序通用的adapter；
如果在当时fastqc分析的时候-a选项没有内容，则默认使用图例中的通用adapter序列进行统计。

10.1.2 图1-1与图1-2中的横坐标与纵坐标分别是什么意思？

横坐标代表reads中的位置，纵坐标代表adapter序列含量的百分比。

10.1.3图1-1与图1-2中最显著的差异是什么？如果两者都是RNA-Seq的数据，哪个可以继续下游分析，哪个不能够进行下游分析？为什么？

图1-1的结果表明少量的序列3'端测到了少部分adapter序列，
且测序时使用的是红色图例代表的adapter；

图1-2的结果表明测序时使用的是红色图例代表的adapter；
除了前20bp，reads后半段序列有很高比例都是adapter，建库可能存在问题。

Part II kmer部分 图 2-1 正常的RNA-Seq建库kmer统计

图 2-2 加入random barcode的RNA-Seq建库kmer统计

图 2-3 kmer的统计显著性分析

10.2 关于kmer的问题：

10.2.1 kmer就是一定长度的序列，比如AATTCCGG就可以叫做8-mer。那么图2-1余图2-2中的横坐标什么意思？纵坐标什么意思？

横坐标代表短序列的长度，纵坐标代表某长度的短序列在所有reads中所出现的频率百分比。

10.2.2 图2-1与图2-2中哪个kmer问题比较严重？为什么？

图2-2的问题更严重，因为该图中kmer的出现位置集中且数量较多，
可能是加入了random barcode，出现了duplication问题。

10.2.3 图2-2中是在reads的5’端加入了约10bp左右的随机序列，结合生物信息学100个基础问题 —— 第9题读懂FastQC报告中的duplicate问题这样做的目的是什么？

一般在进行RNA-Seq测序时是不会进行remove duplication；
但是一些比较特殊的建库流程比如说单细胞RNA-Seq测序时PCR扩增轮数较多，可能出现大量的duplication；

对于这些建库方法，通常需要添加random barcode，
然后需要根据random barcode进行remove duplication。

10.3 思考题：

10.3.1 图2-3是FastQC生成的kmer是否显著的统计报告。其中的每一列是什么意思？这个统计显著性检验计算的p-value是使用什么方法计算的？

第一列：kmer内容
第二列：kmer在序列中某一位置出现的观测数量
第三列：二项分布统计检验P-value
第四列：kmer在某一位置的观察值与理论值的比值
第五列：观察值与理论值比值最高值出现的位置