Taro

  • 首页

  • 关于

  • 标签

  • 分类

  • 归档

  • 搜索

WGS分析笔记(1)- 原始数据以及质控

发表于 2019-06-18 | 更新于 2019-06-22 | 分类于 我的WGS
  • 二代测序方式分为三种:

    • single read 单端测序
    • paired-end read 双端测序
    • mate-pair read 配对测序
  • 虽然有三种方式,但是大多数数据(包括本课题)为双端测序,至于三者之间的差别、优缺点以及适用场合可以自行搜索。

  • 接下来讲的内容都是基于双端测序的。

  • 测序仪原始下机的数据我们称为raw data,二代测序是将DNA片段打断了再测的,每个测序片段我们称为read,质控完的数据称为clean data。既然是双端测序,那么文件就是成对出现的,分别记录reads两端的信息:一般的命名是*.fq1.gz、*.fq2.gz(’ * ‘ 表示通配符),这是一个fastq文件,通常以fq或fastq作为后缀,具体内容下方会介绍。这里我给大家展示我们课题的一个真实数据,给大家看看大小以及命名方式。

    原始数据展示

  • 我们可以看到,首先我把我的用户名打码了,呵呵,开个玩笑。

  • 首先这个样本总共有五个文件,这就是公司给的原始数据,我原封不动的上传到了我的服务器上,也没改名也没怎样,我在图片上画了五个蓝色的框,下面是分别代表的意思:

1
2
3
4
5
W2018001:是样本的名字,其实我一开始提交给公司的名字是2018001,为什么多了个W我也不知道啊。
NZTD...:这一串是公司自动生成的编号,他们内部使用的,不知道也不需要知道啥意思,不同的公司,可能没有这个
HCV...:flowcell_ID的信息,一般情况下一个样本是一样的
L1:flowcell_lane的信息,图中有L1和L2,这也是为什么一个数据他给我四个文件的原因,他在不同的lane上测的,这里还要注意的一点就是,这个lane的值可能是同一个,就是即使碰到都是L1也是可能的
1,2:分别代表reads的两端
  • 我们可以看到,这个30X的WGS的数据量还是很大的,所以为什么他会是压缩文件,不同公司给的命名可能不一样,不过大体不会相差太多,具体命名的含义也可以问问公司的。
  • 这里给到的是一个样本2对文件,实际情况是,可能有多对或者只有一对文件,对于这样的情况,我们的处理方式一般是看作一个文件处理,就是merge一下,有见过有些课题组是在这一步直接merge,我的处理方式是后期bam文件再merge。
  • 好了,拿到这么一个数据,我们做的第一件事情应该是校验数据的完整性!!!这点很重要,即使一般情况下,传输都不会出错,但这一步还是要做的,怎么去做呢,就看这个截图里的第一个文件MD5.txt,其实这就是一个文本文件,打开了看看是这样的

MD5信息

  • 就是这样的四行信息,也就是一个文件名对应一个校验码,最简单的校验方式,linux系统自带的MD5校验命令:
1
$ md5sum -c MD5.txt
  • 这个还是需要一点点时间的,大家要耐心等待。展示一下我这里的结果吧(其实在上传完数据的时候我就校验过)

校验结果


  • 好了,校验完了,没问题,那么就给大家介绍一下文件的格式。

  • 文件是什么呢,其实就是文本文件,可以直接查看,以下是示例:
    fastq格式

  • 一条read一条记录,一条记录占四行,第一行是注释,第二行是序列,就是我们所说的ATCG碱基序列,第三行是‘+’,第四行是对应的每个碱基的测序质量,也就是fastq中的“q”。每条记录之间是没有空格的。

  • 贴一下关于第一行的解释

EAS139 The unique instrument name
136 Run ID
FC706VJ Flowcell ID
2 Flowcell lane
2104 Tile number within the flowcell lane
15343 ‘x’-coordinate of the cluster within the tile
197393 ‘y’-coordinate of the cluster within the tile
1 Member of a pair, 1 or 2 (paired-end or mate-pair reads only)
Y Y if the read fails filter (read is bad), N otherwise
18 0 when none of the control bits are on, otherwise it is an even number
ATCACG Index sequence
  • 第一个是测序机器ID,独一无二的,你可以去illumina官网查到的(我没查过)

  • 第二个是这个机器跑的次数,一般机器都是有使用寿命的,所以次数越多肯定结果越差,那么一般在什么区间合适呢,我的老师给的建议是200-9999,为啥还有下限呢,那是因为,机器刚开始跑的那几次,需要磨合嘛,不是很准。看来这个展示数据不是很好啊,吓得我赶紧去看了下我的数据,还好,在212次(懒得贴图了,hexo弄个图比简书还麻烦),在范围内。

  • 倒数第三个,表示数据有没有被过滤过,一般我们的原始数据肯定是N的,要是给的原始数据是Y,你就得好好问问公司原因了。

  • 其余具体的我就不解释了,有问题可以评论交流。

  • 第二行要注意的是,可能有N的出现,那是因为有些碱基没被测出来。

  • 第四行是ASCII码表示的碱基质量,这样就能保证质量是用一个字符表示的,和碱基一一对应。公式如下:

                     P=10^-[(n-33)/10]

  • 举个栗子:比如“?”在ASCII码表上对应的编号是63,那么n就是63,减去33以后就是30,也就是我们说的Q30了,所以Q = n - 33,P算出来就是0.001,这个就是错误率,反过来,准确率就是99.9%,Q30就是准确率99.9%,同理Q20就是准确率99%。

  • ok,介绍完格式,介绍两款质控的软件,fastqc和fastp

