數(shù)據(jù)預(yù)處理綜述

數(shù)據(jù)預(yù)處理綜述摘要：當(dāng)今社會(huì)生物信息學(xué)已成為整個(gè)生命科學(xué)發(fā)展的重要組成部分，成為生命科學(xué)研究的前沿。隨著測序技術(shù)的不斷進(jìn)步，獲取基因序列的時(shí)間不斷縮短，測序分析中的關(guān)鍵步驟之一的數(shù)據(jù)預(yù)處理也變得尤為重要。本文對(duì)基因測序的主要兩種方法，數(shù)據(jù)預(yù)處理的概念及方法等方面進(jìn)行了論述。隨著技術(shù)的不斷革新我們對(duì)生物信息學(xué)的掌握將更加深入更加靈活，數(shù)據(jù)預(yù)處理技術(shù)的要求也越來越高，它在功能基因的準(zhǔn)確發(fā)現(xiàn)與識(shí)別、基因與蛋白質(zhì)的表達(dá)與功能研究方面都將發(fā)揮關(guān)鍵的作用。關(guān)鍵詞：sanger測序法，Illumina，SequencingbySynthesis，F(xiàn)ASTQC，Trimmomatic1主要的測序方法重點(diǎn)描述sanger法和以Illumina/SolexaGenomeAnalyzer的測序。Sanger法是根據(jù)核苷酸在某一固定的點(diǎn)開始，隨機(jī)在某一個(gè)特定的堿基處終止，并且在每個(gè)堿基后面進(jìn)行熒光標(biāo)記，產(chǎn)生以A、T、C、G結(jié)束的四組不同長度的一系列核苷酸，然后在尿素變性的PAGE膠上電泳進(jìn)行檢測，從而獲得可見的DNA堿基序列。原理:是利用一種DNA聚合酶來延伸結(jié)合在待定序列模板上的引物。直到摻入一種鏈終止核苷酸為止。每一次序列測定由一套四個(gè)單獨(dú)的反應(yīng)構(gòu)成，每個(gè)反應(yīng)含有所有四種脫氧核苷酸三磷酸(dNTP)，并混入限量的一種不同的雙脫氧核苷三磷酸(ddNTP)o由于ddNTP缺乏延伸所需要的3-OH基團(tuán)，使延長的寡聚核苷酸選擇性地在G、A、T或C處終止。終止點(diǎn)由反應(yīng)中相應(yīng)的雙脫氧而定。每一種dNTPs和ddNTPs的相對(duì)濃度可以調(diào)整，使反應(yīng)得到一組長幾百至幾千堿基的鏈終止產(chǎn)物。它們具有共同的起始點(diǎn)，但終止在不同的的核苷酸上，可通過高分辨率變性凝膠電泳分離大小不同的片段，凝膠處理后可用X-光膠片放射自顯影或非同位素標(biāo)記進(jìn)行檢測。DNA的復(fù)制需要：DNA聚合酶，雙鏈DNA模板，帶有3'-OH末端的單鏈寡核苷酸引物，4種dNTP(dATP、dGTP、dTTP和dCTP)。聚合酶用模板作指導(dǎo)，不斷地將dNTP加到引物的3'-OH末端，使引物延伸，合成出新的互補(bǔ)DNA鏈。如果加入一種特殊核苷酸，雙脫氧核苷三磷酸(ddNTP)，因它在脫氧核糖的3’位置缺少一個(gè)羥基，故不能同后續(xù)的dNTP形成磷酸二酯鍵。如，存在ddCTP、dCTP和三種其他的dNTP(其中一種為a-32P標(biāo)記)的情況下，將引物、模板和DNA聚合酶一起保溫，即可形成一種全部具有相同的5'-引物端和以ddC殘基為3’端結(jié)尾的一系列長短不一片段的混合物。經(jīng)變性聚丙烯酰胺凝膠電泳分離制得的放射性自顯影區(qū)帶圖譜將為新合成的不同長度的DNA鏈中C的分布提供準(zhǔn)確信息，從而將全部C的位置確定下來。類似的方法，在ddATP、ddGTP和ddTTP存在的條件下，可同時(shí)制得分別以ddA、ddG和ddT殘基為3'端結(jié)尾的三組長短不一的片段。將制得的四組混合物平行地點(diǎn)加在變性聚丙烯酰胺凝膠電泳板上進(jìn)行電泳，每組制品中的各個(gè)組分將按其鏈長的不同得到分離，制得相應(yīng)的放射性自顯影圖譜。從所得圖譜即可直接讀得DNA的堿基序列。與DNA復(fù)制不同的是sanger測序中的引物是單引物或者是單鏈。第二代DNA序列測序技術(shù)(以Illumina/SolexaGenomeAnalyzer測序?yàn)槔?核心思想：邊合成邊測序(SequencingbySynthesis)，即通過捕捉新合成的末端的標(biāo)記來確定DNA的序列基本原理：Illumina/SolexaGenomeAnalyzer測序的基本原理是邊合成邊測序。在Sanger等測序方法的基礎(chǔ)上，通過技術(shù)創(chuàng)新，用不同顏色的熒光標(biāo)記四種不同的dNTP，當(dāng)DNA聚合酶合成互補(bǔ)鏈時(shí)，每添加一種dNTP就會(huì)釋放出不同的熒光，根據(jù)捕捉的熒光信號(hào)并經(jīng)過特定的計(jì)算機(jī)軟件處理，從而獲得待測DNA的序列信息。操作流程：1)測序文庫的構(gòu)建(LibraryConstruction):首先準(zhǔn)備基因組DNA(雖然測序公司要求樣品量要達(dá)到200ng,但是GnomeAnalyzer系統(tǒng)所需的樣品量可低至100ng,能應(yīng)用在很多樣品有限的實(shí)驗(yàn)中)，然后將DNA隨機(jī)片段化成幾百堿基或更短的小片段，并在兩頭加上特定的接頭(Adaptor)。如果是轉(zhuǎn)錄組測序，則文庫的構(gòu)建要相對(duì)麻煩些，RNA片段化之后需反轉(zhuǎn)成cDNA，然后加上接頭，或者先將RNA反轉(zhuǎn)成cDNA，然后再片段化并加上接頭。片段的大小(Insertsize)對(duì)于后面的數(shù)據(jù)分析有影響，可根據(jù)需要來選擇。對(duì)于基因組測序來說，通常會(huì)選擇幾種不同的insertsize，以便在組裝(Assembly)的時(shí)候獲得更多信息。錨定橋接(SurfaceAttachmentandBridgeAmplification):Solexa測序的反應(yīng)在叫做flowcell的玻璃管中進(jìn)行，flowcell又被細(xì)分成8個(gè)Lane，每個(gè)Lane的內(nèi)表面有無數(shù)的被固定的單鏈接頭。上述步驟得到的帶接頭的DNA片段變性成單鏈后與測序通道上的接頭引物結(jié)合形成橋狀結(jié)構(gòu)以供后續(xù)的預(yù)擴(kuò)增使用。預(yù)擴(kuò)增(DenaturationandCompleteAmplification):添加未標(biāo)記的dNTP和普通Taq酶進(jìn)行固相橋式PCR擴(kuò)增，單鏈橋型待測片段被擴(kuò)增成為雙鏈橋型片段。通過變性，釋放出互補(bǔ)的單鏈，錨定到附近的固相表面。通過不斷循環(huán)，將會(huì)在Flowcell的固相表面上獲得上百萬條成簇分布的雙鏈待測片段。單堿基延伸測序(SingleBaseExtensionandSequencing):在測序的flowcell中加入四種熒光標(biāo)記的dNTP、DNA聚合酶以及接頭引物進(jìn)行擴(kuò)增，在每一個(gè)測序簇延伸互補(bǔ)鏈時(shí)，每加入一個(gè)被熒光標(biāo)記的dNTP就能釋放出相對(duì)應(yīng)的熒光，測序儀通過捕獲熒光信號(hào)，并通過計(jì)算機(jī)軟件將光信號(hào)轉(zhuǎn)化為測序峰，從而獲得待測片段的序列信息。從熒光信號(hào)獲取待測片段的序列信息的過程叫做BaseCalling，Illumina公司BaseCalling所用的軟件是IlluminasGenomeAnalyzerSequencingControlSoftwareandPipelineAnalysisSoftware。讀長會(huì)受到多個(gè)引起信號(hào)衰減的因素所影響，如熒光標(biāo)記的不完全切割。隨著讀長的增加，錯(cuò)誤率也會(huì)隨之上升。數(shù)據(jù)分析(DataAnalyzing):這一步嚴(yán)格來講不能算作測序操作流程的一部分，但是只有通過這一步前面的工作才顯得有意義。測序得到的原始數(shù)據(jù)是長度只有幾十個(gè)堿基的序列，要通過生物信息學(xué)工具將這些短的序列組裝成長的Contigs甚至是整個(gè)基因組的框架，或者把這些序列比對(duì)到已有的基因組或者相近物種基因組序列上，并進(jìn)一步分析得到有生物學(xué)意義的結(jié)果。2數(shù)據(jù)預(yù)處理的步驟及方法：1)Fastqc當(dāng)二代測序的原始數(shù)據(jù)拿到手之后，第一步要做的就是看一看原始reads的質(zhì)量。常用的工具就是fastqc我們?cè)诜?wù)器上用命令行來運(yùn)行fastqc：fastqc[-ooutputdir][--(no)extract][-ffastq|bam|sam][-ccontaminantfile]seqfile1..seqfileN-o用來指定輸出文件的所在目錄，注意是不能自動(dòng)新建目錄的。輸出的結(jié)果是.zip文件，默認(rèn)自動(dòng)解壓縮，命令里加上--noextract則不解壓縮。-f用來強(qiáng)制指定輸入文件格式，默認(rèn)會(huì)自動(dòng)檢測。-c用來指定一個(gè)contaminant文件，fastqc會(huì)把overrepresentedsequences往這個(gè)contaminant文件里搜索。contaminant文件的格式是"Name\tSequences"，#開頭的行是注釋。加上-q會(huì)進(jìn)入沉默模式，即不出現(xiàn)下面的提示：Startedanalysisoftarget.fqApprox5%completefortarget.fqApprox10%completefortarget.fq如果輸入的fastq文件名是target.fq，fastqc的輸出的壓縮文件將是target.fq_fastqc.zip。

