第十章 使用序列數(shù)據(jù)

生物信息學(xué)的核心問題之一是處理大量的（通常定義糟糕或模糊）文件格式洒缀。久而久之惩猫，一些特定的簡(jiǎn)單的人類可讀的格式成為了事實(shí)上的標(biāo)準(zhǔn)。

• 彼得·科克等, 2010年
  好的程序員知道怎樣寫程序宵睦。杰出的程序員知道怎樣重寫（并重新使用）吩抓。
• 大教堂和集市 埃里克·斯蒂芬·雷蒙德

在生物信息學(xué)中，核苷酸（和蛋白質(zhì)）序列普遍以兩種純文本格式存儲(chǔ)：FASTA和FASTQ鸥诽，發(fā)音分別為fast-ah（或fast-A）和fast-Q商玫。 我們將在本節(jié)討論這兩種格式和它們各自的局限性，然后了解一些處理這些數(shù)據(jù)格式的工具牡借。這是一個(gè)簡(jiǎn)短的章節(jié)拳昌，但是本節(jié)有一個(gè)重要的事項(xiàng)需要注意，也就是處理這些特定生物信息學(xué)格式的共同誤區(qū)钠龙。即使像文件格式這樣的細(xì)枝末節(jié)中存在的簡(jiǎn)單錯(cuò)誤炬藤，也需要很多時(shí)間和精力去發(fā)現(xiàn)和修復(fù)，所以請(qǐng)盡早記住這些細(xì)節(jié)碴里。

FASTA格式

FASTA格式源于由William R. Pearson和David J. Lipman開發(fā)的FASTA比對(duì)套裝軟件沈矿。FASTA格式用于存儲(chǔ)各種序列數(shù)據(jù)，不需要每個(gè)堿基對(duì)的質(zhì)量分?jǐn)?shù)咬腋。這些序列數(shù)據(jù)包括參考基因組文件羹膳，蛋白質(zhì)序列，DNA編碼序列（CDS）根竿，轉(zhuǎn)錄本序列等陵像。FASTA也可以用來存儲(chǔ)多序列比對(duì)數(shù)據(jù)就珠，但是我們這里不會(huì)討論這個(gè)FASTA的格式的特殊變體。我們?cè)谇懊娴恼鹿?jié)遇到了FASTA格式醒颖，但在本節(jié)中妻怎，我們將詳細(xì)介紹這種格式，了解該格式常見的缺陷泞歉，并介紹一些處理這種格式的工具逼侦。

FASTA文件由序列條目組成，每個(gè)條目都包含兩部分內(nèi)容：序列描述和序列數(shù)據(jù)腰耙。描述行以一個(gè)大于號(hào)（>）開始榛丢，包含序列標(biāo)識(shí)符和其他（可選）信息。 序列數(shù)據(jù)從描述行的下一行開始沟优，直到另一個(gè)描述行（以 > 開頭的行）出現(xiàn)或文件結(jié)束涕滋。 本章GitHub目錄中的 egfr_ FL ank.fasta 文件是一個(gè)FASTA例子文件：

$ head -10 egfr_flank.fasta
> ENSMUSG00000020122 | ENSMUST00000138518
CCCTCCTATCATGCTGTCAGTGTATCTCTAAATAGCACTCTCAACCCCCGTGAACTTGGT
TATTAAAAACATGCCCAAAGTCTGGGAGCCAGGGCTGCAGGGAAATACCACAGCCTCAGT
TCATCAAAACAGTTCATTGCCCAAAATGTTCTCAGCTGCAGCTTTCATGAGGTAACTCCA
GGGCCCACCTGTTCTCTGGT 
> ENSMUSG00000020122 | ENSMUST00000125984
GAGTCAGGTTGAAGCTGCCCTGAACACTACAGAGAAGAGAGGCCTTGGTGTCCTGTTGTC
TCCAGAACCCCAATATGTCTTGTGAAGGGCACACAACCCCTCAAAGGGGTGTCACTTCTT
CTGATCACTTTTGTTACTGTTTACTAACTGATCCTATGAATCACTGTGTCTTCTCAGAGG
CCGTGAACCACGTCTGCAAT

FASTA格式的簡(jiǎn)單性和靈活性卻帶來了令人遺憾的缺點(diǎn)：FASTA格式是一個(gè)松散定義的特殊格式（雖然這種現(xiàn)象在生物信息學(xué)中十分常見）睬辐。因此挠阁，你可能會(huì)遇到FASTA格式的各種變體，這些變體會(huì)造成微小的錯(cuò)誤溯饵，除非你的程序足夠強(qiáng)大侵俗。這就是為什么最好使用現(xiàn)有的FASTA/FASTQ解析庫而不是執(zhí)行自定義的解析庫；現(xiàn)有的庫已經(jīng)過開源社區(qū)審查（后續(xù)將更多討論這一部分內(nèi)容）丰刊。

關(guān)于FASTA格式最令人頭疼的是描述行中標(biāo)識(shí)符的格式?jīng)]有通用的規(guī)范隘谣。例如，如下FASTA描述行是否指代同一條目啄巧？

>ENSMUSG00000020122|ENSMUST00000138518
> ENSMUSG00000020122|ENSMUST00000125984
>ENSMUSG00000020122|ENSMUST00000125984|epidermal growth factor receptor
>ENSMUSG00000020122|ENSMUST00000125984|Egfr
>ENSMUSG00000020122|ENSMUST00000125984|11|ENSFM00410000138465

