2013年左右秩仆，二代測(cè)試開始崛起码泛，不到一年大有橫掃測(cè)序行業(yè)的勢(shì)頭。記得當(dāng)時(shí)二代測(cè)序服務(wù)公司的業(yè)務(wù)員到我們實(shí)驗(yàn)室推銷業(yè)務(wù)澄耍，宣傳冊(cè)十分“高大上”的樣子噪珊，一問價(jià)格，內(nèi)心直哆嗦齐莲。但當(dāng)時(shí)實(shí)驗(yàn)室有ABi 3500的機(jī)器痢站，正玩的飛起，根本沒在意二代測(cè)序选酗，不過世間萬物逃不過“真香”阵难，研三下學(xué)期已經(jīng)對(duì)二代十分感興趣，后來還做了一段二代測(cè)序試劑研發(fā)芒填。
一轉(zhuǎn)眼呜叫，從上家公司離職一年半了，二代數(shù)據(jù)分析的方法和流程也忘得差不多了氢烘。最近解讀一代的結(jié)果時(shí)怀偷，需要收集突變信息家厌，就想到了能不能用二代call SNP的方法來跑跑一代的數(shù)據(jù)播玖。跑出來的結(jié)果還是沒有肉眼直觀，好在程序能復(fù)用饭于，也可以復(fù)習(xí)下二代軟件的使用蜀踏，記錄一下過程。
基本流程
- ab1數(shù)據(jù)轉(zhuǎn)換成fastq文件
- 質(zhì)控fastq文件
- 構(gòu)建參考序列index
- 序列比對(duì)
- 突變位點(diǎn)分析
軟件準(zhǔn)備
- extract_fastq.exe
該程序?qū)儆?a target="_blank">Staden Package軟件包的一個(gè)組件掰吕，Staden Package是sanger測(cè)序結(jié)果分析果覆、拼接、編輯的強(qiáng)大工具殖熟，感興趣的自己研究局待。Staden Package有windows版和linux版，由于本系列涉及工具大多在linux平臺(tái)，windows的就不介紹了钳榨。下載源碼自己編譯舰罚，命令三部曲：
./configure
make
make install
安裝成功后，執(zhí)行“extract_fastq -h”會(huì)顯示幫助信息薛耻，用法為：
extract_fastq XXX.ab1 > out.fastq
windows版使用也需要在cmd命令行啟動(dòng)营罢。 - cutadapt
cutadapt是python語言編寫的去接頭程序，如果你已經(jīng)安裝python和pip饼齿，那就好辦了饲漾，直接："pip install cutadapt --upgrade"就搞定了。cutadapt的使用比較簡(jiǎn)單缕溉，具體方法可參考官方文檔：https://cutadapt.readthedocs.io/en/stable/guide.html考传，可以記住幾個(gè)常用參數(shù)：
cutadapt [參數(shù)] -o out.fastq in.fastq
-a "ATCG..." 去3‘端接頭
-g "ATCG..." 去5’端接頭
-q 20,20 切除質(zhì)量Q值低于20的堿基，如果僅有一個(gè)數(shù)字僅去除3‘端低質(zhì)量堿基
-t 4 線程數(shù)
cutadapt可用于單端fastq文件質(zhì)控也可用于雙端配對(duì)fastq質(zhì)控证鸥，我們用ab1轉(zhuǎn)換成的fastq一般按單端處理伙菊，如果是雙向測(cè)通也可以當(dāng)作雙端測(cè)序文件，命令格式如下：
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配對(duì)）的3‘端接頭（大寫代表第二的fastq文件的參數(shù)）
-o 輸出fastq1的處理結(jié)果
-p 輸出fastq2的處理結(jié)果
cutadapt處理文件時(shí)會(huì)給出統(tǒng)計(jì)信息敌土，可以看到處理了 多少reads镜硕，砍掉了多少堿基等。 - bwa
大名鼎鼎的bwa返干，網(wǎng)上教程很多兴枯，下載和使用說明在這里：http://bio-bwa.sourceforge.net/，下載后解壓矩欠，同樣是編譯安裝三部曲：
./configure
make
make install
經(jīng)過extract_fastq.exe和cutadapt的處理财剖，ab1文件已經(jīng)成為可用的fastq文件，在序列比對(duì)之前需要用bwa index命令構(gòu)建參考序列的索引文件癌淮，命令如下：
bwa index -a is -p 參考序列名 XXX.fasta
短序列和小型基因組用"-a is"參數(shù)躺坟，大型基因組用"-a bwtsw"參數(shù)，-p 指定index文件的名稱乳蓄，不指定就會(huì)以fasta文件全名做名稱咪橙。 - samtools
李恒博士的另一大作，與bwa配套使用虚倒，下載及使用說明網(wǎng)址：http://samtools.sourceforge.net/美侦，下載安裝過程同bwa。常用命令有：
samtools view 查看bam文件
samtools sort 將bam文件的比對(duì)結(jié)果按在參考序列上的位置從前往后排序魂奥，samtools depth菠剩，samtools mpileup等依賴于該命令。
samtools depth 查看測(cè)序文件對(duì)參考序列的覆蓋深度
samtools mpileup 輸出測(cè)序文件對(duì)參考序列的突變情況
操作步驟
假設(shè)有參考序列A耻煤，保存成fasta格式：A.fasta/A.fa具壮，
Sanger測(cè)序結(jié)果：A-1_xxx.ab1准颓、A-2_xxx.ab1、A-3_xxx.ab1棺妓、...瞬场。
ab1文件保存了測(cè)序的序列、測(cè)序質(zhì)量涧郊、峰圖等信息贯被，它有著特定的格式（ab1文件格式說明看這里），除序列外一般無法直接提取其他信息妆艘，我們可以使用專業(yè)軟件將ab1轉(zhuǎn)換為文本彤灶。
①將所有ab1文件轉(zhuǎn)換為fastq格式并寫入到一個(gè)文件中：
touch A.fastq #創(chuàng)建一個(gè)空文件用來接收轉(zhuǎn)換結(jié)果
ls *.ab1 | sort | xargs 你自己的軟件安裝路徑/extract_fastq.exe >> A.fastq
②去除低質(zhì)量測(cè)序堿基，堿基Q值與準(zhǔn)確率P的對(duì)應(yīng)關(guān)系是：Q=-10×lg(1-P)批旺，Q值越高準(zhǔn)確率越高幌陕。Q=10，準(zhǔn)確率=90%汽煮；Q=20搏熄，準(zhǔn)確率=99%；Q=30暇赤，準(zhǔn)確率=99.9%心例。一般99%的準(zhǔn)確度足夠用于分析，質(zhì)控命令為：
cutadapt -q 20,20 -o A.cut.fastq A.fastq
③構(gòu)建參考序列index：
bwa index -a is -p A A.fasta
④序列比對(duì)鞋囊，使用bwa mem命令止后，該命令采用BWA-MEM算法，適合70bp-1Mbp的長(zhǎng)序列比對(duì)溜腐，命令及參數(shù)為：
bwa mem -R '@RG\tID:A_all\tPL:Sanger\tLB:20200711\tSM:A' A
A.cut.fastq | samtools view -S -b - > A.bam
⑤將比對(duì)結(jié)果進(jìn)行排序：
samtools sort -o A.sorted.bam --reference A.fasta A.bam
要處理bam文件译株，需要對(duì)其格式有所了解，可以參考這篇文章：http://www.reibang.com/p/a584d31418f3挺益，或者直接查看格式文檔：https://samtools.github.io/hts-specs/SAMv1.pdf歉糜。其中對(duì)CIGRA式的理解十分重要，M代表匹配到參考序列但不代表與參考序列一致望众，可能是突變匪补，突變?cè)斍闀?huì)在MD:Z:xxxxx字段記錄，詳查SAM文件說明黍檩。D代表缺失突變叉袍，I代表插入突變，S代表不匹配序列刽酱。
bam文件中對(duì)統(tǒng)計(jì)突變信息有用的就是第4，6瞧捌，13列棵里，可以單獨(dú)輸出到txt文件中進(jìn)一步分析：
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結(jié)果需要了解其輸出格式润文，可以參考這篇文章https://www.cnblogs.com/ywliao/articles/7657761.html。一般我們只關(guān)注前6列：
第1列是參考序列名稱殿怜，第2列是在參考序列上的位置典蝌，第3列是參考序列的堿基，第4列是測(cè)序深度头谜，第5列是每個(gè)測(cè)序位點(diǎn)的情況骏掀，第6列是該位點(diǎn)測(cè)序結(jié)果的質(zhì)量值。
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列最為關(guān)鍵柱告，一般有多大深度該列就有多少種情況截驮，對(duì)于突變我們只關(guān)注上面的四種情況：
- 有ATCG字符，代表某次測(cè)序該位置為其他堿基际度，如第一行葵袭，參考序列為G，8次測(cè)序乖菱，其中一次為T坡锡。
- 負(fù)數(shù)加堿基，代表該位置之后存在缺失窒所，如第二行鹉勒，參考序列第37bp后，缺失一個(gè)T堿基吵取，兩次測(cè)序存在缺失贸弥，其他6次正常。
- 正數(shù)加堿基海渊，代表該位置之后存在插入绵疲，如第三行，參考序列第610bp后臣疑，有一次測(cè)序結(jié)果插入1個(gè)C堿基
- 小寫堿基盔憨，代表負(fù)鏈突變，大寫代表正鏈讯沈，如第四行郁岩，參考序列第2165bp后，一次測(cè)序正鏈缺失缺狠，一次測(cè)序負(fù)鏈同一位置缺失问慎，這表明有兩次測(cè)序在該位置重疊，且一為正向挤茄，一為反向如叼。
結(jié)果分析
根據(jù)目的提取或者處理結(jié)果，需要依賴一定的編程能力穷劈，首推python笼恰。經(jīng)過前面的處理踊沸，可以得到每個(gè)樣本的突變位置及突變類型，利用python很容易從A.bam.txt文件或A.mutant.txt提取這樣的信息社证。比如克隆子測(cè)序逼龟，可以得到每個(gè)克隆子的基因型。
用二代測(cè)序的分析流程追葡，處理一代的結(jié)果腺律，好處是可復(fù)用，可以把所有命令寫入一個(gè)shell腳本宜肉，每次僅一行調(diào)用命令即可匀钧，缺點(diǎn)是不直觀。