解壓后，查看html格式的結(jié)果報(bào)告。結(jié)果分為如下幾項(xiàng):SummaryBasicoteitisticsF1巳］?bam巳m巳qu巳n口巳QualityF己im己qli己ncEQualitymtsur己mFeeb日日巳日巳qu巳n仁巳。口口1:巳rrtF己rb理日3GCccj門i已口tF巳］?s巳qli巳n己巳GCcunt巳rrtE巳QLien□巳Length口imteibutlizinS己qu己nc；己口ueli.&eti口口Lly己Im9OverrepresentedsequencesKinerContent結(jié)果分為綠色的"PASS"，黃色的"WARN"和紅色的"FAIL"。其中各項(xiàng)的意義如下:BasicStatistics1BasicstatisticsBasicStatisticsasure Value..fqFilename..fqFiletype Conventional tase i:^allsTOC\o"1-5"\h\zEncoding IlluiDina 1.5TotalSequences 19794124FilteredSequences0Sequencelength 100%GC； 452Perbasesequencequality:quality就是Fred值，-10*log10(p)，p為測錯(cuò)的概率。所以一條

reads某位置出錯(cuò)概率為0.01時(shí)，其quality就是20。橫軸代表位置，縱軸quality。紅色表

示中位數(shù)，黃色是25%-75%區(qū)間，觸須是10%-90%區(qū)間，藍(lán)線是平均數(shù)。若任一位置的下四分位數(shù)低于10或中位數(shù)低于25,報(bào)"WARN"；若任一位置的下四分位數(shù)低于5或中位數(shù)低于20，報(bào)"FAIL".,QuaityfctinQFacraf-s□!b-asicf（Ikmrtzi1.2encoding}12?4￡&7S910il1213 17IS202122232E2&272E393C3233 353637如約40i^&stlonInfwd(bp)3PerSequenceQualityScores:每條reads的quality的均值的分布，橫軸為quality，縱軸是reads數(shù)目。當(dāng)出現(xiàn)上圖的情況時(shí)，我們就會(huì)知道有一部分reads具有比較差的質(zhì)量。當(dāng)峰值小于27（錯(cuò)誤率0.2%）時(shí)報(bào)"WARN"，當(dāng)峰值小于20（錯(cuò)誤率1%）時(shí)報(bào)"FAIL”。4PerBaseSequenceContent:對(duì)所有reads的每一個(gè)位置，統(tǒng)計(jì)ATCG四種堿基（正常情況）