沒有標(biāo)識(shí)符的標(biāo)準(zhǔn)方案寻歧，我們不能使用簡(jiǎn)單的精確匹配確定一個(gè)標(biāo)識(shí)符是否和FASTA條目的標(biāo)題行匹配。相反秩仆，我們需要依靠FASTA文件描述行和標(biāo)識(shí)符之間的模糊匹配码泛。這可能會(huì)使得判斷標(biāo)準(zhǔn)變得相當(dāng)凌亂：匹配模式的寬容度應(yīng)該是多少？如果正則表達(dá)式過于寬泛是否會(huì)匹配到錯(cuò)誤的序列澄耍？基本上來說噪珊，模糊匹配是一個(gè)脆弱的策略。

幸運(yùn)的是齐莲，這個(gè)問題有更好的解決方案（也很簡(jiǎn)單）：與其事后依靠模糊匹配來糾正不一致的命名痢站，還不如一開始就擁有嚴(yán)格的命名規(guī)范并始終保持一致。然后选酗，在運(yùn)行任何外部來源的數(shù)據(jù)時(shí)哪痰，需要通過幾個(gè)合理性檢查，以確保數(shù)據(jù)遵循對(duì)應(yīng)的格式声滥。 這些檢查不需要太復(fù)雜（檢查重復(fù)的名稱，手動(dòng)檢查一些條目嫩舟，檢查 > 和標(biāo)識(shí)符之間錯(cuò)誤的空格，檢查不同文件之間名字的重合等）怀偷。

如果你需要整理外部數(shù)據(jù)家厌，請(qǐng)始終保留原始文件并編寫腳本將更正的版本寫入新的文件。這樣椎工，腳本可以輕松地重新運(yùn)行饭于，用于處理你獲得的新版本的原始數(shù)據(jù)集（但你仍然需要檢查每一項(xiàng)內(nèi)容－不要盲目信任數(shù)據(jù)！）维蒙。

常見的命名規(guī)范是將描述行利用第一個(gè)空格分為兩部分：標(biāo)識(shí)符和注釋掰吕。此格式的序列通常如下所示：

> gene_00284728 length = 231; type = dna 
GAGAACTGATTCTGTTACCGCAGGGCATTCGGATGTGCTAAGGTAGTAATCCATTATAAGTAACATGCGCGGAATATCCG 
GAGGTCATAGTCGTAATGCATAATTATTCCCTCCCTCAGAAGGACTCCCTTGCGAGACGCCAATACCAAAGACTTTCGTA 
GCTGGAACGATTGGACGGCCCAACCGGGGGGAGTCGGCTATACGTCTGATTGCTACGCCTGGACTTCTCTT 

這里 gene_00284728 是標(biāo)識(shí)符，length = 231; type = dna 是注釋颅痊。此外殖熟，ID應(yīng)該是唯一的。盡管不是標(biāo)準(zhǔn)規(guī)范斑响，但是將第一個(gè)空格之前的內(nèi)容作為標(biāo)識(shí)符菱属，之后的內(nèi)容作為非必要注釋在生物信息學(xué)程序中是很常見的（例如，BEDtools舰罚，Samtools和BWA都是這樣處理的）纽门。有了這個(gè)規(guī)范以后，通過標(biāo)識(shí)符查找一個(gè)特定的序列會(huì)很容易营罢，我們將在本章末尾討論如何高效地利用這種方式處理經(jīng)過索引的FASTA文件赏陵。

FASTQ格式

FASTQ格式通過給序列中每個(gè)堿基添加數(shù)字質(zhì)量分?jǐn)?shù)的形式擴(kuò)展了FASTA文件。FASTQ格式被廣泛用于存儲(chǔ)高通量測(cè)序數(shù)據(jù)饲漾，該數(shù)據(jù)利用每個(gè)堿基的質(zhì)量分?jǐn)?shù)來表示每個(gè)堿基測(cè)序的可信度蝙搔。然而和FASTA格式一樣，F(xiàn)ASTQ格式也擁有變體和缺陷考传，會(huì)使看似簡(jiǎn)單的格式變得很難處理吃型。

FASTQ格式如下所示：

@ DJB775P1:248:D0MDGACXX:7:1202:12362:49613 ①
TGCTTACTCTGCGTTGATACCACTGCTTAGATCGGAAGAGCACACGTCTGAA ②
+ ③
JJJJJIIJJJJJJHIHHHGHFFFFFFCEEEEEDBD?DDDDDDBDDDABDDCA ④
@ DJB775P1:248:D0MDGACXX:7:1202:12782:49716 
CTCTGCGTTGATACCACTGCTTACTCTGCGTTGATACCACTGCTTAGATCGG 
+ 
IIIIIIIIIIIIIIIHHHHHHFFFFFFEECCCCBCECCCCCCCCCCCCCCCC 

① 描述行，以@符號(hào)開始伙菊。包含標(biāo)識(shí)符記錄和其他信息败玉。
② 序列數(shù)據(jù)，可以是一行或多行镜硕。
③ 以 + 開頭的行运翼，位于序列數(shù)據(jù)之后，表示序列的結(jié)束兴枯。在舊版本的FASTQ文件中血淌，通常會(huì)這里重復(fù)描述行，但這是多余的，會(huì)導(dǎo)致不必要的FASTQ文件過大悠夯。
④ 質(zhì)量數(shù)據(jù)癌淮，也可以是一行或多行，但必須和序列長(zhǎng)度相同沦补。每個(gè)堿基質(zhì)量數(shù)值以ASCII字符編碼乳蓄，編碼規(guī)則我們將在后文討論（“堿基質(zhì)量”見344頁）。

