用R语言可视化Sanger测序峰图

写在前面的
最近我师妹天天倒腾一个看峰图的软件。我很是纳闷，她这个软件倒腾过去倒腾过来干啥。跟她打听之后我才知道原来师妹上周和上上周做了很多咽拭子样本的Sanger测序，这不公司返回测序结果了，她忙着给每一个样本的测序结果进行整理、作图以发送给老师。由于样本多，每个样本并非双端都测通了。所以师妹在此花了很多时间。我后来告诉她，R语言可以画峰图，而且效果比她用的那个软件好。因此，这次推送我专门介绍一种基于R语言可视化Sanger测序峰图的方法。

Sanger测序背景

        Sanger测序也叫双脱氧链终止法测序。其原理为向合成中的核酸引入一种双脱氧核苷三磷酸（ddNTP），ddNTP可以阻止核酸扩增。每一次序列测定由一套四个单独的反应构成，每个反应含有所有四种大量的dNTP，和一种少量的ddNTP。由于ddNTP缺乏延伸所需要的3-OH基团，使延长的寡聚核苷酸选择性地在G、A、T或C处终止。终止点由反应中相应的ddNTP而定。每一种dNTPs和ddNTPs的相对浓度可以调整，使反应得到一组长几百至几千碱基的链终止产物。它们具有共同的起始点，但终止在不同的的核苷酸上，可通过高分辨率变性凝胶电泳分离大小不同的片段。
        DNA的复制需要：DNA聚合酶，双链DNA模板，带有3′-OH末端的单链寡核苷酸引物，4种dNTP（dATP、dGTP、dTTP和dCTP）。聚合酶用模板作指导，不断地将dNTP加到引物的3′-OH末端，使引物延伸，合成出新的互补DNA链。当遇到一种特定ddNTP时，由于在脱氧核糖的3’位置缺少一个羟基，故不能同后续的dNTP形成磷酸二酯键。如，存在ddCTP、dCTP和三种其他的dNTP（其中一种为α-32P标记）的情况下，将引物、模板和DNA聚合酶一起保温，即可形成一种全部具有相同的5′-引物端和以ddC残基为3’端结尾的一系列长短不一片段的混合物。经变性聚丙烯酰胺凝胶电泳分离制得的放射性自显影区带图谱将为新合成的不同长度的DNA链中C的分布提供准确信息，从而将全部C的位置确定下来。类似的方法，在ddATP、ddGTP和ddTTP存在的条件下，可同时制得分别以ddA、ddG和ddT残基为3‘端结尾的三组长短不一的片段。将制得的四组混合物平行地点加在变性聚丙烯酰胺凝胶电泳板上进行电泳，每组制品中的各个组分将按其链长的不同得到分离，制得相应的放射性自显影图谱。从所得图谱即可直接读得DNA的碱基序列。与DNA复制不同的是sanger测序中的引物是单引物或者是单链。下图对Sanger测序的基本原理进行了一个总结，如果读者是本科生可以多看看此部分，因为这是期末考试或者考研的重点。如果读者是研究生可以仔细了解一下公司返回的Sanger测序文件。由于目前此技术已经商业化，公司返回的Sanger测序结果是一个记录了测序核酸每个位点碱基读取信号值。当使用DNAstar软件查看峰图时即可看见图2的效果。

图1.Sanger测序基本原理

图2.使用DNAstar软件查看的Sanger测序峰图

使用R语言可视化Sanger测序峰图

        由于公司返回的Sanger测序结果是记录了每个位点碱基读取信号的文件，既然是文件就可以被解析。因此，使用R语言可视化峰图是可行的。我们只需要让R语言读取测序文件即可。一般情况下，如果测序的片段较长，公司会返回测序片段的正向引物测序和反向引物测序的结果。下面我会逐一介绍如何使用R语言可视化峰图及R语言相比于DNAstar的优越性。

任意方向测序结果可视化

        无论是正向引物测序结果还是反向引物测序结果，都可以使用R语言的程序包sangerseqR简单处理。处理的代码如下。（这些代码是直接复制粘贴就能使用的，读者只需替换代码中路径和文件名即可，另外读者还需安装sangerseqR的程序包）

setwd("D:/师妹的病毒测序数据/")
library(sangerseqR)
Virussangerseq<-read.abif("0004_31420082001160_(144-R131)_[144-F].ab1")
Virussangerseq<-sangerseq(Virussangerseq)
hetcalls<-makeBaseCalls(Virussangerseq,ratio = 0.10)
chromatogram(hetcalls, width = 100, height = 2,showcalls = "both")

程序运行好后会输出如下图片

图3.使用R语言简单处理Sanger测序峰图

反向引物处理

        直接对Sanger测序峰图可视化很简单，这似乎和用DNAstar并无区别。但是有一种情况可以凸显R语言优势。当正向引物扩增不成功，即公司返回的测序结果中只有反向引物成功测序了，此时你又想把测序结果调整为与DNA链方向一直的结果时，R语言能很好满足这种需求。代码如下，读者只需粘贴使用即可（路径和文件名需要适当修改调整）

setwd("D:/师妹的病毒测序数据/")
library(sangerseqR)
Virussangerseq<-read.abif("0004_31420082001160_(144-R131)_[144-R].ab1")
Virussangerseq<-sangerseq(Virussangerseq)
SeqRC=Virussangerseq
SeqRC@primarySeq<-complement(Virussangerseq@primarySeq)
SeqRC@secondarySeq<-complement(Virussangerseq@secondarySeq)
SeqRC@traceMatrix<-Virussangerseq@traceMatrix[nrow(Virussangerseq@traceMatrix):1,4:1]
a<-t(apply(Virussangerseq@peakPosMatrix,1,function(x){x[1]<-abs(nrow(Virussangerseq@traceMatrix)-x[1])+1;x[2]<-abs(x[2]-nrow(Virussangerseq@traceMatrix))+1;return(x)}))
SeqRC@peakPosMatrix<-a[nrow(a):1,]
SeqRC@peakAmpMatrix<-Virussangerseq@peakAmpMatrix[nrow(Virussangerseq@peakAmpMatrix):1,4:1]
hetcalls<-makeBaseCalls(SeqRC,ratio = 0.10)
chromatogram(hetcalls, width = 100, height = 2,showcalls = "both")

图4.使用R语言反向互补反向引物峰图

写在后面的
本期推送属于技术分享。因为我考虑到Sanger测序主要在做分子生物学相关实验的学生群体中使用最多。而做实验的学生大多数不会R语言。因此我在这篇推送中直接放出了可以使用的代码，也没有对代码的含义进行解释。这是考虑到以实验为主的学生没有相关基础。而做生物信息学的学生平时主要接触的高通量测序，这种技术资料这些学生群体应该不会特别关注。