1
2
3
4
5
6
7
8
fastqc:
$ wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.7.zip
$ unzip fastqc_v0.11.7.zip
$ cd FastQC/
$ chmod 755 fastqc
fastp:
$ wget http://opengene.org/fastp/fastp
$ chmod 755 fastp
  • 以上是我的安装方式,fastqc是一个很常用的软件,我想只要不是小白都用过。对于ubuntu用户请用我提供的方式进行安装,用apt-get install fastqc可能在使用的时候出现如下报错

fastqc错误

  • 至于fastp,是用来清洗质量不好的reads的,其实类似的软件很多,包括trimmomatic等,之所以选择这个是因为用过那么多软件后,发现fastp无论在安装还是使用上都很顺手,仅此而已。
  • 好了,现在让我们来使用fastqc对raw data的质量进行分析并查看
1
2
3
4
$ fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN
#这是fastqc的使用方法,其实很简单,对于不懂的软件,我一般的推荐是,先看官方说明,再逛逛论坛,这里不展开说这个软件怎么使用了,下面是我的使用代码
$ fastqc -o your/path/to/out -t 4 ${datadir}/*.fq.gz
#是不是超级简单,-o是报告输出的目录,这个目录必须是存在的,-t是线程数
  • 结果如下,要放在早几年,估计是很难见到这么干净的原始数据了,为此我还专门请教了一些老师,这么干净的数据的可信程度。得益于illumina公司对仪器的升级改善,数据的质量也是越来越好了。

fastqc.png

  • 但是即便如此,我们还是要对raw data进行清洗的,数据的质量直接决定了后面分析的准确性,是一切分析的前提,我们一直都在说QC,QC也是基础,但是它的重要性不容忽视。下面是我要查看的质控的内容:

    ​ • read各个位置的碱基质量值分布

    ​ • 碱基的总体质量值分布

    ​ • read各个位置上碱基分布比例,目的是为了分析碱基的分离程度

    ​ • GC含量分布

    ​ • read各位置的N含量

    ​ • read是否还包含测序的接头序列

  • 下面是我的对于clean data的指标:

    ​ 1.将含有接头的reads对去除,或者cut接头

    ​ 2.测序错误率分布检查,一般情况下,每个碱基位置的测序错误率都应该低于1%

    ​ 3.GC含量,理论上A=T,C=G,前几个碱基可能会有波动,GC含量整体在41.5%左右

    ​ 4.平均质量值分布:Q30平均比例≥80%,平均错误率≤0.1%

  • 为了达到这样的目标,我的筛选是条件是:

    ​ 去除含接头(adapter)的一对reads;可以酌情考虑去除接头保留reads

    ​ 当某一端read中的N含量超过该条read长度比例10%,去除该对reads

    ​ 当某一端read中低质量(Q≤20)碱基数超过该read长度比例的50%时,去除该对reads

1
2
3
4
5
6
7
$ fastp -i 1.fq.gz -I 2.fq.gz \
--adapter_sequence ${Adapter_R1} \
--adapter_sequence_r2 ${Adapter_R2} \#不同的试剂盒、仪器,会有不同的接头,这个可以咨询公司,也可以选择去掉这两个参数,影响不大
-o your/path/of/data/cleaned.1.fq.gz \
-O your/path/of/data/cleaned.2.fq.gz \#输出的clean data的位置和名称
-c -q 20 -u 50 -n 15 -5 20 -3 20 -w 16 \
-h your/path/of/report/clean.html -j your/path/of/report/clean.json #这俩参数分别是不同格式报告的输出位置和文件名
  • 别的不多说,解释一下我用到的几个参数:

    ​ -c:对overlap区域进行纠错,适用于PE

    ​ -q:设置低质量的标准,默认是15,多数公司这里设为5

    ​ -u:低质量碱基所占比例,默认40代表40%,只要有一条read不满足条件就成对丢掉

    ​ -n:过滤N碱基过多的reads,15代表个数,因为一般PE150的reads长度是150

    ​ -w:线程数,这里要注意!范围是1-16,建议是8,亲测,8和16差不多

    ​ -5 -3:根据质量值来截取reads,分别对应5‘端和3’端,得到reads长度可能不等

    具体请参考官网说明

  • 这个时候再看clean data结果,就无需fastqc了,因为fastp也会生成一份报告。

fastp报告部分截图

  • 看到结果,其实可以发现,数据还是保留了大部分的,所以完全不用担心cut太狠,导致数据不够。我的原则是,在数据量面前,更应该保证的是数据的准确性,即使会有一点损失,也是值得的,毕竟research更注重的也是准确性,相对于降低假阴性,我更倾向于提高真阳性。
  • 那么得到clean data以后就可以进行下一步mapping了,mapping内容下次再与大家分享。
  • 水平有限,要是存在什么错误请指出,可发送邮箱至shiyuant@outlook.com!请大家多多批评指正,相互交流,共同成长,谢谢!!!

WGS分析笔记(0)- 我的数据

发表于 2019-06-17 | 更新于 2019-06-18 | 分类于 我的WGS
  • 我研究的课题主要是通过二代重测序获得基因组数据进而进行疾病分析

概述

  • 二代测序技术(next generation sequencing,NGS),又称为高通量测序技术(high-throughput sequencing,HTS)或深度测序(deep sequencing),通过一次对几十万到几百万条核酸分子进行序列测定,进而对一个物种的转录组和基因组进行细致全貌的分析
  • 一般我们对人类基因组的研究策略主要是三个:
    • 全基因组重测序(whole genome re-sequencing,WGS)
    • 全外显子组重测序(whole exome re-sequencing,WES)
    • 目标区域测序(target region re-sequencing,TR)
  • 对于这三个策略如何应用,这完全看课题需要,没有人说过这三个策略必须分开进行,实际上目前很多课题组在进行研究的时候,都会先分析WES和panel的数据,再考虑进行WGS的测序分析,也有些课题组通过panel筛选病例,再进行WES或WGS的分析。以上只是两个简单的思路,适用于不同的课题需求。那到底怎么去选取研究策略,还是根据课题的需求,这充分说明了,课题设计的重要性,从我踩过的坑来看,一个好的课题设计,能让你后续的分析少碰很多壁,而且我觉得课题是循序渐进的,具体怎么去做要根据预实验好好的验证
  • WGS和WES怎么选择呢,也是看研究疾病的情况,比较而言,WES的成本会低,而目前报道来看80%的疾病都是由于exon区域的突变贡献的,WGS的成本相对较高,深度也不及WES,类似嵌合体什么的,靠WGS是检测不出来的,但是WGS的覆盖广,noncoding区域一直是有待研究的区域,WGS也有利于SV的检测,目前大多数类似文章的研究还是WES为主,但是主流的声响一直在推广WGS
  • 而我的课题就是属于没有被好好设计的那种,那么我的研究对象主要是(全是)WGS的数据了,既然这已经改变不了了,那就好好分析,以找出差异基因,争取早日脱离苦海
  • 那么之后我要分析的内容其实是以WGS的数据为对象,当然,这并没有太大影响,因为很多分析是大同小异的
  • 这里贴一张某测序公司培训时给的WGS数据分析流程,我的样本测序就是交给这家公司来完成的,这张照片比较早,大概是2018年4月时候的培训,不知道这家公司现在还是不是这么分析的,但这并不重要,目前大体的分析流程大同小异,一般发表的文章也会提供分析方法,也是值得大家借鉴的

分析流程

  • 而我后面的分析方法和流程与这里也是有差别的,但我为什么没有贴我自己的流程图呢,主要是因为我的流程细节部分还在不断的修改过程中,我后面的文章会按照我的流程进行一步步的介绍,也会抛出我遇到的问题,希望能得到大家的指点,大家共同交流,批评指正!!!

我搭了一个自己的主页

发表于 2019-06-17 | 更新于 2019-06-22
  • 2019年6月17日,天气还是挺热的,我在废了几小时搭建完这个hexo主页以后,决定开始把我的简书文章搬过来了
  • 其实我从来没有写这种东西的习惯,在一个朋友的建议下,我开始写,想把我学到的做过的东西记下来,一方面方便以后自己查阅,一方面,如果能帮助到其他人也是好的,毕竟我在学习之初,还是很痛苦的。
  • 我本身是比较懒的一个人,所以我一开始选择简书,一个直接往上写东西就行的平台,毕竟对我来说,搭建一个主页好累啊,但是简书实在太过分了,莫名其面锁了我的文章,在我申诉以后还不搭理我,这我就很气了,我的系列文章,你中间给我锁一篇算怎么回事!!!
  • 三次申诉无果后,我打算自己建立一个主页,毕竟写技术文章这种是我个人觉得是纯粹找diss的,因为我的技术也就那样,所以装乎我是不会去,而这个主页有个好处是,我没开通评论模块(实在是太麻烦了),这样大家也diss不了我了,但是如果我的文章出现错误,大家完全可以邮件联系我(shiyuant@outlook.com),我会第一时间答复并且修改
  • 希望能和大家多多交流,共同学习进步
12
TongShiyuan

TongShiyuan

13 日志
5 分类
8 标签
RSS
GitHub 简书 E-Mail
© 2019 TongShiyuan
博客全站共33.2k字
|