與FASTA文件一樣夕膀，F(xiàn)ASTQ文件中一個(gè)常見的規(guī)范是將描述行利用第一個(gè)空格分解為兩部分：標(biāo)識(shí)符記錄和注釋虚倒。

如何正確地解析FASTQ文件其實(shí)是非常棘手的。一個(gè)常見的誤區(qū)是將每一個(gè)以 @ 符號(hào)開始的行都當(dāng)作說明行产舞。然而魂奥，@ 本身也是一個(gè)有效的質(zhì)量符號(hào)。FASTQ的序列行和質(zhì)量行可以延續(xù)到下一行易猫，所以以 @ 符號(hào)開頭的行可能是質(zhì)量行耻煤，而不是標(biāo)題行。所以准颓，寫一個(gè)始終把以 @ 開頭的行作為標(biāo)題行的解析程序會(huì)導(dǎo)致不正確的解析哈蝇。但是，我們可以利用質(zhì)量分?jǐn)?shù)的字符數(shù)必須等于序列的字符數(shù)這一準(zhǔn)則來準(zhǔn)確解析這種格式瞬场，這其實(shí)就是后文中 readfq 解析器的工作原理买鸽。

計(jì)算FASTA/FASTQ條目的輸入和輸出

作為純文本格式涧郊，可以輕松使用Unix工具處理FASTQ和FASTA文件贯被。一個(gè)常見的生物信息學(xué)命令行是這樣的：

$ grep -c "^>" egfr_flank.fasta 
5

正如本書47頁“Pipes實(shí)戰(zhàn)：利用grep和pipes創(chuàng)建簡(jiǎn)單的程序”所說，你必須引用** > **字符來防止shell把它當(dāng)作一個(gè)重定向運(yùn)算符（并覆蓋你的FASTA文件Ｗ彼摇）彤灶。這是一個(gè)保險(xiǎn)的計(jì)算FASTA文件數(shù)目的方式，因?yàn)殡m然格式定義寬松批旺，但是每個(gè)序列都具有一個(gè)描述行幌陕，只有這些描述行以 > 符號(hào)開始 。

我們可能會(huì)嘗試使用類似的方法處理FASTQ文件汽煮，用 @ 符號(hào)代替 >：

$ grep -c "^@" untreated1_chr4.fq 
208779

這意味著 untreated1_chr4.fq 有208779個(gè)條目搏熄。但是，仔細(xì)檢查 untreated1_chr4.fq 文件暇赤，你會(huì)發(fā)現(xiàn)每個(gè)FASTQ條目占用四行心例，但總行數(shù)是：

$ wc -l untreated1_chr4.fq 
817420 untreated1_chr4.fq 

而817420/4 = 204355，這和 grep -c 命令得到的結(jié)果有很大不同鞋囊！這是發(fā)生了什么呢止后？請(qǐng)記住，@ 是一個(gè)有效的質(zhì)量符號(hào)，質(zhì)量行可以以這個(gè)符號(hào)開始译株，可以使用 grep "^@" untreated1_chr4.fq | less 命令行來查看這樣的例子瓜喇。

如果你很確定你的FASTQ文件每個(gè)序列條目占用四行，可通過 wc -l 命令估計(jì)文件行數(shù)并除以4的方式來估計(jì)序列條目歉糜。如果你不確定某些FASTQ條目是否包含多行乘寒，一個(gè)更強(qiáng)大的計(jì)算序列條目的方法是bioawk：

$ bioawk -cfastx 'END {print NR}' untreated1_chr4.fq 
204355 

核苷酸編碼

隨著基礎(chǔ)FASTA / FASTQ格式的覆蓋，讓我們來看看以這些格式出現(xiàn)的核苷酸和堿基質(zhì)量分?jǐn)?shù)的標(biāo)準(zhǔn)編碼匪补。 顯然肃续，編碼核苷酸很簡(jiǎn)單：A，T叉袍，C始锚，G分別代表核苷酸腺嘌呤，胸腺嘧啶喳逛，胞嘧啶和鳥嘌呤瞧捌。小寫的堿基經(jīng)常被用來代表“軟覆蓋”重復(fù)序列或低復(fù)雜度序列（RepeatMasker 和 Tandem Repeats Finder 程序均采用這種格式）。重復(fù)序列和低復(fù)雜度的序列也可以被“硬覆蓋”润文，核苷酸序列被替換為N（有時(shí)是X）姐呐。

簡(jiǎn)并（或模糊）核苷酸編碼被用于表示兩個(gè)或多個(gè)堿基。例如典蝌，N可以用于表示任何堿基曙砂。國際純粹與應(yīng)用化學(xué)聯(lián)合會(huì)（IUPAC）擁有一組標(biāo)準(zhǔn)化的核苷酸編碼，既包括明確的也包括簡(jiǎn)并的（見表10-1）骏掀。

表 10-1 IUPAC 核苷酸代碼

有些生物信息學(xué)程序?qū)τ诤?jiǎn)并核苷酸的處理可能會(huì)有所不同。例如葵袭，BWA讀長(zhǎng)比對(duì)工具可以將參考基因組中的簡(jiǎn)并核苷酸字符轉(zhuǎn)換為隨機(jī)堿基（Li and Durbin, 2009）涵妥，但是由于使用了隨機(jī)種子集合，所以不會(huì)在重新索引比對(duì)數(shù)據(jù)的時(shí)候產(chǎn)生兩種不同的結(jié)果坡锡。

