2013年左右，二代测试开始崛起，不到一年大有横扫测序行业的势头。记得当时二代测序服务公司的业务员到我们实验室推销业务，宣传册十分“高大上”的样子，一问价格，内心直哆嗦。但当时实验室有ABi 3500的机器，正玩的飞起，根本没在意二代测序，不过世间万物逃不过“真香”，研三下学期已经对二代十分感兴趣，后来还做了一段二代测序试剂研发。
一转眼，从上家公司离职一年半了，二代数据分析的方法和流程也忘得差不多了。最近解读一代的结果时，需要收集突变信息，就想到了能不能用二代call SNP的方法来跑跑一代的数据。跑出来的结果还是没有肉眼直观，好在程序能复用，也可以复习下二代软件的使用，记录一下过程。
基本流程
- ab1数据转换成fastq文件
- 质控fastq文件
- 构建参考序列index
- 序列比对
- 突变位点分析
软件准备
- extract_fastq.exe
该程序属于Staden Package软件包的一个组件，Staden Package是sanger测序结果分析、拼接、编辑的强大工具，感兴趣的自己研究。Staden Package有windows版和linux版，由于本系列涉及工具大多在linux平台，windows的就不介绍了。下载源码自己编译，命令三部曲：
./configure
make
make install
安装成功后，执行“extract_fastq -h”会显示帮助信息，用法为：
extract_fastq XXX.ab1 > out.fastq
windows版使用也需要在cmd命令行启动。 - cutadapt
cutadapt是python语言编写的去接头程序，如果你已经安装python和pip，那就好办了，直接："pip install cutadapt --upgrade"就搞定了。cutadapt的使用比较简单，具体方法可参考官方文档：https://cutadapt.readthedocs.io/en/stable/guide.html，可以记住几个常用参数：
cutadapt [参数] -o out.fastq in.fastq
-a "ATCG..." 去3‘端接头
-g "ATCG..." 去5’端接头
-q 20,20 切除质量Q值低于20的碱基，如果仅有一个数字仅去除3‘端低质量碱基
-t 4 线程数
cutadapt可用于单端fastq文件质控也可用于双端配对fastq质控，我们用ab1转换成的fastq一般按单端处理，如果是双向测通也可以当作双端测序文件，命令格式如下：
cutadapt -a ADAPTER_FWD -A ADAPTER_REV -o out.1.fastq -p out.2.fastq reads.1.fastq reads.2.fastq
-a 去fastq1的3‘端接头
-A 去fastq2（与fastq1配对）的3‘端接头（大写代表第二的fastq文件的参数）
-o 输出fastq1的处理结果
-p 输出fastq2的处理结果
cutadapt处理文件时会给出统计信息，可以看到处理了 多少reads，砍掉了多少碱基等。 - bwa
大名鼎鼎的bwa，网上教程很多，下载和使用说明在这里：http://bio-bwa.sourceforge.net/，下载后解压，同样是编译安装三部曲：
./configure
make
make install
经过extract_fastq.exe和cutadapt的处理，ab1文件已经成为可用的fastq文件，在序列比对之前需要用bwa index命令构建参考序列的索引文件，命令如下：
bwa index -a is -p 参考序列名 XXX.fasta
短序列和小型基因组用"-a is"参数，大型基因组用"-a bwtsw"参数，-p 指定index文件的名称，不指定就会以fasta文件全名做名称。 - samtools
李恒博士的另一大作，与bwa配套使用，下载及使用说明网址：http://samtools.sourceforge.net/，下载安装过程同bwa。常用命令有：
samtools view 查看bam文件
samtools sort 将bam文件的比对结果按在参考序列上的位置从前往后排序，samtools depth，samtools mpileup等依赖于该命令。
samtools depth 查看测序文件对参考序列的覆盖深度
samtools mpileup 输出测序文件对参考序列的突变情况
操作步骤
假设有参考序列A，保存成fasta格式：A.fasta/A.fa，
Sanger测序结果：A-1_xxx.ab1、A-2_xxx.ab1、A-3_xxx.ab1、...。
ab1文件保存了测序的序列、测序质量、峰图等信息，它有着特定的格式（ab1文件格式说明看这里），除序列外一般无法直接提取其他信息，我们可以使用专业软件将ab1转换为文本。
①将所有ab1文件转换为fastq格式并写入到一个文件中：
touch A.fastq #创建一个空文件用来接收转换结果
ls *.ab1 | sort | xargs 你自己的软件安装路径/extract_fastq.exe >> A.fastq
②去除低质量测序碱基，碱基Q值与准确率P的对应关系是：Q=-10×lg(1-P)，Q值越高准确率越高。Q=10，准确率=90%；Q=20，准确率=99%；Q=30，准确率=99.9%。一般99%的准确度足够用于分析，质控命令为：
cutadapt -q 20,20 -o A.cut.fastq A.fastq
③构建参考序列index：
bwa index -a is -p A A.fasta
④序列比对，使用bwa mem命令，该命令采用BWA-MEM算法，适合70bp-1Mbp的长序列比对，命令及参数为：
bwa mem -R '@RG\tID:A_all\tPL:Sanger\tLB:20200711\tSM:A' A
A.cut.fastq | samtools view -S -b - > A.bam
⑤将比对结果进行排序：
samtools sort -o A.sorted.bam --reference A.fasta A.bam
要处理bam文件，需要对其格式有所了解，可以参考这篇文章：https://www.jianshu.com/p/a584d31418f3，或者直接查看格式文档：https://samtools.github.io/hts-specs/SAMv1.pdf。其中对CIGRA式的理解十分重要，M代表匹配到参考序列但不代表与参考序列一致，可能是突变，突变详情会在MD:Z:xxxxx字段记录，详查SAM文件说明。D代表缺失突变，I代表插入突变，S代表不匹配序列。
bam文件中对统计突变信息有用的就是第4，6，13列，可以单独输出到txt文件中进一步分析：
samtools view A.sorted.bam | cut -f 4,6,13 > A.bam.txt
⑥提取序列的突变信息：
samtools mpileup -aa -t DP,AD --reference A.fasta A.sorted.bam > A.mutant.txt
要解读samtools mpileup结果需要了解其输出格式，可以参考这篇文章https://www.cnblogs.com/ywliao/articles/7657761.html。一般我们只关注前6列：
第1列是参考序列名称，第2列是在参考序列上的位置，第3列是参考序列的碱基，第4列是测序深度，第5列是每个测序位点的情况，第6列是该位点测序结果的质量值。
A 31 G 8 ..T..... XXXXXX^X
A 37 T 8 ....-1A...-1A. SSSLSSLS
A 610 G 8 ......+1C.. XXXXXHXX
A 2165 G 17 .......-1A,,,,,,-1a,,., NXNNINNXXXXXRXXNX
其中第5列最为关键，一般有多大深度该列就有多少种情况，对于突变我们只关注上面的四种情况：
- 有ATCG字符，代表某次测序该位置为其他碱基，如第一行，参考序列为G，8次测序，其中一次为T。
- 负数加碱基，代表该位置之后存在缺失，如第二行，参考序列第37bp后，缺失一个T碱基，两次测序存在缺失，其他6次正常。
- 正数加碱基，代表该位置之后存在插入，如第三行，参考序列第610bp后，有一次测序结果插入1个C碱基
- 小写碱基，代表负链突变，大写代表正链，如第四行，参考序列第2165bp后，一次测序正链缺失，一次测序负链同一位置缺失，这表明有两次测序在该位置重叠，且一为正向，一为反向。
结果分析
根据目的提取或者处理结果，需要依赖一定的编程能力，首推python。经过前面的处理，可以得到每个样本的突变位置及突变类型，利用python很容易从A.bam.txt文件或A.mutant.txt提取这样的信息。比如克隆子测序，可以得到每个克隆子的基因型。
用二代测序的分析流程，处理一代的结果，好处是可复用，可以把所有命令写入一个shell脚本，每次仅一行调用命令即可，缺点是不直观。
网友评论