生物信息学与基因组学1 生物信息学科研 博士

测序的世界

2018-06-19  本文已影响251人  刘小泽

测序的世界很奇妙，不同的数据处理可能得出不同的结论，入门生信首先要做的就是了解你的数据
还等什么？跟我一起来探索吧～

测序原理：

早期如何测序呢？

如果给你这么一串珠子，让你统计绿色，白色，红色，其他这四种类型的数量，那你是不是要从某个位置一个一个数？没错，你已经开始对这串珠子进行了肉眼测序。

肉眼测序.png

其实对于DNA测序，思路也是这样。

早在1954年，Whitfeld等就提出了测定多聚核糖核苷酸链的降解法，该方法利用磷酸单酯酶的脱磷酸作用和高碘酸盐的氧化作用从链末端逐一分离寡核糖核苷酸并测定其种类。目的就是想通过这种一个一个“数”的方法来得到DNA的碱基顺序。

这里再补充一个小知识，DNA有多大？它的直径也就2nm，两个核苷酸之间的小沟0.34nm。肉眼观察肯定不行，那么显微镜呢？光学显微镜可见光的波长约为700nm，单单一个直径就看不到，更别说看更小的分子了DNA对于它来说根本不可见。要想用普通显微镜观察不可能，更强大的电子显微镜呢？别忘了，DNA分子并不是平躺着让你去观察，他有四级结构，不断扭曲、折叠，简单的二维图片是很难区分开来的，能看到也仅仅只能看到一小部分。就像上面的这张珠子图，拍一张照片，就能准确区分被压着的珠子和它上面的珠子顺序吗？不能够！

人类基因组共有31.6亿个DNA碱基对，23对染色体，测一条染色体就需要测1-2亿个碱基对，像这种数数的办法肯定行不通，太耗费人力。那么有没有什么办法可以实现半自动化呢？

视觉不行，另辟蹊径，听过“盲人摸象”的故事吧，70年代的Sanger发明了“双脱氧终止反应法”，他就是利用了双脱氧核苷酸 ddNTP去摸索DNA分子。例如，一条序列5‘ --> 3' 实际为TGACTTCG但我们之前不知道，这样操作：

操作过程可能会有些无聊，但是不难懂，希望你能理解其中的意思，因为后来测序都是受此启发，这是鼻祖

A: 其实也不用刻意去除，最后利用凝胶电泳和放射自显影只能看到带有荧光标记的ddNTP，他们的排列顺序先利用电泳条带前后关系确定下，再用A-T, T-A, C-G, G-C关系反转一下，就能知道我们的测序序列

当然，千万别小瞧了一代测序（Sanger测序），到今天这种方法还是经常使用。他之所以不是现在的主流，不是因为它测不准（恰恰相反，第一代测序技术的主要特点就是测序读长可达1000bp，准确性高达99.999%，二三代所不能及），而是因为它的通量低，成本高。
要知道，目前二代测序的错误率还在
目前一代测序在验证序列（就是平时送公司测序返回来自己blast的那些）以及验证基因组组装完整性方面都是金标准。

二代测序：

一代测序的原理我们搞清楚了，没有PCR，利用聚合酶延伸链。这就有一个问题，酶的活性会下降，所以一代最长能测1000bp，再多了就要重头开始一遍导致成本高；另外它一次只测一条，也就是所谓的通量低。但是呢，他的确有值得借鉴的地方，就是准确度很高，99.999%。人类基因组计划HGP，1990年-2003年完成，利用的就是改进的Sanger测序法

后来技术不断改进，Roche公司的454技术、illumina公司的Solexa/Hiseq技术和ABI公司的SOLID技术标志第二代测序技术诞生。其中Roche公司的454测序系统是第二代测序技术中第一个商业化运营的测序平台。

其中Illumina市场规模占到75%以上，主打Hiseq，下面👇就主要介绍他的PE（Pair End双端）测序原理：

边合成边测序