堿基質(zhì)量

FASTQ條目的每個(gè)序列堿基在質(zhì)量行都具有相應(yīng)的數(shù)字質(zhì)量分?jǐn)?shù)蓬网。每個(gè)堿基的質(zhì)量分?jǐn)?shù)都編碼為一個(gè)ASCII字符。質(zhì)量行看起來就像一串隨機(jī)字符鹉勒，如第四行所示：

@ AZ1:233:B390NACCC:2:1203:7689:2153 
GTTGTTCTTGATGAGCCATGAGGAAGGCATGCCAAATTAAAATACTGGTGCGAATTTAAT 
+ 
CCFFFFHHHHHJJJJJEIFJIJIJJJIJIJJJJCDGHIIIGIGIJIJIIIIJIJJIJIIH 

（這個(gè)FASTQ條目在本章的 README 文件中帆锋，如果你想按照這里所說的跟著一起看看。）

記住贸弥，ASCII字符僅在內(nèi)部表示為0到127之間的整數(shù)（詳情參見 man ascii 命令）窟坐。因?yàn)椴⒉皇撬械腁SCII字符都可以在屏幕上顯示出來（例如，字符呼應(yīng) "\ 07" 會(huì)發(fā)出“叮”的聲音）哲鸳，質(zhì)量分?jǐn)?shù)僅限于可打印的ASCII字符臣疑，范圍從33到126（代表空格的字符32被省略）。

所有編程語言都具有將字符轉(zhuǎn)換為十進(jìn)制ASCII碼和將十進(jìn)制ASCII碼轉(zhuǎn)換為字符的功能徙菠。在Python中讯沈，分別對(duì)應(yīng)函數(shù) ord() 和 chr()。讓我們使用 ord() 函數(shù)在Python交互式解釋器中將質(zhì)量字符轉(zhuǎn)換為十進(jìn)制ASCII碼列表：

>>> qual = "JJJJJJJJJJJJGJJJJJIIJJJJJIGJJJJJIJJJJJJJIJIJJJJHHHHHFFFDFCCC" 
>>> [ord(b) for b in qual] 
[74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 71, 74, 74, 74, 74, 74, 73,  
73, 74, 74, 74, 74, 74, 73, 71, 74, 74, 74, 74, 74, 73, 74, 74, 74, 74, 74,  
74, 74, 73, 74, 73, 74, 74, 74, 74, 72, 72, 72, 
67, 67, 67] 

遺憾的是婿奔，將這些ASCII數(shù)值轉(zhuǎn)換成有意義的質(zhì)量分?jǐn)?shù)可能非常棘手缺狠，因?yàn)槟壳按嬖谌N不同的質(zhì)量體系：Sanger，Solexa和Illumina（見表10-2）萍摊。負(fù)責(zé)Biopython挤茄，BioPerl和BioRuby這類項(xiàng)目的開放性生物信息學(xué)基金會(huì)（OBF），將它們命名為 fastq-sanger冰木，fastq-solexa 和 fastaq-illumina穷劈。值得慶幸的是，生物信息學(xué)領(lǐng)域最終似乎確定統(tǒng)一使用Sanger編碼（也就是這里顯示的質(zhì)量行格式）踊沸，因此我們將使用此方案逐步完成轉(zhuǎn)換過程歇终。

表 10-2 FASTQ質(zhì)量體系（來源 Cock et al., 2010，已獲得許可）

首先逼龟，我們需要減去一個(gè)偏移值來將這個(gè)Sanger質(zhì)量分?jǐn)?shù)轉(zhuǎn)換為PHRED質(zhì)量分?jǐn)?shù)评凝。PHRED是由Phil Green撰寫的早期堿基調(diào)用程序，用于解析他自己寫的熒光追蹤數(shù)據(jù)腺律。通過查看表10-2奕短，我們注意到Sanger格式的偏移量為33，因此我們從每個(gè)質(zhì)量分?jǐn)?shù)中減去33：

>>> phred = [ord(b)-33 for b in qual]
>>> phred
[41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 38, 41, 41, 41, 41, 41, 40,
40, 41, 41, 41, 41, 41, 40, 38, 41, 41, 41, 41, 41, 40, 41, 41, 41, 41, 41,
41, 41, 40, 41, 40, 41, 41, 41, 41, 39, 39, 39, 39, 39, 37, 37, 37, 35, 37,
34, 34, 34]

將概率轉(zhuǎn)化為質(zhì)量崖飘，我們使用這個(gè)函數(shù)的倒數(shù)：

在我們的例子中榴捡，我們需要前一個(gè)方程式。將其應(yīng)用于我們的PHRED質(zhì)量分?jǐn)?shù)：

>>> [10**(-q/10) for q in phred]
[1e-05, 1e-05, 1e-05, 1e-05, 1e-05, 1e-05, 1e-05, 1e-05, 1e-05, 1e-05, 1e-05,
1e-05, 0.0001, 1e-05, 1e-05, 1e-05, 1e-05, 1e-05, 0.0001, 0.0001, 1e-05,
1e-05, 1e-05, 1e-05, 1e-05, 0.0001, 0.0001, 1e-05, 1e-05, 1e-05, 1e-05,
1e-05, 0.0001, 1e-05, 1e-05, 1e-05, 1e-05, 1e-05, 1e-05, 1e-05, 0.0001,
1e-05, 0.0001, 1e-05, 1e-05, 1e-05, 1e-05, 0.0001, 0.0001, 0.0001, 0.0001,
0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001]

