可视化测序质量是进行质量控制的重要参照,使得我们的质量控制可以有的放矢,那么这一节集结了小编对于质量可视化的一些思考,希望可以帮助您呀~
fastq文件是sequence文件的重要格式之一,首先我们先下载一个fastq文件以便于以后操作
Fastq文件下载指南
NCBI数据下载随记-CSDN博客这是小编写过的另外一篇从SRA数据库里面下载数据的指南,这里我们选择直接从终端下载数据
那么首先得到一个Bioproject的项目id,比如说我这里已经准备去下载PRJNA257197此项目的数据
esearch -db sra -query PRJNA257197
首先可以进行查询获得对此项目最基本的认知了解,当然了可以更一步深入
esearch -db sra -query PRJNA257197 | efetch -format runinfo > runinfo.csv
我们可以把runinfo转变为xlxs格式便于查看,这里就不再操作了,我们直接截取Run项目id(SRR**)
cat runinfo.csv | cut -f 1 -d ',' | tail +2 | head -10
这里我们就先截取一部分的SRR(因为我提前查看过,所以知道SRR是从第二行开始的,所以直接用tail截断了,如果没有查看的话,可以使用grep SRR 平替,效果一样)
然后我们得到一个SRR列表,咱们就下载第一个SRR(因为小编电脑网速实在感人,所以sorry~)
(这一篇只是起到指导作用,所以应该无伤大雅,读者可以多下载一些无妨)
fastq-dump -X 10000 --split-files SRR1972917
#只下载前10000reads,但是分开下载pair-end数据
我们接着就将我们得到的其中之一的SRR1972917_1.fastq进行Quality Control的第一步——可视化测序质量,但是首先让我们来看一下我们下载的文件的格式和属性是什么
fastq文件属性分析
首先我们先用nano简单看一下我们下载的fastq文件,可以认为基本的单位如下
那么第一行和第三行我们可以暂且忽略(后面在可视化环节会有作用),第二行代表的是基因的序列,而我们需要重点关注的事第四行,表示点对点的质量参数,如果是Sanger / Illumina 1.8+ 格式,那么比如我们的质量值是20,那么其+33对应的ASCII字符即为表示质量的参数。那么通过ASCII值对应的规律,我们得到的质量值又代表了什么呢?我们将质量值除以10比如这里是2,然后取负,并把其作为10的幂,得到0.01,那么就说明这个碱基测序出错的概率是0.01,当然我们不需要自己算,在后面的可视化可以非常清晰地看到
然后其实每一行数据都是一个read,这是illumina测序方式(这里就不多介绍了,可以暂且理解为DNA片段,但是不仅仅如此),我们这里看到这里每一个read都是101bp的长度,其实这应该是对原始数据操作之后的结果,不然实际实验测序当中是很难做到这样的。
Visualize Sequencing Data Quality
然后我们使用简单的一个命令,可以对我们下载的数据进行可视化(使用FastQC程序,其实这个名字有点误导效果,因为实际上这个程序仅仅提供可视化,不会提供Quality Control服务)
fastqc SRR1972917_1.fastq
然后就会生成一个html文件,里面是对fastq文件的质量可视化
下面是对每一参数和图像的详细解释
Basic Statistics
那么这个数据量小编认为还是很多的,首先告知了我们这个序列的测序方式,这里是Sanger/ Illumina1.9,那么总共的read数量是10000条(因为我们就下载了前10000条数据),每一个序列(这里sequence和read可以认为等价的)的长度是101bp,那将sequence数量和长度相乘就可以到总的数据量1Mbp。%GC含量具有显著的生物学意义,可以协助判断数据质量
Sequence Length Distribution
虽然显示的顺序不是这样的,但是小编建议可以先看到这个数据(对于后面的分析还是很重要的)
横轴表示的是序列的长度,纵轴表示的是在该序列长度下序列的数量
在前面已经铺垫过了,我们下载的这个数据每一个Sequence的都是101bp,所以只有一个峰,但是如果我们遇到了不同长度的序列的时候,这幅图就会提示我们后面的图像其他的一下信息(后面有介绍原因)
Per Base Sequence Quality
首先看横轴,横轴这里是101bp(有没有想到代表着什么?),表示每个Sequence的对应位点,纵轴是其质量值,表示每个序列在这个位点的质量值分布。这个图的信息量很多,咱们逐一分析
首先每个Sequence的对应位点是什么意思,我们刚才提到过我们的序列都是101bp的,所以这个1其实代表了这10000个序列的第一个碱基位,然后纵轴的箱线图就表示了这10000个碱基的质量值分布,黄色代表了下四分位数到上四分位数的区间,如果黄色箱子越矮就说明数据分布就越集中,而黑色线代表的是变异值,这是根据上四分位数和下四分位数计算出来的,如果不在这个黑色线的范围内,就说明这个数据就是变异值,换句话说黑色线包裹的范围就是可以接受的碱基质量,如果这个黑色线全部在绿色范围就说明这个碱基位的质量很高。
红色线表示中位数(其实蓝色也是表示中位数,只是由于渲染的原因使得存在偏差,可以通过蓝色线看走势,红色线看值)
是时候呼应伏笔了,那么如果我们获得的数据不是一致的长度怎么办?比如说有的序列是50bp,有的序列是100bp,那么在0-50bp的范围内,这两个序列的碱基都会纳入考量,但是在50-100bp的范围内,50bp就没有纳入考量了。所以我们先关注序列长度分布是非常重要的,这样我们才可以知道不同的序列信息究竟在哪个范围体现。
(但是秉持着深入理解的学习理念,大家有没有想过为什么要采取这样的碱基质量‘登记方法’,为什么每一Sequence的碱基位要对齐呢?这是因为测序的原理,每一个read的测序是同时的(或者几乎同时,会受酶活性,序列自身性质等等影响),那么可以认为在测对应碱基位受到外部因素的影响几乎是同步且接近的,那么我们采取这样的记录方式,就可以清楚地知道哪一个碱基位出现了较大误差,这样对于误差的统计才会有意义)
Per Sequence Quality Scores
对这幅图的理解很简单,横轴是序列的平均的质量值,纵轴是相应质量值的序列个数,那么如果在越高质量值的序列数量越多(出现峰值是其中标志),就说明这个序列的质量一般是越高的
Per Base Sequence Content
这个表示在每一个碱基位的ATCG各个碱基的比例,通过是否存在某一个碱基比例突然变高,可以反应测序是否具有碱基倾向性
Per Sequence GC Content
(红线)表示在特定的GC比例对应序列的个数,从理论上应该是单峰分布(蓝线),如果对于蓝线偏离程度越大,那么就说明数据质量一般越差(特别是多峰情况,可能是不同物种的DNA杂合在一起测序了)
Per Base N Content
N表示不确定的碱基,这里表示每个位点不确定碱基的相对比例。这一幅图是很友好了(每一个位点几乎比例均为0),但是说实话,如果质量很低和不确定也没有很大差别了
Sequence Duplication Levels
这一幅图的横轴是序列的重复次数,纵轴是特定重复次数对应的序列比例,一般来说理想情况下一个序列只有一次重复(其本身)所以集中在1处最好,但是如果存在PCR偏好或者特定mRNA高度表达,那么就会出现特定序列重复多次,在后面的level出现峰值,所以最上面的71.28%代表的就是没有出现重复的序列,或者说是有效序列
Overrepresented Sequence
这是全搜索所有序列(这里取的Sequence长度是50bp,注意这里的Sequence只是一段序列,不是read)得到的结果,表示这些序列在所有序列当中重复的次数和比例,这个数据通常可以反应接头污染(见后)
Adaptor Content
这个图显示的就是出现接头污染的比例,不同的颜色表示不同的接头(这会因为不同的测序手法和平台而有不同)
总结
我们对上述数据的分析,我们可以大致得到这个数据的质量如何,并且可以初步指示出现了哪些状况导致了这些测序数据上的质量下降,那么这也为我们进行Quality Control奠定了很好的基础。
小编水平有限,如有不足错误之处请多指教友善交流~