<图一>第一步： 构建DNA文库
超声波将DNA分子打断成300-800bp长序列片段（人类基因组打成300-500bp），用酶补平为平末端，然后3‘端加一个A碱基（因为接头的3‘端有一个突出的T），再在两端加上互补配对的adapter，再通过PCR扩增达到一定浓度，构成单链DNA文库。
【接头设计很巧妙，两个作用：
1. 实现桥式扩增，高效；2. 可以实现双端测序】

好奇，接头怎么加上去的呢？来源https://www.fimm.fi/en/services/technology-centre/sequencing/next-generation-sequencing/dna-library-preparation

Adding Adapter

另外，利用低循环扩增技术在接头处进行修饰加上一些周边。

Adapter修饰

<图二>第二步： 上样 【重点搞清楚lane上有哪两种接头，待测序列含有哪两种接头，这很重要！】

<图三>第三步：桥式PCR

桥式形成

<图四>第四步：测序

如何确定上面👆形成的cluster的碱基排序顺序呢？illumina采取了“一次加一个荧光碱基，用完失效”的办法。官网给出的解释如下图：【有没有感觉和Sanger的方法很像？illumina的测序就是在Sanger基础上加上了桥式PCR，能克服Sange低通量的缺点】 测序
边合成变测序～测序

一轮测序是这样完成的：

双端测序之Forward Strand

先是primer结合到靠近p5的sequencing primer binding site1上，再加入特殊的dNTP【它的3‘ 羟基被叠氮基团替代，因此每次只能添加一个dNTP；还含有荧光基团，能激发不同颜色】；

在dNTP被添加到合成链上后，所有未使用的游离dNTP和DNA聚合酶会被洗脱掉；

再加入激发荧光缓冲液，用激光激发荧光信号，光学设备记录荧光信号的记录，计算机将光学信号转化为测序碱基

这一个循环就能测定flowcell上成千上万的cluster，这就实现了高通量

再向下一轮：

再加入化学试剂淬灭荧光信号并使dNTP 3’ 叠氮基团变成羟基，这样能继续向下进行再加一个，并且保证这个不再发出荧光。如此重复直至所有链的碱基序列被检测出。得到了Forward Strand序列

因为一个cluster的序列是一样的，所以理论上cluster的荧光颜色应该一致，下面是来自网站http://tucf-genomics.tufts.edu/home/faq的扫描图片。本来仪器得出的是黑白的，颜色是后来计算机计算分析后加上去的

扫描图片

Index测序： 上面的循环结束后，read product被冲掉，index1 primer和链上的index1 互补配对，进行index1的检测。测完后，洗脱产物，得到index1 的序列。接下来p5与lane上的p5‘配对，测得了index2，并洗脱

Index测序

双端测序之Reverse Strand： 洗脱掉index2 产物后，还是一个桥式扩增，得到双链，再变性得到原始Forward strand 和 新的Reverse Strand， 除去测完的Forward strand。然后和测Forward一样，也是先连接primer，只是连接的位点是Primer Binding Site2，测完后得到reverse strand序列

一个小补充：知道了PE Seq,那么单端SE 测序怎么测的呢？single-end只将index，Primer binding site以及P7/P5添加到 fragamented DNA片段的一端，另一端直接连上P5/P7，将片段固定在Flowcell上桥式PCR生成DNA簇，然后单端测序读取序列

一点小问题:illunima一次只添加一个dNTP，确实能够保证准确测量同聚物的长度问题（同聚物就是3‘端一个一个加dNTP聚合而成的聚合物，因为还带接头，所以不能直接说成DNA分子）。
那么为何illumina会限制合成的链的长度呢，不能像Sanger法一样，最长测1k？
原因就出在二代测序多出来的PCR过程：每一个位点都要测许多次，比如一段时间后的PCR得到的每个cluster都各包含200条完全相同的序列，那就需要对这200条序列的同一个位点进行测序。
第一轮我们来测第一个位点（假设位点1是A）正常来讲，200条序列都应该加A碱基，但是不巧只有199个在位点1都加了碱基A，有一条序列没有加上，所以就出现了199个红色1个灰色【当然目前还构不成影响】；
第二轮（假设位点2是G）大家应该都加G测得绿色，但是之前的那个没有加上A的，他要对之前的失误进行补偿，因此别的序列加G的时候，它加上了本该上次就加的A，它得到了红色，这个红色在一大群的绿色中就是作为杂信号存在的。依次向下，测序长度越长，杂信号越多，最后可能标准信号和杂信号各一半，这样系统无法判断，只能给N，而N多了对于后续的分析处理很麻烦，去了吧丢失数据，不去吧又是冗余。