的分布橫軸為位置，縱軸為百分比。正常情況下四種堿基的出現(xiàn)頻率應(yīng)該是接近的，而且沒有位置差異。因此好的樣本中四條線應(yīng)該平行且接近。當(dāng)部分位置堿基的比例出現(xiàn)bias時(shí)，即四條線在某些位置紛亂交織，往往提示我們有overrepresentedsequence的污染。當(dāng)所有位置的堿基比例一致的表現(xiàn)出bias時(shí)，即四條線平行但分開，往往代表文庫有bias（建庫過程或本身特點(diǎn)），或者是測序中的系統(tǒng)誤差。當(dāng)任一位置的A/T比例與G/C比例相差超過10%，報(bào)"WARN"；當(dāng)任一位置的A/T比例與G/C比例相差超過20%，報(bào)"FAIL"。5PerBaseGCContent:對(duì)所有reads的每個(gè)位置，統(tǒng)計(jì)GC含量。如果建庫足夠均勻，reads的每個(gè)位置應(yīng)當(dāng)是沒有差異的，所以GC含量的線應(yīng)當(dāng)平行于X軸，反映樣品（基因組、轉(zhuǎn)錄組等）的GC含量。當(dāng)部分位置GC含量出現(xiàn)bias時(shí)，往往提示我們有overrepresentedsequence的污染。當(dāng)所有位置的GC含量一致的表現(xiàn)出bias時(shí)，往往代表文庫有bias（建庫過程或本身特點(diǎn)），或者是測序中的系統(tǒng)誤差。當(dāng)任一位置的GC含量偏離均值的5%時(shí)，報(bào)"WARN"；當(dāng)任一位置的GC含量偏離均值的10%時(shí)，報(bào)"FAIL”。6PerSequenceGCContent：統(tǒng)計(jì)reads的平均GC含量的分布。紅線是實(shí)際情況，藍(lán)線是理論分布（正態(tài)分布，均值不一定在50%，而是由平均GC含量推斷的）。曲線形狀的偏差往往是由于文庫的污染或是部分reads構(gòu)成的子集有偏差（overrepresentedreads）。形狀接近