公式。關(guān)于這個(gè)格式的更多詳情梁沧，請(qǐng)參閱Cock et al., 2010檀何。

示例：檢查和去除低質(zhì)量堿基

請(qǐng)注意前一個(gè)例子中堿基的準(zhǔn)確度是如何下降的；這是Illumina測(cè)序的特征誤差分布。本質(zhì)上來說频鉴，我們通過這種測(cè)序技術(shù)獲得的讀長(zhǎng)的堿基測(cè)序錯(cuò)誤概率逐漸增加（朝向3'端）栓辜。這可能對(duì)下游分析產(chǎn)生深遠(yuǎn)的影響！當(dāng)處理測(cè)序數(shù)據(jù)時(shí)垛孔，你應(yīng)該總是

• 了解測(cè)序技術(shù)的錯(cuò)誤分布和局限性（例如藕甩，是否受到GC含量的影響）
• 考慮上述內(nèi)容將如何影響你的分析

所有這一切都是實(shí)驗(yàn)特異性的，而且需要精心的規(guī)劃周荐。

把我們關(guān)于堿基準(zhǔn)確度的Python列表作為一個(gè)學(xué)習(xí)工具是很用的狭莱，可以看到如何將概率轉(zhuǎn)換為質(zhì)量，但它不能幫助我們了解數(shù)百萬序列的質(zhì)量概況概作。從這個(gè)意義上來說腋妙，一圖抵千字－有軟件可以幫助我們看到讀長(zhǎng)中堿基的質(zhì)量分布。最流行的是Java程序 FastQC讯榕，這款軟件很容易運(yùn)行并輸出有用的圖形和質(zhì)量指標(biāo)辉阶。如果你喜歡使用R語言，你可以使用Bioconductor中的軟件包 qrqc（由鄙人所寫）瘩扼。我們將在示例中使用 qrqc 軟件包谆甜，這樣我們就可以自己了解如何可視化這些數(shù)據(jù)。

我們先來安裝運(yùn)行這個(gè)例子所有必要的程序集绰。首先规辱，利用如下命令行在R中安裝 qrqc 軟件包：

> library(BiocInstaller)
> biocLite('qrqc')

接下來，讓我們來安裝兩個(gè)程序用來去除低質(zhì)量堿基：sickle 和 seqtk 栽燕。seqtk 是由Heng Li編寫的一個(gè)通用序列工具包 - 擁有在序列末端去除低質(zhì)量堿基的子命令（還有許多有用的其他功能）罕袋。利用Homebrew軟件包管理系統(tǒng)，sickle 和 seqtk 很容易在Mac OS X操作系統(tǒng)上安裝（例如碍岔，利用brew install seqtk 和 brew install sickle 命令行）浴讯。

在安裝好這些程序之后，讓我們來開始修整在GitHub存儲(chǔ)庫中本章目錄下的 untreated1_chr4.fq FASTQ文件蔼啦。這個(gè)FASTQ文件由Bioconductor中 pasillaBamSubset 軟件包的 untreated1_chr4.bam BAM文件生成（更多信息請(qǐng)參見本章目錄中的 README 文件）榆纽。為了保持事情簡(jiǎn)單，我們將使用每個(gè)程序的默認(rèn)設(shè)置捏肢。從 sickle 軟件包開始：

$ sickle se -f untreated1_chr4.fq -t sanger -o untreated1_chr4_sickle.fq
FastQ records kept: 202866
FastQ records discarded: 1489

sickle 通過 -f 指定輸入文件奈籽，通過 -t 指定質(zhì)量分?jǐn)?shù)類型，通過 -o 得到修整后的文件鸵赫。

現(xiàn)在衣屏，讓我們來運(yùn)行 seqtk trimfq 命令，這個(gè)命令需要一個(gè)參數(shù)并以標(biāo)準(zhǔn)形式輸出修整后的序列：

$ seqtk trimfq untreated1_chr4.fq > untreated1_chr4_trimfq.fq 

讓我們?cè)赗中比較上述結(jié)果辩棒。我們使用 qrqc 軟件包根據(jù)位置收集這些文件中的堿基質(zhì)量分布狼忱，然后使用 ggplot2 軟件包可視化這些結(jié)果膨疏。我們可以逐一加載這些軟件包，但一個(gè)不錯(cuò)的工作流程是利用 lapply() 函數(shù)自動(dòng)執(zhí)行：

# trim_qual.R -- 檢查修剪前后堿基質(zhì)量
library(qrqc)

# FASTQ文件
fqfiles <- c(none="untreated1_chr4.fq", 
             sickle="untreated1_chr4_sickle.fq", 
             trimfq="untreated1_chr4_trimfq.fq")

# 使用 qrqc 軟件包的 readSeqFile 函數(shù)加載每個(gè)文件
# 我們只需要堿基質(zhì)量钻弄，所以關(guān)掉 
# readSeqFile 函數(shù)的一些其他功能成肘。
seq_info <- lapply(fqfiles, function(file) {
                   readSeqFile(file, hash=FALSE, kmer=FALSE)
                   }) 

# 將堿基質(zhì)量提取為數(shù)據(jù)框，并追加 
# 一列使用修剪工具（或沒有使用）的情況斧蜕。這部分
# 在后面的作圖中使用双霍。
quals <- mapply(function(sfq, name) {
                                     qs <- getQual(sfq)
                                     qs$trimmer <- name
                                     qs
                                    }, seq_info, names(fqfiles), SIMPLIFY=FALSE)