Q:有同学问了: 上面说的，依次向下，测序长度越长，杂信号越多。这个说不通啊，依次向下以后，永远都只有1/200个杂信号啊，怎么会越来越多呢？
A:上面提到的只是假设只是一条序列出现失误，其实测序过程中随着链的延长大家的错误都是不断增加的，因为酶活性在降低。对于一个测序位点，第一轮可能有一条序列错的，那么肯定他以后都会错吧；到了第二轮可能有两条或者更多的序列出现没加的现象，而他们也会一步错，步步错。

因此，错误来源就是某个位点有序列没跟上队伍，这次不加，以后都是错的，这样的序列多了，也就导致整体错误率升高

数据产生：

大体上我们平常使用的测序仪就是这样（以Hiseq 2000为例），后续升级版主要提高通量

Hiseq 2000测序仪 相机识别 illumina测序仪产出统计

数据初步分析：

使用fastqc进行质量分析，这是一款Java软件，支持多线程。写这篇文章的时候版本是v0.11.7。

软件前期准备：

  1. 官网下载好用filezilla导入linux服务器

  2. 直接在服务器中wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.7.zip

最后这个chmod有必要说一下，这个权限管理命令

chmod 用3个数字来表达对 用户（文件或目录的所有者），用户组（同组用户），其他用户 的权限：
　　如：chmod 755 fastqc
　　数字7是表达同时具有读，写，执行权限：(7 = 4 + 2+ 1)
　　读取--用数字4表示；
　　写入--用数字2表示；
　　执行--用数字1表示; 三者皆否：0

后期软件使用：
   基本格式 （各种参数+多个文件～支持多线程）
   
   fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] seqfile1 .. seqfileN
   
   -o --outdir FastQC生成的报告文件的储存路径
   --extract 使用这个参数是让程序不打包【默认会打包成一个压缩文件】
   -t --threads 选择程序运行的线程数，每个线程会占用250MB内存（一般与文件数量一致就好）
   -q --quiet 安静运行模式【不选这个选项，程序会实时报告运行的状况】
结果分析：
# 顺便说一下这里用到了curl -O（保留远程文件的文件名） -L（对于自动跳转的网页，curl 就会跳转到新的网址）
# 当然用wget也可以，至于二者区别，可以参考https://www.cnblogs.com/lsdb/p/7171779.html
    curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_1.fastq.gz
    curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_2.fastq.gz
    curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1977249_1.fastq.gz
    curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1977249_2.fastq.gz
    
# 下载完可以检查一下数据完整性，当然不是必须
    md5sum *.gz
    
# 质控四个文件，我们可以采用四线程
# 大概用时   real   0m23.344s （如果你也想统计时间，就在命令前加time）
    fastqc *.gz -t 4
#将结果.html用filezilla导出，浏览器查看

  1. Basic Statistics 基本信息

    Basic Statistics

    • Encoding: 测序平台编号，现在Sanger/ Illumina 1.8以上都是Phred 33编码

    • Total sequences: reads数量（reads就是高通量测序平台产生的序列标签，翻译为读段！）

    • Sequence length: 测序长度

    • %GC: GC含量： 需要重点关注，可以帮助区别物种，人类细胞42%左右

  2. Per base sequence quility：每个测序read各碱基质量【十分重要！】

Per base sequence quility
  1. Per sequence quility scores：每条序列 质量统计
Per sequence quility scores
  1. Per base sequence content：read各个位置碱基比例分布
Per base sequence content
  1. Per sequence GC content： 序列平均GC分布
Per sequence GC content
  1. Per base N content： N含量分布
Per base N content
  1. Sequence length distribution： 序列长度统计
Sequence length distribution
  1. Sequence duplication level：统计序列完全一样的reads的频率
