3、RNAseq(3)--对RNAseq测序数据的质量控制（fastqc）

Posted 2023-04-19

参考技术A 使用命令 fastqc -o <output dir> <seqfile1,seqfile2..> 来进行质量报告。

需要注意的是./fastqc前面的"."不可以省略

每个fastqc文件会获得一个质量分析报告，来描述此次RNA-seq的测序质量。 获取质量报告如图：

从read水平来总览，判断测序质量。 Encoding ：测序平台的版本，因为不同版本的 error p的计算方法不一样。 Total sequence：测序深度。一共测序的read数。是质量分析的主要参数。 Sequence length：测序长度。 %GC：GC碱基含量比，一般是物种特异性，比如人类是42%左右。

横坐标： 第1-100个测序得到的碱基 纵坐标： 测序质量评估。这里的Q=-10*lg10(error P),即20%代表1%的错误读取率，30%代表0.1%的错误读取率 箱型图： 红色线，是某个顺序下测序碱基所有测序质量的中位数。黄色块，是测序质量在25%-75%区域。蓝色线，平均数。 一般要求： 测序箱型图10%的线大于Q=20。Q20过滤法。

横坐标：同上。 纵坐标：tail的index编号。 目的：防止测序过程中某些tail受不可控因素测序质量低。 标准：蓝色表示质量高，浅色或暖色表示质量低，后续的分析可以去除低质量tail。

从read的总体测序质量分布来判定此次的测序质量，是质量分析的重要标准之一。 横坐标：表示read的测序质量Q=-10*lg10(error P)。 纵坐标：表示在该Q值下的read 的数量 标准：需要集中在高分区

横坐标：1-100的测序碱基位置 纵坐标：碱基百分比 标准：理论上，ATCG碱基的分布应该差别不大，即四条线应该大致平行状态。如果AT或CG差异超过10%，此项检测是危险的。一般是测序机器前几个碱基测序时候因为状态调整导致测序略有偏差，如果前几个碱基偏差较大，可以在后期将前几个碱基切掉。造成这个偏差较大的原因重要是由于测序数据中的adapter没有clean干净。所以在接下来的数据处理中需要对测序数据进行clean data的处理，下一个章节会详细介绍clean data的步骤。

统计read的碱基长度，本例理论上测序应该全是150bp。 横坐标：是read的碱基长度 纵坐标：是该长度下的read数量

横坐标：每个read的平局GC含量占比 纵坐标：一定GC比下的read数 标准：蓝色是理论值，红色是真实值。两者接近是比较好的状态。如果有双峰，可能混有了其他物种的DNA序列。

一般测序在初步生成fastq文件时候，adapter会被去除，但是有的会没有去除或者遗漏部分adapter。所以这一步是检测RNA-seq测序过程中adapter是否去除。如果没有去除会严重影响后续的比对工作。没有去除的adapter在质量处理环节会被处理掉。

参考链接：
https://zhuanlan.zhihu.com/p/61847802

RNA-seq原理

参考技术A 测序技术发展：

1977Sanger测序--1996焦磷酸测序--2003cmPCR--2003ZMW---2012纳米孔测序

RNA-seq的一些技术限制，测序误差主要由生物学误差(生物学重复，比如取30只小鼠采样)和技术性误差(技术性重复，比如对1只小鼠采样3次)造成，如果想要得到的数据为无偏的，那么生物学重复最重要，因为生物个体代表着样本，而技术手段只会造成不可控干扰。总的来说，只做技术性重复的实验结果偏差最大，技术性重复+生物学重复的实验结果偏差也可能较大，除非生物学重复远大于技术性重复(因为当生物学重复次数不足时，技术性重复能扩大样本单一的影响)，无论如何，多做生物学重复，这有助于你的结论被其他人复现。

原理详解：

A 为了保证细胞在标记的过程中是单独分开的，10X开发了微流体设备(microfuidic device)进行预处理，设备有三个上样孔，分别加入你的1.样本细胞悬液(Sample) 2.凝胶小球(Beads) 3.分离液(Oil)，下图为具体设备的示意图。

