RNA-seq测序质控

最后发布时间 : 2023-10-30 00:11:17 浏览量 :

质控的原因及相关软件

在A survey of best practices for RNA-seq data analysis里面，提到了人类基因组应该有70%~90%的比对率，并且多比对read（multi-mapping reads）数量要少。另外比对在外显子和所比对链（uniformity of read coverage on exons and the mapped strand）的覆盖度要保持一致。因此，可以对之前得到的BAM比对文件进行质检。
对BAM文件进行QC的软件包括：
Qualimap：对二代数据进行质控的综合软件
Picard：综合质控学习软件。
RSeQC是发表于2012年的一个RNA-Seq质控工具，属于python包。提供了一系列有用的小工具能够评估高通量测序尤其是RNA-seq数据。比如一些基本模块：检查序列质量、核酸组分偏性、PCR偏性、GC含量偏性，还有RNA-seq特异性模块：评估测序饱和度、映射读数分布、覆盖均匀性、链特异性、转录水平RNA完整性等。

转录组测序数据饱和度检验

DNA序列：测序不饱和，影响拼接
RNA-seq：定量不准
宏基因组：不能反映物种组成
该分析类似于群落研究中的稀疏曲线，用于评估测序深度是否足够。

该图中横坐标表示有效比对reads的百分比，纵轴表示该取样条件下表达量与终值的偏差比例，数值越趋近于1则表示表达量越趋于饱和。

测序数据量对于NGS分析非常重要，测序数据量过低，不能有效覆盖基因组完整信息，测序数据量过高，会造成冗余，不够经济。为了验证当前测序量能否满足需求，或者说加大测序量是否能够进一步挖掘信息，通常需要进
行饱和度分析。
充足的有效数据是信息分析准确的必要条件。转录组测序检测到的基因数目与测序数据量成正相关性，即测序数据量越大，检测到的基因数目越多。但一个物种的基因数目是有限的，而且基因转录具有时间特异性和空间特异性，所以随着测序量的增加，检测到的基因数目会趋于饱和。
为了评估数据是否充足，需要查看随着测序数据量的增加，新检测到的基因是否越来越少或没有，即检测到的基因数目是否趋于饱和。
使用各样品的Mapped Reads对检测到的基因数目的饱和情况进行模拟，绘制曲线图如下：

注：通过将Mapped Reads等量地分成100份，逐渐增加数据查看检测到的基因数量来绘制饱和度曲线。横坐标为Reads数目（以10^6为单位），纵坐标为检测到的基因数量（以10^3为单位）。

如果上图中各曲线越往右延伸，斜率越小，逐渐趋于平坦，说明各样品随着测序数据量的增加，新检测到的基因越来越少，趋于饱和，有效测序数据量充足。

图中不同颜色线条代表该样品中不同表达水平的基因表达量饱和度曲线，越早到达平台期则表示该样品中该基因表达量越高；反之则表示该基因表达量较低。

如图中多数基因的表达饱和度曲线最终并未达到1，则表明该样品的测序量不足以代表真实的基因表达水品，应加大对该样品的测序深度。

表达量饱和度

使用 RSeQC 分析表达量饱和度，即：对测序结果按 5%，10%，15%...100% 的比例进行抽样，对不同的抽样比例分别计算所有基因表达量（即每个基因计算 20次），然后与实际表达量（假设 100% 抽样情况下得出的为实际表达量）进行比较，获得相对误差。饱和度分析主要目的是评估所测数据量是否足够用于正确计算基因表达量。理论上在较少数据量情况下基因的表达量与实际表达量偏差较大，而当数据量达到饱和阈值后，数据量再增长，基因表达量也近乎不变，此时基因表达量不再受数据量的影响。

注：横坐标是重采样的比率，纵坐标为该比率下基因的表达量与实际表达量的相对误差。基因按表达量水平被划入四个组中，Q1：表达量水平居于倒数25%的基因；Q2：表达量水平居于倒数25%-50%的基因；Q3：表达量水平居于前25%-50%的基因；Q4：表达量水平居于前25%的基因。

rRNA污染率评估

由于用于转录组建库的Total RNA中大部分为rRNA，因此在建库过程中，很难完全去除rRNA，有其是原核生物转录组测序。

可以通过分析测序结果中rRNA的比例，以评估转录组文库构建是否合格，通常的做法是，将指控后的高质量reads与参考基因组进行比对，并更具比对结果的注释信息统计测序数据中属于rRNA的比例，一般认为rRNA含量低于15%为正常。

冗余序列分析（重复序列统计）

冗余序列主要来自于建库过程中的PCR扩增富集过程，由于PCR可能具有偏好性，因此可能导致部分基因扩增的不均匀，从而使测序结果中该基因序列的冗余性增高。