Sequence duplication level
  1. Overrepresented sequences：大量重复序列
Overrepresented sequences

  1. Adapter content: 接头含量

    Adapter content
  1. （还有一类这里没体现）Kmer content： 重复短序列
Kmer content

数据过滤：

主要针对接头序列和低质量序列

java -jar <path to trimmomatic.jar> PE [-threads <threads 线程数] 
        [-phred33 | -phred64] 质量体系，默认64，但我们现在一般都是33
        [-trimlog <logFile>] log文件
        <input 1> <input 2> 双端测序原始fq文件
        <paired output 1> <unpaired output 1>  双端1输出文件 、 过滤掉的文件
        <paired output 2> <unpaired output 2>  双端2同理
        <step > 
        <详细讲一下step：>
        1. ILLUMINACLIP： 根据上面qc部分的测试数据，我们设置TruSeq2-PE.fa:2:40:15
         TruSeq2-PE.fa是接头序列；
         2是比对时接头序列时所允许的最大错误数；
         40是要求PE的两条read同时和adapter序列比对，匹配度加起来超40%，那么就认为这对PE reads含  有adapter，并在对应的位置需要进行切除；
        【那么为何不是匹配100%？因为测序时并不是把 adapter全测了，只测了一部分】
         15 指的是只要某条read的某部分与adapter超过了15%的相似度就认为包含，就要去除
        2. SLIDINGWINDOW： 滑动窗口长度 我们设置为4:20，代表窗口长度为4，窗口中的平均质量值至少为20，否  则会开始切除
        3. LEADING，指的是read开头的碱基是否要被切除的质量阈值，这里设为2
        4. TRAILING，指的是read末尾的碱基是否要被切除的质量阈值，这里设为2
        5. MINLEN，指的是read被切除后至少需要保留的长度，如果低于该长度，会被丢掉，这里设置25
        ​

#一个☝️范例(这个测试数据是phred 64的，所以使用默认值就好)：
        java -jar trimmomatic-0.36.jar PE -threads 8 \
        -trimlog logfile \
        reads_1.fq.gz reads_2.fq.gz \
        out.read_1.fq.gz out.trim.read_1.fq.gz \
        out.read_2.fq.gz out.trim.read_2.fq.gz \  ILLUMINACLIP:/path/Trimmomatic/adapters/TruSeq2-PE.fa:2:40:15 \
        SLIDINGWINDOW:4:20 \
        LEADING:2 TRAILING:2 \
        MINLEN:25</pre>

# 当然，这只是对一个样本的双端测序文件进行操作，那么问题来了：
        # 如果你的样本比较多呢？还要手动一个个输入吗？
        # 虽然说vim中使用快速替换 :%s/AA/BB/g 可以全局替换AA成BB，但还是有点麻烦
        # 能不能让程序来一个自动化呢？是可以的！可以看作是上面👆脚本的改进版～
        # 在脚本中构建for循环
        # vi trim.sh
        ​
         #开头这样写而不写*.gz的目的是避免了样本名称的重复
        for filename in *_1.fastq.gz  
        do
         #对于filename（%）向右匹配_的全部内容(*),然后这部分去掉，留下这之前的
        base=${filename%_*}  
        java -jar trimmomatic-0.36.jar PE \       # 加上反斜线能让程序整体更清晰
         -threads 8 \
         -trimlog logfile \
         ${base}_1.fastq.gz ${base}_2.fastq.gz \
         out.${base}_1.fastq.gz out.trim.${base}_1.fastq.gz \
         out.${base}_2.fastq.gz out.trim.${base}_2.fastq.gz \
         ILLUMINACLIP:/home/u1239/biosoft/Trimmomatic-0.36/adapters/TruSeq2-PE.fa:2:30:10 \
         SLIDINGWINDOW:5:20 \
         LEADING:5 TRAILING:5 \
         MINLEN:50
        done
多样本trim

多样本质控：

fastqc对于一两个样品还能吃得消，但样本多了，我们如何同时查看对比他们的信息呢？

这里就可以使用界面更加优美的multiqc软件了, 他就是fastqc结果的整合我会从一个初学者角度来一步步进行，其中包括了中间犯的错误及改正～，如果你也在运行的过程发现了类似的错误，可以参考一下。

