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

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

解壓后，查看html格式的結果報告。結果分為如下幾項:SummaryBasicoteitisticsF1巳］?bam巳m巳qu巳n口巳QualityF己im己qli己ncEQualitymtsur己mFeeb日日巳日巳qu巳n仁巳?？诳?:巳rrtF己rb理日3GCccj門i已口tF巳］?s巳qli巳n己巳GCcunt巳rrtE巳QLien□巳Length口imteibutlizinS己qu己nc；己口ueli.&eti口口Lly己Im9OverrepresentedsequencesKinerContent結果分為綠色的"PASS"，黃色的"WARN"和紅色的"FAIL"。其中各項的意義如下: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為測錯的概率。所以一條

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

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

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

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

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