正態(tài)但偏離理論分布的情況提示我們可能有系統(tǒng)偏差。偏離理論分布的reads超過15%時(shí)，報(bào)"WARN"；偏離理論分布的reads超過30%時(shí)，報(bào)"FAIL"。7PerBaseNContent:當(dāng)測序儀器不能辨別某條reads的某個(gè)位置到底是什么堿基時(shí)，就會(huì)產(chǎn)生“N”。對(duì)所有reads的每個(gè)位置，統(tǒng)計(jì)N的比率。正常情況下N的比例是很小的，所以圖上常常看到一條直線，但放大Y軸之后會(huì)發(fā)現(xiàn)還是有N的存在，這不算問題。當(dāng)Y軸在0%-100%的范圍內(nèi)也能看到“鼓包”時(shí)，說明測序系統(tǒng)出了問題。當(dāng)任意位置的N的比例超過5%，報(bào)"WARN"；當(dāng)任意位置的N的比例超過20%，報(bào)"FAIL”。SequenceLengthDistribution:reads長度的分布。當(dāng)reads長度不一致時(shí)報(bào)"WARN"；當(dāng)有長

度為0的read時(shí)報(bào)“FAIL”。DuplicateSequences:統(tǒng)計(jì)序列完全一樣的reads的頻率。測序深度越高，越容易產(chǎn)生一定程度的duplication，這是正常的現(xiàn)象，但如果duplication的程度很高，就提示我們可能有bias的存在（如建庫過程中的PCRduplication）。橫坐標(biāo)是duplication的次數(shù)，縱坐標(biāo)是duplicatedreads的數(shù)目，以u(píng)niquereads的總數(shù)作為100%。上圖的情況中，相當(dāng)于uniquereads數(shù)目?20%的reads是觀察到兩個(gè)重復(fù)的，?7%是觀察到三次重復(fù)的，依此類推。可以想象，如果原始數(shù)據(jù)很大（事實(shí)往往如此），做這樣的統(tǒng)計(jì)將非常慢，所以fastqc中用fq數(shù)據(jù)的前200,000條reads統(tǒng)計(jì)其在全部數(shù)據(jù)中的重復(fù)情況。重復(fù)數(shù)目大于等于10的reads被合并統(tǒng)計(jì)，這也是為什么我們看到上圖的最右側(cè)略有上揚(yáng)。大于75bp的reads只取50bp（不知道怎么選的）進(jìn)行比較。但由于reads越長越不容易完全相同（由測序錯(cuò)誤導(dǎo)致），所以其重復(fù)程度仍有可能被低估。當(dāng)非unique的reads占總數(shù)的比例大于20%時(shí)，報(bào)"WARN"；當(dāng)非unique的reads占總數(shù)的比例大于50%時(shí)，報(bào)"FAIL“。OverrepresentedSequences:如果有某個(gè)序列大量出現(xiàn)，就叫做over-represented。fastqc的標(biāo)準(zhǔn)是占全部reads的0.1%以上。和上面的duplicateanalysis一樣，為了計(jì)算方便，只取了fq數(shù)據(jù)的前200,000條reads進(jìn)行統(tǒng)計(jì)，所以有可能over-representedreads不在里面。而且大于75bp的reads也是只取50bp。如果命令行中加入了-ccontaminantfile，出現(xiàn)的over-representedsequence會(huì)從contaminant_file里面找匹配的hit（至少20bp且最多一個(gè)mismatch），可以給我們一些線索。當(dāng)發(fā)現(xiàn)超過總reads數(shù)0.1%的reads時(shí)報(bào)”WARN"，當(dāng)發(fā)現(xiàn)超過總reads數(shù)1%的reads時(shí)報(bào)”FAIL"。11OverrepresentedKmers:如果某k個(gè)bp的短序列在reads中大量出現(xiàn)，其頻率高于統(tǒng)計(jì)期望的話，fastqc將其記為over-representedk-mer。默認(rèn)的k=5，可以用-k--kmers選項(xiàng)來調(diào)節(jié)，范圍是2-10。出現(xiàn)頻率總體上3倍于期望或是在某位置上5倍于期望的k-mer被認(rèn)為是over-representedofastqc除了列出所有over-representedk-mers，還會(huì)把前6個(gè)的perbasedistribution畫出來。當(dāng)有出現(xiàn)頻率總體上3倍于期望或是在某位置上5倍于期望的k-mer時(shí)，報(bào)”WARN“；當(dāng)有出現(xiàn)頻率在某位置上10倍于期望的k-mer時(shí)報(bào)"FAIL"。)22222222)222222222）Trimmomatic它是一個(gè)針對(duì)Illumina高通量測序的readstrim的工具。即能夠針對(duì)paired-end也能弄singleended.它能夠利用FASTQ文件（phred+33或者是phred+64堿基質(zhì)量格式，取決于Illumina測序的機(jī)器）.對(duì)于single-ended，一個(gè)輸入文件和一個(gè)輸出文件，加上參數(shù)。對(duì)于paired-end數(shù)據(jù)，兩個(gè)輸入文件，4個(gè)輸出文件，分別為2個(gè)是'paired'，2個(gè)是'unpaired'（一個(gè)為forward的,一個(gè)為reverse的）。Trimmomatic用兩種策略來去除adapter:PalindromeandSimplesimpletrimming是利用每一個(gè)adapter序列去跟reads匹配，如果匹配上，就刪除read的這部分序列。Palindrometrimming是在adapter序列中re