分界线1: 尝试自己安装

下载地址：

https://files.pythonhosted.org/packages/fa/3e/fbdcbaebadad2110eed3b4212ced58617cbd9badbfc728c5eabc53916b84/multiqc-1.5.tar.gz

解压缩后，会发现和以前安装的源代码文件不同，他没有直接显示可执行文件，也没有配置文件。这是因为multiqc是一个python软件，这里先看一下setup.py：

MultiQC setup.py

直接使用python setup.py会报错，因为可能你的服务器并没有setuptools，这是python依赖的第三方库。

先安装setuptools：

下载地址：

https://files.pythonhosted.org/packages/1a/04/d6f1159feaccdfc508517dba1929eb93a2854de729fa68da9d5c6b48fa00/setuptools-39.2.0.zip

解压缩完setuptools会看到有这些文件，这也是一般的python软件常见的：包括了setup.py，easy_install.py等

一般的Python程序 
# 主要使用setup.py
python setup.py build # 编译
python setup.py install #安装

在第二步安装的时候报错,原因是不能对/usr/local/lib下的python进行操作，因为不管/usr/bin还是/usr/local/bin，都是可读不可改，如果自己家目录环境变量中查不到python的路径，那么安装过程中会自动调取更上层目录的python使用。这往往会引发一系列问题。如果要自己更改的话，需要自己home目录下有python

报错信息
好吧，那么接下来我们再安装自己的python，毕竟自己的软件用起来方便：
# 下载地址 
wget https://www.python.org/ftp/python/3.6.5/Python-3.6.5.tgz
# .tgz 就等同tar.gz
# 解压后会看到配置文件configure，设置成自己的软件环境变量即可
./configure --prefix=/YOUR PATH/
make
make install
# 可能结尾会有报错，但是不影响使用，出现了可执行文件python就成功了
# 将python复制到你的环境变量中就能直接调取
# 最后使用which python检测是否安在了自己的目录下

自己的python安装完成了，进入到解压好的multiqc文件夹下, python setup.py build 编译会运行成功，然后python setup.py install

安装成功

接着，multiqc就会出现在了自己的软件环境变量中了，输入multiqc就会看到

multiqc

分界线2: 同一件事，换个玩法

目的：安装multiqc
途径：conda自动安装
缘由：好久不用conda安装软件，一直坚持源代码，因为好掌控。但是有的软件依赖的包很多，自己又很难发现这些潜在的关系，所以想重新试试conda，但并不是傻瓜式的使用，而是让conda听我的话～授之以鱼，不如授之以渔

Here we go！

安装： 安装过程需要注意

  1. 可以自定义新文件夹，默认是在home/miniconda3/

  2. ⚠️不要将conda添加到PATH中，我们只需要把它视作一个保姆babysitter，而不要当一个管家housekeeper

  3. 安装完conda，把miniconda/bin下的conda复制到你的/biotools/PATH（这个环境变量因人而异）

添加源： 清华源和中科大源都不错，别忘了再添加一个bioconda源

新建conda专属下载目录： 你可以在你的biotools目录下新建一个conda软件存放目录，例如conda_install。然后把这个文件夹添加到环境变量。以后你用conda安装的所有软件都存放在这里，

conda专属目录

⚠️和你自己安装的软件要区分开， 然后利用conda install -p /PATH/conda_install multiqc 就可以实现multiqc的安装了

写在后面

  1. multiqc的使用很简单，然后结果还是交互式的，能导出很多格式

  2. 这次实现了人工和自动二者的有机结合：复杂的，多依赖性的软件还是要靠conda解决。但是如果一个软件用两种方法都安装了，怎么选择？这个利用bashrc中设置的顺序：哪个目录优先，就先搜索哪个目录

结果一目了然：trim掉低质量序列和接头后，新的数据质量值都在q30以上
MultiQC 1 MultiQC 2

欢迎关注我们的公众号～_～　　
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到Bioplanet520@outlook.com

bioinfoplanet
上一篇下一篇

猜你喜欢

热点阅读