# 將列表中單獨(dú)的數(shù)據(jù)框合并為單個(gè)數(shù)據(jù)框
d <- do.call(rbind, quals) 

# 堿基質(zhì)量可視化
p1 <- ggplot(d) + geom_line(aes(x=position, y=mean, linetype=trimmer))
p1 <- p1 + ylab("mean quality (sanger)") + theme_bw()
print(p1) 

# 使用 qrqc 軟件包的 qualPlot 函數(shù)與列表生成面板圖
# 只顯示 10% 到 90% 的分位數(shù)和低通曲線
p2 <- qualPlot(seq_info, quartile.color=NULL, mean.color=NULL) + theme_bw()
p2 <- p2 + scale_y_continuous("quality (sanger)")
print(p2) 

圖 10-1 修剪前以及利用 sickle 和 seqtk trimfq 處理后的讀長(zhǎng)平均堿基質(zhì)量

圖 10-2 修剪前以及利用 sickle 和 seqtk trimfq 處理后每個(gè)堿基位置的 10% 到 90% 分位數(shù)和低通曲線

在一行序列中，我們可以從末尾去除低質(zhì)量的堿基－修剪命令并不困難掀宋。更重要的步驟是通過比較去除低質(zhì)量堿基之前和之后的文件來可視化這些修剪程序?qū)ξ覀兊臄?shù)據(jù)做了怎樣的修改深纲。檢查程序如何更改我們的數(shù)據(jù)而不是盲目相信它們才是強(qiáng)大的生物信息學(xué)的重要組成部分，也是遵守黃金法則（不要信任你的工具）的表現(xiàn)劲妙。在這個(gè)例子中湃鹊，檢查一小部分?jǐn)?shù)據(jù)只需要不超過20行代碼（忽略空白行和提高可讀性的注釋）和幾分鐘的時(shí)間，但關(guān)于這些程序?qū)ξ覀兊臄?shù)據(jù)做了什么和程序之間的差異镣奋，這些步驟給了我們有價(jià)值的理解币呵。如果我們?cè)敢猓覀円部梢愿牟煌脑O(shè)置來同時(shí)運(yùn)行這兩個(gè)質(zhì)量修飾程序侨颈，并比較這些設(shè)置對(duì)結(jié)果的影響余赢。大部分嚴(yán)謹(jǐn)?shù)纳镄畔W(xué)都是這樣的過程：運(yùn)行程序，對(duì)比原始數(shù)據(jù)與輸出結(jié)果哈垢，運(yùn)行程序妻柒，比較輸出結(jié)果，如此往復(fù)耘分。

FASTA/FASTQ解析示例：核苷酸計(jì)數(shù)

編寫你自己的FASTA/FASTQ解析器其實(shí)并不是很困難举塔，事實(shí)上，這是一個(gè)非常有教育性的編程練習(xí)陶贼。 但是當(dāng)我們需要使用解析器來處理真實(shí)的數(shù)據(jù)時(shí)啤贩，最好還是使用可靠的現(xiàn)有庫。記住本章開頭引用的話：杰出的程序員知道何時(shí)重用代碼拜秧。有很多開源的，免費(fèi)提供的章郁，可靠的解析器已經(jīng)經(jīng)過了開源社區(qū)審核枉氮，并且旨在加速解析過程志衍。我們將使用Heng Li的 readfq 應(yīng)用程序，因?yàn)樗梢越馕鯢ASTA和FASTQ文件聊替，使用起來很簡(jiǎn)單楼肪，是獨(dú)立應(yīng)用（意味著它不需要安裝任何依賴項(xiàng)）。Biopython和BioPerl這兩個(gè)流行的程序庫中有很好的可選FASTA/FASTQ解析器惹悄。

我們將使用Python版本的 readfq 程序春叫，readfq.py。你可以從GitHub上下載這個(gè)Python文件（本書的GitHub資源庫中也有這個(gè)文件的副本）泣港。雖然我們可以使用 from readfq import readfq 命令來執(zhí)行 readfq.py 中的FASTA/FASTQ解析程序暂殖，但是對(duì)于單文件腳本其實(shí)更簡(jiǎn)單，你可以直接將程序復(fù)制粘貼到你的腳本中当纱。在本書中呛每，我們使用 from readfq import readfq 命令來避免占據(jù)例子文件的空間。

readfq 程序的 readfq() 函數(shù)使用簡(jiǎn)單坡氯。readfq() 采用一個(gè)文件對(duì)象（例如晨横，一個(gè)已經(jīng)通過
open('filename.fa') 或 sys.stdin 打開的文件名稱），并會(huì)產(chǎn)生FASTQ/FASTA條目箫柳，直到它到達(dá)文件或流的末尾手形。每個(gè)FASTQ/FASTA條目由 readfq() 函數(shù)返回為包含條目描述，序列和質(zhì)量的元組悯恍。如果 readfq() 正在讀取一個(gè)FASTA文件叁幢，質(zhì)量行是無。這就是利用 readfq() 函數(shù)讀取FASTQ/FASTA序列行的所有內(nèi)容坪稽。

我們編寫一個(gè)簡(jiǎn)單的程序來計(jì)算文件中每個(gè)IUPAC核苷酸的數(shù)量：

#!/usr/bin/env python 
# nuccount.py -- 一個(gè)文件中的所有核苷酸數(shù)目 
import sys
from collections import Counter ①
from readfq import readfq

IUPAC_BASES = "ACGTRYSWKMBDHVN-." ②

# 初始化計(jì)數(shù)器
counts = Counter() ③

