在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比对文件进行质检。 DNA序列:测序不饱和,影响拼接 RNA-seq:定量不准 宏基因组:不能反映物种组成 该图中横坐标表示有效比对reads的百分比,纵轴表示该取样条件下表达量与终值的偏差比例,数值越趋近于1则表示表达量越趋于饱和。 测序数据量对于NGS分析非常重要,测序数据量过低,不能有效覆盖基因组完整信息,测序数据量过高,会造成冗余,不够经济。为了验证当前测序量能否满足需求,或者说加大测序量是否能够进一步挖掘信息,通常需要进 注:通过将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到同一位置的为重复序列 正常的结果为整体趋势呈线性平滑延伸扩散,随着横轴数值的增高,对应纵轴的数值下降。 显示不匹配位点在reads位置的统计 GC含量统计 read_GC.py -i Col-16_1_unique.bam -o Col-16_1_read_GC得到插入片段大小的平均值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 . cutadaptUser 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 图片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(基因间区)。 |