当我们把样本细胞悬液加入设备时，每一个细胞会与凝胶小球单独结合，然后被分离液包裹，形成一个油包水的密闭小液滴(droplet)。进一步地，细胞和凝胶小球相遇不久后会裂解，释放出里面的各种物质，RNA(mRNA、tRNA、rRNA)，蛋白质，脂质，DNA等。实际上Beads上联接了不同的接头，其中有一个接头包含ploy(dT)序列，在细胞裂解后释放的核酸中，只有mRNA带有polyA tail，于是Beads的poly(dT)接头就可以从众多的裂解产物里捕获到mRNA(实际上drop-seq采用3'端测序，就是为了检测polyA tail)。

Master Mix中带有反转录试剂，当mRNA被捕获后，就可以从它的3‘端开始作为模板，进行反转录出cDNA的第一条链，这第一条链就沿着poly(dT)序列延申，长在了beads上，形成了图一7中的STAMPs，接着我们把反转录出来的cDNA序列洗脱，以cDNA的第一条链为模板，进行PCR，合成cDNA的第二条链，然后就是我们熟悉的cDNA扩增以及illumina测序。

如何确定测序序列来自哪个细胞？single cell的RNA-seq和bulk的RNA-seq的最大区别是什么？是barcode，或者说是cell barcode(实际上DNA自带barcode，cell barcode是人为控制的)。每一种single cell的beads上都有着相同的cell barcode(beads与beads间的cell barcode是不同的)，假设每个beads只捕获一个cell，那么则每个cell都被cell barcode 单独标记了。

如何保证每个beads只捕获一个cell？第一是控制cell和beads的流速，第二是beads的数目远远超过cell的数目，即绝大多数的beads都是空的，只有少数的才捕获到了cell。但是还是有个别的droplet里面会两个或者更多的细胞，这就需要进行质控(QualityControl)。

接下来可以参照10X Genomics的说明书详解single cell RNA-seq的barcode。

实际上beads上一开始只接了Read1、Barcode、Poly(dT)。

名词解释：

Poly(dT): 用来和mRNA的polyA结合，捕获mRNA

UMI: 用来标记不同的PCR产物(用于count计数)。为了减少由于复制引起的误差(重复抽样导致重复计数)，人们在一些单细胞测序的步骤中增加了UMI(unique molecular identifiers)，UMIs 是由 4-10 个随机核苷酸组成的序列，在 mRNA 反转录后，进入到文库中，每一个 mRNA，随机连上一个 UMI，因此可以计数不同的 UMI，最终计数 mRNA 的数量。

10X Barcode: 用来标记不同的single cell

Sample Index: 用来标记不同的sample

P5和P7: 用来进行illumina的桥式PCR测序

Truseq Read 1、2: 用来进行连接beads，cDNA的PCR扩增和加P7接头

在这些序列中，P5、P7、Truseq Read 1、2 的序列是已知的。

其他的序列是怎么一步一步添加上去的？

具体步骤：

利用Poly(dT)来捕获mRNA，在mRNA的5'端插入TSO(Template Switch Oligo模板切换低聚糖)引物，然后从mRNA的polyA开始反转录，直至mRNA的DNA序列被转录完成，然后在beads序列的3'端插入CCC，再对mRNA的TSO进行反转录，至此完成了cDNA的第一条链(序列顺序和mRNA逆序)。上述步骤很重要，因为中间cDNA的序列我们是不知道的(仪器测序长度有限)，如果不加上这个接头，就没有办法设计引物来合成cDNA的第二条链。

将mRNA溶解，对cDNA的第一条链加入UMI引物，以cDNA的第一条链为模板合成cDNA的第二条链。最后使用PCR(聚合酶链式反应)对cDNA(拷贝DNA)进行扩增(为了富集)。

PCR原理

因为II代测序(NGS)的illumina测序不能测很长的seq，约为200-700bp，所以不能测得mRNA全长，因此需要进一步把合成的cDNA利用酶打断到illumina能测的长度(长度有些随机，比如300bp的cDNA能通过头尾150bp完整测序，但700bp的cDNA只能通过头尾150bp测序+参考基因组推断出来)。然后在cDNA的3'端插入Truseq Read2引物(和Truseq Read1引物匹配为头尾，中间序列就是reads)、P5、P7。

最后的测序数据(reads)从Truseq Read1后的10X Barcode开始，一直到Truseq Read2为止。

PCR扩增是对cDNA单链进行复制，后面的桥式PCR是对完整的样本进行复制(增加数据深度)，总的来说各个cDNA呈均匀分布，然后进行抽样。

RNA-seq duplications有PCR duplication(最主要)、cluster duplication、optical duplication。

实际上仪器会对核苷酸进行染色，然后判断颜色确定ATCG碱基，因此有很多原因会导致机器误判，和后续QC有关。

1.某些核苷酸对颜色附着不明显

2.大片区域颜色相同(相同类型核苷酸)，而其中仅有几个颜色不同的点(不同类型的核苷酸)