for name, seq, qual in readfq(sys.stdin): ④
    # 對(duì)于每一個(gè)序列條目顷帖，向計(jì)數(shù)器添加所有的堿基數(shù)量
    counts.update(seq.upper()) ⑤

＃打印結(jié)果 
for base in IUPAC_BASES:
print base + "\t" + str(counts[base]) ⑥

① 我們使用 collections 模塊中的 Counter 計(jì)數(shù)核苷酸。Counter 的工作原理很像Python的 dict 類（技術(shù)上來說渤滞，它是 dict 的一個(gè)子類）贬墩，具有使計(jì)數(shù)更容易的附加功能。</br>
② 這個(gè)全局變量定義了所有的IUPAC核苷酸妄呕，我們稍后將用它來打印順序一致的堿基（因?yàn)橄馪ython字典一樣陶舞，Counter 對(duì)象不保持順序）。在腳本的頂部放置類似IUPAC_BASES的大寫常量讓讀者清楚是一個(gè)很好的做法绪励。</br>
③ 創(chuàng)建一個(gè)新的 Counter 對(duì)象肿孵。</br>
④ 此行使用 readfq 模塊中的 readfq 函數(shù)從文件句柄參數(shù)讀取FASTQ/FASTQ條目（在這種情況下唠粥，sys.stdin，標(biāo)準(zhǔn)輸入）轉(zhuǎn)換成 name停做，seq 和 qual 變量（通過Python中被稱為元組拆封的功能）晤愧。</br>
⑤ Counter.update() 方法讀取任何可迭代的對(duì)象（在此情況下，序列堿基的字符串）蛉腌，并將它們添加到計(jì)數(shù)器中官份。我們也可以使用for循環(huán)遍歷 seq 變量的每個(gè)字符，利用 counts[seq.upper()] 命令遞增計(jì)數(shù)烙丛。請(qǐng)注意舅巷，我們用 upper() 函數(shù)把所有的字符轉(zhuǎn)換為大寫，因此小寫的“軟覆蓋”堿基也會(huì)被計(jì)數(shù)蜀变。</br>
⑥ 最后悄谐，我們迭代所有的IUPAC堿基并輸出它們的數(shù)量。 </br>

這個(gè)版本需要通過標(biāo)準(zhǔn)輸入獲得輸入文件库北，所以保存此文件并使用 chmod + X nuccount.py 命令添加執(zhí)行權(quán)限后爬舰，我們可以通過下述命令運(yùn)行：

$ cat contam.fastq | ./nuccount.py 
A   103 
C   110 
G   94 
T   109 
R   0 
Y   0 
S   0 
W   0 
K   0 
M   0 
B   0 
D   0 
H   0 
V   0 
N   0 
-   0 
.   0 

請(qǐng)注意，我們不一定要讓Python腳本可執(zhí)行寒瓦；另一個(gè)選擇是簡(jiǎn)單地用 python nuccount.py 命令啟動(dòng)它 情屹。 無論哪種方式都會(huì)得到相同的結(jié)果，所以這最終是一種風(fēng)格的選擇杂腰。但是垃你，如果你確實(shí)想讓它變成一個(gè)可執(zhí)行的腳本，記住要做以下的事情：

• 包含標(biāo)識(shí)行 #!/usr/bin/env python
• 使用 chmod + X <scriptname.py> 命令使腳本可執(zhí)行

這個(gè)腳本有很多可以改進(jìn)的地方：添加對(duì)每個(gè)序列堿基組成統(tǒng)計(jì)的支持喂很，采用文件參數(shù)而不是通過標(biāo)準(zhǔn)輸入惜颇，“軟覆蓋”堿基（小寫）計(jì)數(shù)，CpG位點(diǎn)計(jì)數(shù)或在文件中發(fā)現(xiàn)非IUPAC核苷酸時(shí)發(fā)出警告少辣。核苷酸計(jì)數(shù)很簡(jiǎn)單－腳本中最復(fù)雜的部分是 readfq() 函數(shù)凌摄。這是重用代碼之美：書寫良好的函數(shù)和庫讓我們不需要重寫復(fù)雜的解析器。相反漓帅，我們使用更廣泛的社區(qū)合作開發(fā)的軟件锨亏，測(cè)試，排除故障忙干，使其更有效率（參見 readfq 軟件在GitHub上的提交歷史作為例子）器予。重用軟件不是欺騙－這是行家寫程序的方式。

索引FASTA文件

我們經(jīng)常需要對(duì)FASTA文件的子序列進(jìn)行有效的隨機(jī)訪問（給定區(qū)域）捐迫。乍一看來乾翔，寫一個(gè)腳本來完成似乎并不困難。例如弓乙，我們可以編寫一個(gè)遍歷FASTA條目的腳本末融，提取與指定區(qū)域重疊的序列钧惧。但是暇韧，在提取一些隨機(jī)序列時(shí)勾习，這不是一個(gè)有效的方法。為了知道為什么懈玻，假設(shè)我們需要訪問小鼠基因組第8號(hào)染色體（123,407,082至123,410,742區(qū)域的）序列巧婶。這種方法會(huì)不必要地在內(nèi)存中解析和加載1號(hào)染色體至7號(hào)染色體的數(shù)據(jù)，即使我們不需要從這些染色體中提取子序列涂乌。從磁盤中讀取全部染色體并將其復(fù)制到內(nèi)存中可能效率相當(dāng)?shù)拖拢覀優(yōu)榱颂崛?.6 kb數(shù)據(jù)必須加載所有125 Mb的8號(hào)染色體艺栈！從FASTA文件中提取大量隨機(jī)子序列可能計(jì)算成本相當(dāng)高昂。