显示PCR重复序列的分布，一种是定义序列一样为重复序列，一种是定位map到同一位置的为重复序列

正常的结果为整体趋势呈线性平滑延伸扩散，随着横轴数值的增高，对应纵轴的数值下降。
如出现峰值，则说明冗余序列含量异常。

不匹配统计

mismatch_profile.py -l 125 -i Col-16_1_unique.bam -o Col-16_1_mismatch_profile

显示不匹配位点在reads位置的统计

GC含量统计

read_GC.py -i Col-16_1_unique.bam -o Col-16_1_read_GC

计算插入片段大小

inner_distance.py -i Col-16_2_unique.bam -o Col-16_2_inner_distance -r /opt/database/Arabidopsis/TAIR10.bed

得到插入片段大小的平均值mean与标准偏差SD。

基因覆盖均一度

基因覆盖度分析是样品中所有基因的5’到 3’区域上序列覆盖情况的综合展现，用于评估测序实验结果的均一性 (或是偏向性)。

geneBody_coverage.py \
    -r  genome.bed  \
    -i hisat2.bam \
    -o output

如果在靠近左端有明显的峰值，说明测序结果有明显的5’偏向性。如果在靠近右端有明显的峰值，说明测序结果有明显的3’偏向性。

若测序结果有明显的偏向性，则表明测序对于目的样品基因的覆盖率不均一，这可能是由于样本本身的序列结构导致PCR富集的偏好，也有可能是用于建库的RNA质量不高。

剪切饱和度分析

unction_saturation.py -i Col-16_2_unique.bam -r /opt/database/Arabidopsis/TAIR10.bed -o Col-16_2_junction_saturation

判断当前测序深度是否足够用来执行选择性可变剪切分析

剪切位点统计

junction_annotation.py -i Col-16_2_unique.bam -r /opt/database/Arabidopsis/TAIR10.bed -o Col-16_2_junction_annotation

测序数据质控

去除adapter
去除低质量reads
去除reads部分低质量的区域

Fastqc

fastqc -t 20  \
../data/test_R*.gz  \
-o .

cutadapt

User guide

cutadapt -j 6 \ # 线程数
  --times 1  \ # 去除一侧adapter
  -e 0.1 \ #adpapter与reads 10%的错误率
  -O 3  \ # reads至少3bp的adapter去除
  --quality-cutoff 25  \ #去除3’低质量的reads
  -m 55 \ # 去除adapter后小于55bp的reads删除，理论上reads不能小于40
  -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \ #reads1的adapter
  -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT \ # reads2的adapter
  -o test_R1_cutadapt.fq.gz \# trim后的reads1
  -p test_R2_cutadapt.fq.gz \# trim后的reads2
  data/test_R1.fq.gz \ # 
  data/test_R2.fq.gz

图片alt

=== Summary ===

Total read pairs processed:          2,000,000
  Read 1 with adapter:                 156,263 (7.8%)
  Read 2 with adapter:                 168,535 (8.4%)

== Read fate breakdown ==
Pairs that were too short:              77,983 (3.9%)
Pairs written (passing filters):     1,922,017 (96.1%)

Total basepairs processed:   400,000,000 bp
  Read 1:   200,000,000 bp
  Read 2:   200,000,000 bp
Quality-trimmed:               8,226,684 bp (2.1%)
  Read 1:     1,908,358 bp
  Read 2:     6,318,326 bp
Total written (filtered):    379,501,887 bp (94.9%)
  Read 1:   189,947,981 bp
  Read 2:   189,553,906 bp

比对区域分布统计

将比对到基因组上的Reads分布情况进行统计，定位区域分为CDS（编码区）、Intron（内含子）、Intergenic（基因间区）和UTR（ 5'和3'非翻译区）。在基因组注释较为完全的物种中，通常比对到CDS（编码区）的reads含量最高比对到Intron（内含子）区域的reads来源于pre-mRNA的残留或由可变剪切过程中发生的内含子保留事件导致的，而比对到Intergenic（基因间区）的Reads则可能转录自新基因或新的非编码RNA。

比对到Exon的reads比例最高，因为转录组测序就是针对外显子转录所得的mRNA。
比对到Intron区域的reads可能来源于pre-mRNA的残留或可变剪接过程中发生的内含子滞留事件。
而比对到Intergenic的reads可能是由于参考基因组注释不完善。

该分析统计比对到基因组上reads的分布情况，定位区域分为Exon(外显子)、Intron(内含子)和Intergenic(基因间区)。
如果结果显示比对到外显子区域的reads比例不高，则可能是由于建库过程中mRNA以外的其它RNA残留较多，也可能是参考基因组本身错误较多。

reads在染色体上分布

http://bioinfo.online/html/article/20231953789.html

转录组简介下游质控