允許容易和快速隨機(jī)訪問的常見計(jì)算策略是索引文件湾盒。索引文件在生物信息學(xué)中是普遍存在的湿右；在未來的章節(jié)中，我們將索引基因組罚勾，以便它們可以用作比對(duì)參考序列毅人，我們將索引包含比對(duì)讀長(zhǎng)的文件，用于快速隨機(jī)訪問尖殃。在本節(jié)中丈莺，我們將介紹一下索引后的FASTA文件是如何允許我們快速輕松地提取子序列的。一般來說送丰，我們經(jīng)常索引整個(gè)基因組缔俄，但為了簡(jiǎn)化本章剩余部分的例子，我們只考慮小鼠基因組的第8號(hào)染色體器躏。從本章的GitHub的文件目錄中下載 Mus_musculus.GRCm38.75.dna.chromosome.8.fa.gz 文件并用 gunzip 命令解壓俐载。

索引壓縮的FASTA文件

雖然 samtools faidx 命令確實(shí)可以處理壓縮的FASTA文件，但是它只適用于BGZF壓縮文件（更深入的
話題登失，我們將在425頁“快速訪問經(jīng)過BGZF和Tabix命令索引的制表符分隔的序列文件”進(jìn)一步討論）遏佣。

我們將使用Samtools對(duì)該文件進(jìn)行索引，它是一種廣泛使用的操作SAM和BAM比對(duì)格式的工具套件（我們將在第11章介紹更多細(xì)節(jié)）壁畸。你可以通過一個(gè)端口或軟件包管理器安裝 samtools （例如贼急，Mac OS X系統(tǒng)的Homebrew或Ubuntu系統(tǒng)的 apt-get）。 很像Git捏萍，Samtools使用子命令來完成的任務(wù)太抓。首先，我們需要使用 faidx 子命令索引FASTA文件：

$ samtools faidx Mus_musculus.GRCm38.75.dna.chromosome.8.fa 

這個(gè)命令將創(chuàng)建一個(gè)名為 Mus_musculus.GRCm38.75.dna.chromosome.8.fa.fai 的索引文件令杈。我們?cè)谔崛∽有蛄袝r(shí)不必?fù)?dān)心這個(gè)索引文件的細(xì)節(jié)特征－ samtools faidx 命令會(huì)注意這些細(xì)節(jié)的走敌。要訪問特定區(qū)域的子序列，我們使用 samtools faidx <in.fa> <region> 命令逗噩，其中 <in.fa> 是FASTA文件（你剛剛索引的）掉丽，<region> 的格式為 chromosome: start-end（染色體：起始位置－終止位置）跌榔。例如：

$ samtools faidx Mus_musculus.GRCm38.75.dna.chromosome.8.fa 8:123407082-1??23410744 
>8:123407082-1??23410744 
GAGAAAAGCTCCCTTCTTCTCCAGAGTCCCGTCTACCCTGGCTTGGCGAGGGAAAGGAAC 
CAGACATATATCAGAGGCAAGTAACCAAGAAGTCTGGAGGTGTTGAGTTTAGGCATGTCT 
[...] 

一定要注意染色體語法的差異（UCSC的Chr 8與ENSEMBL的8的格式）。如果 samtools faidx 命令沒有序列返回捶障，這可能就是原因僧须。

samtools faidx 命令允許一次訪問多個(gè)區(qū)域，因此我們可以這樣做：

$ samtools faidx Mus_musculus.GRCm38.75.dna.chromosome.8.fa  \
     8:123407082-1??23410744 8:123518835-123536649 
>8:123407082-1??23410744 
GAGAAAAGCTCCCTTCTTCTCCAGAGTCCCGTCTACCCTGGCTTGGCGAGGGAAAGGAAC 
CAGACATATATCAGAGGCAAGTAACCAAGAAGTCTGGAGGTGTTGAGTTTAGGCATGTCT 
[...] 
>8:123518835-123536649 
TCTCGCGAGGATTTGAGAACCAGCACGGGATCTAGTCGGAGTTGCCAGGAGACCGCGCAG 
CCTCCTCTGACCAGCGCCCATCCCGGATTAGTGGAAGTGCTGGACTGCTGGCACCATGGT 
[...] 

什么讓訪問索引文件更快项炼？

在第3章（見45頁“全能UNIX流程：速度與美貌共存”）担平， 我們討論了為什么從磁盤讀取和寫入數(shù)據(jù)相比存儲(chǔ)在內(nèi)存中的數(shù)據(jù)非常緩慢。我們可以通過使用指向文件中某些區(qū)塊位置的索引避免不必要地從磁盤中讀取整個(gè)文件锭部。在我們使用的FASTA文件的情況下暂论，索引基本上存儲(chǔ)量每個(gè)序列在文件中的起始位置（以及其他必要的信息）。

當(dāng)我們查找像8號(hào)染色體（123,407,082-123,410,744）這樣的區(qū)域時(shí)拌禾，samtools faidx 命令使用索引文件中的信息快速計(jì)算那些堿基在文件中的具體位置取胎。然后，使用一個(gè)稱為文件搜索的操作湃窍，程序跳轉(zhuǎn)到文件中的精確位置（稱為平移）闻蛀，并開始讀取序列。擁有預(yù)先計(jì)算的文件偏移量和跳轉(zhuǎn)到這些確切位置的能力是快速訪問索引文件特定區(qū)域的原因坝咐。