測序的世界很奇妙岖常，不同的數(shù)據(jù)處理可能得出不同的結(jié)論歧寺，入門生信首先要做的就是了解你的數(shù)據(jù)
還等什么翅雏？跟我一起來探索吧～

測序原理：

早期如何測序呢纵散？

如果給你這么一串珠子梳码，讓你統(tǒng)計綠色，白色伍掀，紅色掰茶，其他這四種類型的數(shù)量，那你是不是要從某個位置一個一個數(shù)蜜笤？沒錯濒蒋，你已經(jīng)開始對這串珠子進行了肉眼測序。

肉眼測序.png

其實對于DNA測序把兔，思路也是這樣沪伙。

早在1954年，Whitfeld等就提出了測定多聚核糖核苷酸鏈的降解法县好，該方法利用磷酸單酯酶的脫磷酸作用和高碘酸鹽的氧化作用從鏈末端逐一分離寡核糖核苷酸并測定其種類围橡。目的就是想通過這種一個一個“數(shù)”的方法來得到DNA的堿基順序。

這里再補充一個小知識缕贡，DNA有多大翁授？它的直徑也就2nm拣播，兩個核苷酸之間的小溝0.34nm。肉眼觀察肯定不行收擦，那么顯微鏡呢诫尽？光學(xué)顯微鏡可見光的波長約為700nm，單單一個直徑就看不到炬守，更別說看更小的分子了DNA對于它來說根本不可見牧嫉。要想用普通顯微鏡觀察不可能，更強大的電子顯微鏡呢减途？別忘了酣藻，DNA分子并不是平躺著讓你去觀察，他有四級結(jié)構(gòu)鳍置，不斷扭曲辽剧、折疊，簡單的二維圖片是很難區(qū)分開來的税产，能看到也僅僅只能看到一小部分怕轿。就像上面的這張珠子圖，拍一張照片辟拷，就能準(zhǔn)確區(qū)分被壓著的珠子和它上面的珠子順序嗎撞羽？不能夠！

人類基因組共有31.6億個DNA堿基對衫冻，23對染色體诀紊，測一條染色體就需要測1-2億個堿基對，像這種數(shù)數(shù)的辦法肯定行不通隅俘，太耗費人力邻奠。那么有沒有什么辦法可以實現(xiàn)半自動化呢？

視覺不行为居，另辟蹊徑碌宴，聽過“盲人摸象”的故事吧，70年代的Sanger發(fā)明了“雙脫氧終止反應(yīng)法”蒙畴，他就是利用了雙脫氧核苷酸 ddNTP去摸索DNA分子贰镣。例如，一條序列5‘ --> 3' 實際為TGACTTCG但我們之前不知道忍抽，這樣操作：

操作過程可能會有些無聊八孝，但是不難懂，希望你能理解其中的意思鸠项，因為后來測序都是受此啟發(fā)，這是鼻祖

• 設(shè)置四個反應(yīng)體系1-4子姜，分別加入引物祟绊、DNA聚合酶楼入、四種dNTP、一定比例的ddNTP（帶有放射性標(biāo)記）例如1中是ddATP牧抽，它就負(fù)責(zé)測定T堿基的位置嘉熊；依次2是ddCTP，3是ddTTP扬舒， 4是ddGTP

• 假如擴增過程中ddATP遇到了T位點阐肤，就結(jié)合并終止（因為ddNTP的2‘和3'都沒有羥基），那么其他的dNTP就無法結(jié)合讲坎。也許你會問孕惜，那么每次如果是這四種ddNTP占了先機先結(jié)合了，其他的都無法再結(jié)合晨炕，結(jié)果不就是一種可能衫画，只能確定四個位點？并不是這樣： 大批量的核苷酸結(jié)合本就是一個隨機問題瓮栗，不一定每次都是ddNTP結(jié)合削罩，有可能是dNTP結(jié)合了好幾個然后ddNTP才回過神，結(jié)合完再宣告一個聚合過程結(jié)束费奸。

  想象ddNTP是VIP會員弥激，dNTP為普通用戶，只要ddNTP不結(jié)合（不占用綠色通道）愿阐，普通的dNTP就有機會秆撮。

• 因此，在這反反復(fù)復(fù)的結(jié)合過程中换况，就會把所有的位點都結(jié)合（因為測序測的不是幾百幾千條序列职辨，而是幾千萬條序列），這其中有的是dNTP結(jié)合戈二，有的是ddNTP結(jié)合舒裤，但是，在一個大規(guī)模樣本中觉吭，ddNTP會結(jié)合所有位點腾供，只是位點結(jié)合次數(shù)多少的問題，可能這一個位點有五千條序列測到這列就停了鲜滩，但是下一個位點又有其他的8千條序列測到伴鳖。 說的簡單一點就是，一段時間內(nèi)大量的ddNTP會結(jié)合完所有測序位點徙硅。

• Q: 最后怎么看結(jié)果榜聂？這里很有意思，一個位點可能會被大量的dNTP和ddNTP結(jié)合嗓蘑，怎么去除dNTP的干擾

A: 其實也不用刻意去除须肆，最后利用凝膠電泳和放射自顯影只能看到帶有熒光標(biāo)記的ddNTP匿乃，他們的排列順序先利用電泳條帶前后關(guān)系確定下，再用A-T, T-A, C-G, G-C關(guān)系反轉(zhuǎn)一下豌汇，就能知道我們的測序序列

• 附上一張圖幢炸，幫助理解
  
  Sanger測序.JPEG
  
  有了早期的第一次測序成功，才有了后來1983年的Kary Mullis 發(fā)明PCR測序儀拒贱，利用PCR才有了我們更加效率的NGS（二代測序）宛徊。進步的是方法，不變的是基本理念逻澳。

當(dāng)然闸天，千萬別小瞧了一代測序（Sanger測序），到今天這種方法還是經(jīng)常使用赡盘。他之所以不是現(xiàn)在的主流号枕，不是因為它測不準(zhǔn)（恰恰相反，第一代測序技術(shù)的主要特點就是測序讀長可達(dá)1000bp陨享，準(zhǔn)確性高達(dá)99.999%葱淳，二三代所不能及），而是因為它的通量低抛姑，成本高赞厕。
要知道，目前二代測序的錯誤率還在
目前一代測序在驗證序列（就是平時送公司測序返回來自己blast的那些）以及驗證基因組組裝完整性方面都是金標(biāo)準(zhǔn)定硝。

二代測序：

一代測序的原理我們搞清楚了皿桑，沒有PCR，利用聚合酶延伸鏈蔬啡。這就有一個問題诲侮，酶的活性會下降，所以一代最長能測1000bp箱蟆，再多了就要重頭開始一遍導(dǎo)致成本高沟绪；另外它一次只測一條，也就是所謂的通量低空猜。但是呢绽慈，他的確有值得借鑒的地方，就是準(zhǔn)確度很高辈毯，99.999%坝疼。人類基因組計劃HGP，1990年-2003年完成谆沃，利用的就是改進的Sanger測序法

后來技術(shù)不斷改進钝凶，Roche公司的454技術(shù)、illumina公司的Solexa/Hiseq技術(shù)和ABI公司的SOLID技術(shù)標(biāo)志第二代測序技術(shù)誕生管毙。其中Roche公司的454測序系統(tǒng)是第二代測序技術(shù)中第一個商業(yè)化運營的測序平臺腿椎。

其中Illumina市場規(guī)模占到75%以上桌硫，主打Hiseq夭咬，下面??就主要介紹他的PE（Pair End雙端）測序原理：

• 必備名詞：

  flowcell： 測序反應(yīng)的載體/容器啃炸，1個flowcell有8個lane
  

  flowcell.png

  lane： 測序反應(yīng)的平行泳道，試劑添加卓舵、洗脫等過程的發(fā)生位置

  tile： 每次熒光掃描的位置南用，肉眼是看不到的

  雙端測序： 可能序列比較長有四五百bp，兩邊各測120-150bp

  junction： 雙端測序中間一些沒有測到的區(qū)域

  flowcell構(gòu)造：一個lane包含兩列（swath）掏湾，每一列有60個tile裹虫，每個tile會種下不同的cluster，每個tile在一次循環(huán)中會拍照4次（每個堿基一次）

• 邊合成變測序（sequence by synthesis, SBS）～合成

邊合成邊測序

<圖一>第一步： 構(gòu)建DNA文庫
超聲波將DNA分子打斷成300-800bp長序列片段（人類基因組打成300-500bp）融击，用酶補平為平末端筑公，然后3‘端加一個A堿基（因為接頭的3‘端有一個突出的T），再在兩端加上互補配對的adapter尊浪，再通過PCR擴增達(dá)到一定濃度匣屡，構(gòu)成單鏈DNA文庫。
【接頭設(shè)計很巧妙拇涤，兩個作用：
1. 實現(xiàn)橋式擴增捣作，高效；2. 可以實現(xiàn)雙端測序】

好奇鹅士，接頭怎么加上去的呢券躁？來源https://www.fimm.fi/en/services/technology-centre/sequencing/next-generation-sequencing/dna-library-preparation

Adding Adapter

另外，利用低循環(huán)擴增技術(shù)在接頭處進行修飾加上一些周邊掉盅。

Adapter修飾

<圖二>第二步： 上樣 【重點搞清楚lane上有哪兩種接頭也拜，待測序列含有哪兩種接頭，這很重要趾痘！】

• flowcell是用于吸附流動DNA片段的槽道慢哈，測序就在此進行浑吟。上面構(gòu)建好的文庫中的待測序列事先配置好一定的濃度贺纲，經(jīng)過這里的時候，會在特異的化學(xué)試劑作用下备典，強力隨機地附著在lane上瓦侮，與上面的短序列配對艰赞。上樣的結(jié)果就是lane吸附住了沖過來的DNA，并且可以在表面進行橋式PCR擴增肚吏。

• 【??這里有三點要注意：
  1. lane與lane之間一般不會相互影響方妖，也就是說一般不會出現(xiàn)lane1固定的DNA又與lane2結(jié)合。
  2. lane上隨機分布兩種接頭罚攀，p5‘（與P5互補）党觅，P7（與P7'互補）雌澄。
  待測序列自帶了p5接頭和p7接頭；
  3. 序列只能一開始是利用p5接頭互補杯瞻，因為p7接頭和lane是一樣的嘛】

<圖三>第三步：橋式PCR

• 第一輪擴增模版：flowcell表面固定的序列 --> 模版鏈

  第一輪結(jié)果：序列補成雙鏈镐牺。

  一個很重要的概念：我們要測序的是模版鏈p5 - p7，開始它與lane接頭配對產(chǎn)生了互補鏈魁莉，后來強堿試劑作用下兩條互補鏈被分開睬涧，由于模版鏈沒有附著在lane上，模版鏈被沖走旗唁，但是互補鏈依然穩(wěn)穩(wěn)固定在lane上畦浓。接下來就要對互補鏈p5‘- p7‘ 進行操作～

  擴增模版

• 去雜：加入NaOH強堿性溶液使雙鏈DNA變性，互補鏈由于和lane上短序列強力連接固定住了检疫；模板鏈?zhǔn)チ穗p鏈氫鍵連接讶请，好似懸空，它會被洗脫

• 橋式形成： 加入緩沖溶液屎媳，互補鏈的p7‘和lane上的p7互補（但還是一個lane中的）就像下圖這樣（摘自illumina官網(wǎng)）目的是快速擴增lane p7接頭連接的鏈夺溢，也就是下圖中的Forward Strand，它和我們的模版鏈?zhǔn)且恢碌慕宋Ｎ覀兒髞頊y序只用這一半

橋式形成

• 橋式PCR： PCR彎成橋狀企垦，一輪橋式擴增一倍

• 循環(huán)： 大約35個循環(huán)后，最終每個DNA片段都將在各自的位置上集中成束晒来，稱為cluster钞诡，這是一群完全相同的序列。目的在于實現(xiàn)放大單一堿基的信號強度湃崩，滿足后期測序需求

• 解鏈： 橋式PCR完成后荧降，形成了很多的橋形的互補雙鏈，再次強堿解鏈攒读。
  這一次不再進行復(fù)制朵诫，而是利用一種酶--甲酰胺基嘧啶糖苷酶（Fpg）選擇性的切掉lane 上p5‘ 連接的鏈，只留下了與lane p7連接的鏈即Forward Strand

<圖四>第四步：測序

測序

邊合成變測序～測序

一輪測序是這樣完成的：

雙端測序之Forward Strand

先是primer結(jié)合到靠近p5的sequencing primer binding site1上，再加入特殊的dNTP【它的3‘ 羥基被疊氮基團替代日缨，因此每次只能添加一個dNTP钱反；還含有熒光基團，能激發(fā)不同顏色】；

在dNTP被添加到合成鏈上后面哥，所有未使用的游離dNTP和DNA聚合酶會被洗脫掉哎壳；

再加入激發(fā)熒光緩沖液，用激光激發(fā)熒光信號尚卫，光學(xué)設(shè)備記錄熒光信號的記錄归榕，計算機將光學(xué)信號轉(zhuǎn)化為測序堿基

這一個循環(huán)就能測定flowcell上成千上萬的cluster，這就實現(xiàn)了高通量

再向下一輪：

再加入化學(xué)試劑淬滅熒光信號并使dNTP 3’ 疊氮基團變成羥基焕毫，這樣能繼續(xù)向下進行再加一個蹲坷，并且保證這個不再發(fā)出熒光驶乾。如此重復(fù)直至所有鏈的堿基序列被檢測出邑飒。得到了Forward Strand序列

因為一個cluster的序列是一樣的，所以理論上cluster的熒光顏色應(yīng)該一致级乐，下面是來自網(wǎng)站http://tucf-genomics.tufts.edu/home/faq的掃描圖片疙咸。本來儀器得出的是黑白的，顏色是后來計算機計算分析后加上去的

掃描圖片

Index測序： 上面的循環(huán)結(jié)束后风科，read product被沖掉撒轮，index1 primer和鏈上的index1 互補配對，進行index1的檢測贼穆。測完后题山，洗脫產(chǎn)物，得到index1 的序列故痊。接下來p5與lane上的p5‘配對顶瞳，測得了index2，并洗脫

Index測序

雙端測序之Reverse Strand： 洗脫掉index2 產(chǎn)物后愕秫，還是一個橋式擴增慨菱，得到雙鏈，再變性得到原始Forward strand 和 新的Reverse Strand戴甩， 除去測完的Forward strand符喝。然后和測Forward一樣，也是先連接primer甜孤，只是連接的位點是Primer Binding Site2协饲，測完后得到reverse strand序列

一個小補充：知道了PE Seq,那么單端SE 測序怎么測的呢？single-end只將index缴川，Primer binding site以及P7/P5添加到 fragamented DNA片段的一端茉稠，另一端直接連上P5/P7，將片段固定在Flowcell上橋式PCR生成DNA簇二跋，然后單端測序讀取序列

一點小問題:illunima一次只添加一個dNTP战惊，確實能夠保證準(zhǔn)確測量同聚物的長度問題（同聚物就是3‘端一個一個加dNTP聚合而成的聚合物，因為還帶接頭，所以不能直接說成DNA分子）吞获。
那么為何illumina會限制合成的鏈的長度呢况凉，不能像Sanger法一樣，最長測1k各拷？
原因就出在二代測序多出來的PCR過程：每一個位點都要測許多次刁绒，比如一段時間后的PCR得到的每個cluster都各包含200條完全相同的序列，那就需要對這200條序列的同一個位點進行測序烤黍。
第一輪我們來測第一個位點（假設(shè)位點1是A）正常來講知市，200條序列都應(yīng)該加A堿基，但是不巧只有199個在位點1都加了堿基A速蕊，有一條序列沒有加上嫂丙，所以就出現(xiàn)了199個紅色1個灰色【當(dāng)然目前還構(gòu)不成影響】；
第二輪（假設(shè)位點2是G）大家應(yīng)該都加G測得綠色规哲，但是之前的那個沒有加上A的跟啤，他要對之前的失誤進行補償，因此別的序列加G的時候唉锌，它加上了本該上次就加的A隅肥，它得到了紅色，這個紅色在一大群的綠色中就是作為雜信號存在的袄简。依次向下腥放，測序長度越長，雜信號越多绿语，最后可能標(biāo)準(zhǔn)信號和雜信號各一半秃症，這樣系統(tǒng)無法判斷，只能給N汞舱，而N多了對于后續(xù)的分析處理很麻煩伍纫，去了吧丟失數(shù)據(jù)，不去吧又是冗余昂芜。

Q:有同學(xué)問了: 上面說的莹规，依次向下，測序長度越長泌神，雜信號越多良漱。這個說不通啊，依次向下以后欢际，永遠(yuǎn)都只有1/200個雜信號啊母市，怎么會越來越多呢？
A:上面提到的只是假設(shè)只是一條序列出現(xiàn)失誤损趋，其實測序過程中隨著鏈的延長大家的錯誤都是不斷增加的患久，因為酶活性在降低。對于一個測序位點，第一輪可能有一條序列錯的蒋失，那么肯定他以后都會錯吧返帕；到了第二輪可能有兩條或者更多的序列出現(xiàn)沒加的現(xiàn)象，而他們也會一步錯篙挽，步步錯荆萤。

因此，錯誤來源就是某個位點有序列沒跟上隊伍铣卡，這次不加链韭，以后都是錯的，這樣的序列多了煮落，也就導(dǎo)致整體錯誤率升高

數(shù)據(jù)產(chǎn)生：

大體上我們平常使用的測序儀就是這樣（以Hiseq 2000為例）敞峭，后續(xù)升級版主要提高通量

Hiseq 2000測序儀

• 簡單了解下數(shù)據(jù)產(chǎn)生大體流程：

  從熒光信號的產(chǎn)生到堿基序列的識別這一過程，主要包括圖象校正（即空間校正）州邢、cluster識別儡陨、熒光校正（即光學(xué)校正）、phasing/prephasing（即化學(xué)校正）量淌、堿基識別、PF（Illumina默認(rèn)的數(shù)據(jù)過濾算法Pass Filtering）嫌褪、質(zhì)量評估等7個步驟

• 照相機是如何識別的呀枢？

  利用了CCD相機(1)對每一個簇(cluster)進行識別，確定其坐標(biāo)笼痛；(2)提取每個簇分別在A裙秋、G、C缨伊、T四個波長的信號強度值摘刑。另外拍照過程相當(dāng)耗時，一次循環(huán)所產(chǎn)生的信號需要40分鐘左右才能拍照收集完畢刻坊。使用相機的掃描功能會更快一些

相機識別

• 數(shù)據(jù)量產(chǎn)出：測序儀搭配了兩個flowcell枷恕，簡稱雙流動槽。比較經(jīng)典的Hiseq2500一次能產(chǎn)出700-800Gb數(shù)據(jù)（此處Gb為測序堿基數(shù)谭胚，不同于字節(jié)數(shù)的Gb）

  關(guān)于數(shù)據(jù)轉(zhuǎn)換徐块，舉個例子比較好理解：Gb是測序中的數(shù)據(jù)量，1 Gigabase= 十億堿基灾而。人類全基因組測序得到了90G的原始數(shù)據(jù)胡控，也就是900億堿基。得到的900億堿基旁趟，也對應(yīng)900億個質(zhì)量值昼激，加起來就是1800億個字符。想一想fastq的格式：第一行是測序說明，一般是45個字符橙困，也就是說敛劝，每一條測序reads中第一行就有大概45個字符。那么多少條reads呢纷宇？根據(jù)PE150計算：測序策略是一條reads包括150bp夸盟，現(xiàn)在900億堿基，就對應(yīng)900億/150=60億條reads 像捶。因此第一行總字符是：60億*45=270億個字符上陕。注意到fastq文件共四行，其中1拓春、2释簿、4行的總數(shù)量分別為270億、900億硼莽、900億庶溶，第三行就是一個+，基本可以忽略不計懂鸵。加起來總共2070億字符偏螺。計算機中，根據(jù)編碼規(guī)則不同匆光，字符與字節(jié)對換關(guān)系不同套像。

  
  

  
  

  Fastq文件是ASCII編碼文件，其中每一個字符就對應(yīng)一個ASCII碼终息，也就等于一個字節(jié)夺巩。計算機的1 GB（Gigabytes） 是1024^3 個字節(jié)
  因此，二者對換關(guān)系就是：全基因組測序的90Gb對應(yīng)（2070x10^8 /1024^3 ）=193GB計算機存儲空間周崭。
  或者更快的計算： 測序報告會給出reads數(shù)柳譬，如果測序策略是PE150，那么占用硬盤空間大小就是n(reads)(150+150+45)/1024^3
  另外续镇，測序儀下機后的數(shù)據(jù)都是用gz壓縮后的文件.fastq.gz美澳，能壓縮2.7倍，大概71G左右磨取。

看一下來自illumina官網(wǎng)的統(tǒng)計數(shù)據(jù)人柿，換算一下大概能知道，高配的HiSeq X10每臺機器每輪測序所產(chǎn)生的數(shù)據(jù)量大約是1.6Tb忙厌，即1600 Gb凫岖，2017年的JP摩根大會上發(fā)布的NovaSeq 6000一次測序保守估計能得到3000G的測序數(shù)據(jù)。如果要做人類基因組重測序逢净，Hiseq系列一般需要3-5天哥放，而Nova只需要40小時歼指！

illumina測序儀產(chǎn)出統(tǒng)計

數(shù)據(jù)初步分析：

使用fastqc進行質(zhì)量分析，這是一款Java軟件甥雕，支持多線程踩身。寫這篇文章的時候版本是v0.11.7。

軟件前期準(zhǔn)備：

• 下載方式有兩種：

1. 官網(wǎng)下載好用filezilla導(dǎo)入linux服務(wù)器

2. 直接在服務(wù)器中wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.7.zip

• 接著安裝unzip fastqc_v0.11.7.zip --> cd FastQC -- > chmod755 fastqc

最后這個chmod有必要說一下社露，這個權(quán)限管理命令

chmod 用3個數(shù)字來表達(dá)對 用戶（文件或目錄的所有者）挟阻，用戶組（同組用戶），其他用戶 的權(quán)限：
　　如：chmod 755 fastqc
　　數(shù)字7是表達(dá)同時具有讀峭弟，寫附鸽，執(zhí)行權(quán)限：(7 = 4 + 2+ 1)
　　讀取--用數(shù)字4表示；
　　寫入--用數(shù)字2表示瞒瘸；
　　執(zhí)行--用數(shù)字1表示; 三者皆否：0

• 設(shè)置完權(quán)限后坷备，還需要將FastQC文件夾（這里請注意是文件夾，而非fastqc這個可執(zhí)行程序）導(dǎo)入環(huán)境變量

  echo 'export PATH=/YOUR/FASTQC PATH/:$PATH' >> ~/.bashrc

  再source ~/.bashrc

• 檢查軟件是否安裝成功 fastqc --help 出現(xiàn)幫助信息就可以使用啦情臭！

后期軟件使用：

   基本格式 （各種參數(shù)+多個文件～支持多線程）
   
   fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] seqfile1 .. seqfileN
   
   -o --outdir FastQC生成的報告文件的儲存路徑
   --extract 使用這個參數(shù)是讓程序不打包【默認(rèn)會打包成一個壓縮文件】
   -t --threads 選擇程序運行的線程數(shù)省撑，每個線程會占用250MB內(nèi)存（一般與文件數(shù)量一致就好）
   -q --quiet 安靜運行模式【不選這個選項，程序會實時報告運行的狀況】

結(jié)果分析：

• 如果你有自己的轉(zhuǎn)錄組或者其他數(shù)據(jù)俯在，可以現(xiàn)在測試了
• 如果沒有竟秫，想學(xué)一下這個軟件流程以及如何解讀結(jié)果，可以下載公共數(shù)據(jù)（下載兩個雙端測序文件共4個）

# 順便說一下這里用到了curl -O（保留遠(yuǎn)程文件的文件名） -L（對于自動跳轉(zhuǎn)的網(wǎng)頁朝巫，curl 就會跳轉(zhuǎn)到新的網(wǎng)址）
# 當(dāng)然用wget也可以鸿摇，至于二者區(qū)別，可以參考https://www.cnblogs.com/lsdb/p/7171779.html
    curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_1.fastq.gz
    curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_2.fastq.gz
    curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1977249_1.fastq.gz
    curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1977249_2.fastq.gz
    
# 下載完可以檢查一下數(shù)據(jù)完整性劈猿，當(dāng)然不是必須
    md5sum *.gz
    
# 質(zhì)控四個文件，我們可以采用四線程
# 大概用時   real   0m23.344s （如果你也想統(tǒng)計時間潮孽，就在命令前加time）
    fastqc *.gz -t 4
#將結(jié)果.html用filezilla導(dǎo)出揪荣，瀏覽器查看

• 首先看到的是一個總結(jié)報告，
  左邊這一欄會告訴你往史，哪些是提出警告的（??表示）仗颈，那些是fail的（?）
  
  QC Summary

• 接下來一部分一部分解析，共10部分

1. Basic Statistics 基本信息

   Basic Statistics

   • Encoding: 測序平臺編號椎例，現(xiàn)在Sanger/ Illumina 1.8以上都是Phred 33編碼

   • Total sequences: reads數(shù)量（reads就是高通量測序平臺產(chǎn)生的序列標(biāo)簽挨决，翻譯為讀段！）

   • Sequence length: 測序長度

   • %GC: GC含量： 需要重點關(guān)注订歪，可以幫助區(qū)別物種脖祈，人類細(xì)胞42%左右

2. Per base sequence quility：每個測序read各堿基質(zhì)量【十分重要！】

Per base sequence quility

• 橫軸：測序序列的1-251個堿基刷晋；

• 縱軸：質(zhì)量得分盖高，score = -10 * log10（error）慎陵，例如錯誤率error為1%，那么算出的score就是20

• 箱線圖boxplot：對每一個堿基的質(zhì)量的統(tǒng)計喻奥。箱子上面的須（up bar）為90%分位數(shù)席纽，下面的須（down bar）為10%分位數(shù)，箱子中的紅線為中位數(shù)即50%分位數(shù)撞蚕，箱子頂（upside）為75%分位數(shù)润梯，箱子低（downside）為25%分位數(shù)。這個boxplot的意義：一是看數(shù)據(jù)是否具有對稱性甥厦；二是看數(shù)據(jù)分布差異纺铭，這里主要利用了第二點。bar的跨度越大矫渔，說明數(shù)據(jù)越不穩(wěn)定彤蔽。

• 藍(lán)色的線將各個堿基的質(zhì)量平均值連接起來

• 解釋一下：圖中藍(lán)線的走勢為何先高后低？因為目前采用的邊合成邊測序使用的是化學(xué)方法促使鏈由5'向3'延伸庙洼，也就是利用了DNA聚合酶顿痪。剛開始測序，合成反應(yīng)還不是很穩(wěn)定油够，但是酶的質(zhì)量還很好蚁袭，所以會在高質(zhì)量區(qū)域內(nèi)有一定的波動（這里的1-30bp），后來穩(wěn)定了石咬，但是隨著時間的推移揩悄，酶的活力逐漸下降，特異性也變差鬼悠，所以越往后出錯幾率越大删性。
  【就像一個司機開車，一開始小心謹(jǐn)慎焕窝，起步慢蹬挺，開的也慢，慢慢提速它掂。后來越開越帶勁巴帮，但是也越來越困，疲勞駕駛?cè)菀壮鍪?/strong>】

• 一般能用的數(shù)據(jù)都要求至少Q(mào)20虐秋，也就是下四分位（10%分位數(shù)）的質(zhì)量值要大于20榕茧。因此這里的189bp后面的需要切除，才能繼續(xù)分析

• 二代測序客给，最好是達(dá)到Q20的堿基要在95%以上（最差不低于90%）用押，Q30要求大于85%（最差也不要低于80%）

3. Per sequence quility scores：每條序列 質(zhì)量統(tǒng)計

Per sequence quility scores

• 橫軸：質(zhì)量值0-40，也即是Q值

• 縱軸：每個質(zhì)量值對應(yīng)的read數(shù)

• 我們的例子中一條read有251bp起愈， 那么其中任意一條的251bp的質(zhì)量平均值就是這條read的質(zhì)量值只恨。只要大部分都高于20說明比較正常

4. Per base sequence content：read各個位置堿基比例分布

Per base sequence content

• 橫軸：各堿基位置译仗；縱軸：堿基百分比

• 四條線四種顏色代表四種堿基在每個位置的平均含量（一個位置會測很多reads，然后求一個平均）

• 一般來講官觅，A=T纵菌， C=G， 但是剛開始測序儀不穩(wěn)定可能出現(xiàn)波動休涤，這是正常的咱圆。一般不是波動特別大的，像這里cut掉前5bp就夠了功氨。另外如果A序苏、T 或 C、G間出現(xiàn)偏差捷凄，只要在1%以內(nèi)都是可以接受的

5. Per sequence GC content： 序列平均GC分布

Per sequence GC content

• 橫軸為平均GC含量忱详； 縱軸為每個GC含量對應(yīng)的序列數(shù)量

• 藍(lán)線為系統(tǒng)計算得到的理論分布；紅線為測量值跺涤，二者越接近越好

• 這里不相符可能有兩個原因：

  1. 前面提到了匈睁，GC可以作為物種特異性根據(jù)，這里出現(xiàn)了其他的峰有可能混入了其他物種的DNA桶错；

  2. 目前二代測序基本都會有序列偏向性(所說的 bias)航唆，也就是某些特定區(qū)域會被反復(fù)測序，以至于高于正常水平院刁，變相說明測序過程不夠隨機糯钙。這種現(xiàn)象會對以后的變異檢測以及CNV分析造成影響

• 如果出現(xiàn)怎么辦？-- 把和我們使用物種GC-content有差異的reads拿出來做blast退腥，來確認(rèn)是否為某些雜菌

6. Per base N content： N含量分布

Per base N content

• N是指儀器不能識別ATCG時給出的結(jié)果任岸，一般不會出現(xiàn)。但是如果出現(xiàn)并且量還很大狡刘，應(yīng)該就是測序系統(tǒng)或者試劑的問題

• 任意位置的N的比例超過5%演闭，報"WARN"；任意位置的N的比例超過20%颓帝，報"FAIL"

7. Sequence length distribution： 序列長度統(tǒng)計

Sequence length distribution

• 理想情況下，測得的序列長度應(yīng)該是相等的窝革。實際上總有些偏差

• 當(dāng)reads長度不一致時報"WARN"购城；當(dāng)有長度為0的read時報“FAIL”

• 這里顯示大部分都落在251bp這個測序長度上，有少量為250或252bp虐译，但這不影響瘪板；如果偏差很大就不可信了

8. Sequence duplication level：統(tǒng)計序列完全一樣的reads的頻率

Sequence duplication level

• 橫坐標(biāo)是duplication的次數(shù)；縱坐標(biāo)是duplicated reads的數(shù)目（紅線）

• 解釋下橫坐標(biāo)為何會有>10, >50等出現(xiàn)：測序的原始數(shù)據(jù)很大漆诽，如果每一條reads都統(tǒng)計侮攀，將耗時很久锣枝。這里軟件只采用了數(shù)據(jù)的前200,000條reads統(tǒng)計其在全部數(shù)據(jù)中的重復(fù)數(shù)目,另外大于75bp的reads只取50bp進行比較。重復(fù)數(shù)大于10的reads被合并統(tǒng)計成了>10兰英，以此類推...

• unique reads總數(shù)（藍(lán)線）作為100%撇叁，上圖中可以看出，大概僅有2%的uniqe reads可以觀察到兩次重復(fù)畦贸。也就是說陨闹，我們這里的非unique reads占總數(shù)比例僅有2%左右。

• 正常情況下的確薄坏，測序深度越高趋厉，越容易產(chǎn)生一定程度的duplication。高程度的duplication level胶坠，提示我們可能有bias的存在（如建庫過程中的PCR duplication）君账。

  另外和做的項目也有關(guān)，一般轉(zhuǎn)錄組測序的結(jié)果中duplication level都比較高沈善，60-70%都正常乡数，這是因為轉(zhuǎn)錄組測的是基因的覆蓋深度，各個基因表達(dá)量不同矮瘟，如果某個基因覆蓋度較高【tip：覆蓋度是指基因/轉(zhuǎn)錄組測序測到的部分占整個組的比例】瞳脓，那么測的部分就越多，相對應(yīng)的duplication也會更高澈侠；對于外顯子組測序來說劫侧，一般覆蓋度比較一致，這里出現(xiàn)了duplication就不太正常哨啃。

• 當(dāng)非unique的reads占總數(shù)的比例大于20%時烧栋，報"WARN"；當(dāng)非unique的reads占總數(shù)的比例大于50%時拳球，報"FAIL“

9. Overrepresented sequences：大量重復(fù)序列

Overrepresented sequences

• 和第8個duplication計算一樣审姓，也是取前200,000進行統(tǒng)計，大于75bp只取50bp祝峻。

• 發(fā)現(xiàn)超過總reads數(shù)0.1%的reads時報”WARN“魔吐，當(dāng)發(fā)現(xiàn)超過總reads數(shù)1%的reads時報”FAIL“

10. Adapter content: 接頭含量

    Adapter content

• 表示序列中兩端adapter的情況

• 軟件內(nèi)置了四種常用的測序接頭序列， fastqc 有一個參數(shù)-a可以自定義接頭序列

• 此圖中使用的illumina universal adapter并未去除莱找，后期再使用trimmomatic/cutadapt去接頭

• 我們在學(xué)習(xí)測序原理的時候知道了酬姆，接頭的作用一個是能夠使得序列結(jié)合到flowcell上寞肖，另外一個多樣本測序時利用接頭旁邊的index加以區(qū)分

• 什么情況下能夠測到接頭呢锹淌？ 一般測序read的長度大于插入片段（待測序列）的長度時會發(fā)生。對于WGS建庫測序來講米罚，一般不會發(fā)生浮定，因為他的待測序列要幾百bp相满，而測序也就150bp算高了层亿。但是對于RNA-seq，一般測序序列比較短立美，尤其是miRNA匿又，只有幾十bp，這是就會測到read末尾的接頭

11. （還有一類這里沒體現(xiàn)）Kmer content： 重復(fù)短序列

Kmer content

• 表示：在序列中某些特征的短序列重復(fù)出現(xiàn)的次數(shù)

• 這個圖是轉(zhuǎn)錄組測序的一個文件悯辙，可以看到6-9bp幾種短序列都出現(xiàn)了好多次琳省。出現(xiàn)的原因可能是：

  1. 沒有去除軟件內(nèi)置的adapter或者沒有使用-a參數(shù)自定義adapter

  2. 序列本身重復(fù)度較高，例如在建庫PCR過程出現(xiàn)序列偏向性bias--> 這在轉(zhuǎn)錄組測序中確實存在

---

數(shù)據(jù)過濾：

主要針對接頭序列和低質(zhì)量序列

• 工具： 有許多工具能干這事躲撰，例如SOAPnuke针贬、cutadapt、untrimmed拢蛋、sickle和seqtk等桦他，其中經(jīng)常用到的是Trimmomatic（也是一個java程序）

• Trimmomatic:

  • 安裝：官網(wǎng)下載我使用的是0.36版本
    Trimmomatic官網(wǎng)
    
    下載后直接解壓安裝，其中有一個trimmomatic-0.36.jar是執(zhí)行文件谆棱，使用時輸入java -jar trimmomatic-0.36.jar就可以快压。另外還有一個adapter文件夾，里面存放了常用的illumina測序儀接頭序列fasta格式,后續(xù)處理接頭需要制定具體文件垃瞧。
    自帶接頭序列

• 關(guān)于adapter： 目前絕大部分的illumina的Hiseq和Miseq系列使用的都是Truseq3蔫劣，過去的GA2測序儀使用的是Truseq2，PE/SE對應(yīng)單端還是雙端測序 个从， 如果使用的不是illumina測的脉幢，可以按照里面的格式自己新建一個接頭文件，但其中的命名要注意嗦锐。詳情見官網(wǎng)主頁
  自定義Adapter說明

• 以PE測序為例說一下命令參數(shù)設(shè)置：

java -jar <path to trimmomatic.jar> PE [-threads <threads 線程數(shù)] 
        [-phred33 | -phred64] 質(zhì)量體系嫌松，默認(rèn)64，但我們現(xiàn)在一般都是33
        [-trimlog <logFile>] log文件
        <input 1> <input 2> 雙端測序原始fq文件
        <paired output 1> <unpaired output 1>  雙端1輸出文件 奕污、 過濾掉的文件
        <paired output 2> <unpaired output 2>  雙端2同理
        <step > 
        <詳細(xì)講一下step：>
        1. ILLUMINACLIP： 根據(jù)上面qc部分的測試數(shù)據(jù)萎羔，我們設(shè)置TruSeq2-PE.fa:2:40:15
         TruSeq2-PE.fa是接頭序列；
         2是比對時接頭序列時所允許的最大錯誤數(shù)碳默；
         40是要求PE的兩條read同時和adapter序列比對贾陷，匹配度加起來超40%，那么就認(rèn)為這對PE reads含  有adapter嘱根，并在對應(yīng)的位置需要進行切除昵宇；
        【那么為何不是匹配100%？因為測序時并不是把 adapter全測了儿子，只測了一部分】
         15 指的是只要某條read的某部分與adapter超過了15%的相似度就認(rèn)為包含，就要去除
        2. SLIDINGWINDOW： 滑動窗口長度 我們設(shè)置為4:20砸喻，代表窗口長度為4柔逼，窗口中的平均質(zhì)量值至少為20蒋譬，否  則會開始切除
        3. LEADING，指的是read開頭的堿基是否要被切除的質(zhì)量閾值愉适，這里設(shè)為2
        4. TRAILING犯助，指的是read末尾的堿基是否要被切除的質(zhì)量閾值，這里設(shè)為2
        5. MINLEN维咸，指的是read被切除后至少需要保留的長度剂买，如果低于該長度，會被丟掉癌蓖，這里設(shè)置25
        ?

#一個??范例(這個測試數(shù)據(jù)是phred 64的瞬哼，所以使用默認(rèn)值就好)：
        java -jar trimmomatic-0.36.jar PE -threads 8 \
        -trimlog logfile \
        reads_1.fq.gz reads_2.fq.gz \
        out.read_1.fq.gz out.trim.read_1.fq.gz \
        out.read_2.fq.gz out.trim.read_2.fq.gz \  ILLUMINACLIP:/path/Trimmomatic/adapters/TruSeq2-PE.fa:2:40:15 \
        SLIDINGWINDOW:4:20 \
        LEADING:2 TRAILING:2 \
        MINLEN:25</pre>

# 當(dāng)然，這只是對一個樣本的雙端測序文件進行操作租副，那么問題來了：
        # 如果你的樣本比較多呢坐慰？還要手動一個個輸入嗎？
        # 雖然說vim中使用快速替換 :%s/AA/BB/g 可以全局替換AA成BB用僧，但還是有點麻煩
        # 能不能讓程序來一個自動化呢结胀？是可以的！可以看作是上面??腳本的改進版～
        # 在腳本中構(gòu)建for循環(huán)
        # vi trim.sh
        ?
         #開頭這樣寫而不寫*.gz的目的是避免了樣本名稱的重復(fù)
        for filename in *_1.fastq.gz  
        do
         #對于filename（%）向右匹配_的全部內(nèi)容(*),然后這部分去掉责循，留下這之前的
        base=${filename%_*}  
        java -jar trimmomatic-0.36.jar PE \       # 加上反斜線能讓程序整體更清晰
         -threads 8 \
         -trimlog logfile \
         ${base}_1.fastq.gz ${base}_2.fastq.gz \
         out.${base}_1.fastq.gz out.trim.${base}_1.fastq.gz \
         out.${base}_2.fastq.gz out.trim.${base}_2.fastq.gz \
         ILLUMINACLIP:/home/u1239/biosoft/Trimmomatic-0.36/adapters/TruSeq2-PE.fa:2:30:10 \
         SLIDINGWINDOW:5:20 \
         LEADING:5 TRAILING:5 \
         MINLEN:50
        done

多樣本trim

• 如果報錯糟港，就說明參數(shù)錯誤例如我在運行的時候就有的問題百思不得解，敲了好幾遍代碼應(yīng)該沒錯院仿，但就是報錯秸抚，換個服務(wù)器就好了，難道是平臺問題意蛀？其實并不是,就是因為一點點小地方,有時只是一個代表換行的反斜線\
  報錯源于參數(shù)

---

多樣本質(zhì)控：

fastqc對于一兩個樣品還能吃得消耸别，但樣本多了，我們?nèi)绾瓮瑫r查看對比他們的信息呢县钥？

這里就可以使用界面更加優(yōu)美的multiqc軟件了, 他就是fastqc結(jié)果的整合我會從一個初學(xué)者角度來一步步進行秀姐，其中包括了中間犯的錯誤及改正～，如果你也在運行的過程發(fā)現(xiàn)了類似的錯誤若贮，可以參考一下省有。

---

分界線1: 嘗試自己安裝

---

下載地址：

https://files.pythonhosted.org/packages/fa/3e/fbdcbaebadad2110eed3b4212ced58617cbd9badbfc728c5eabc53916b84/multiqc-1.5.tar.gz

解壓縮后，會發(fā)現(xiàn)和以前安裝的源代碼文件不同谴麦，他沒有直接顯示可執(zhí)行文件蠢沿，也沒有配置文件。這是因為multiqc是一個python軟件匾效，這里先看一下setup.py：

MultiQC setup.py

直接使用python setup.py會報錯舷蟀，因為可能你的服務(wù)器并沒有setuptools，這是python依賴的第三方庫。

先安裝setuptools：

下載地址：

https://files.pythonhosted.org/packages/1a/04/d6f1159feaccdfc508517dba1929eb93a2854de729fa68da9d5c6b48fa00/setuptools-39.2.0.zip

解壓縮完setuptools會看到有這些文件野宜，這也是一般的python軟件常見的：包括了setup.py扫步，easy_install.py等

一般的Python程序

# 主要使用setup.py
python setup.py build # 編譯
python setup.py install #安裝

在第二步安裝的時候報錯,原因是不能對/usr/local/lib下的python進行操作，因為不管/usr/bin還是/usr/local/bin匈子，都是可讀不可改河胎，如果自己家目錄環(huán)境變量中查不到python的路徑，那么安裝過程中會自動調(diào)取更上層目錄的python使用虎敦。這往往會引發(fā)一系列問題游岳。如果要自己更改的話，需要自己home目錄下有python

報錯信息

好吧其徙，那么接下來我們再安裝自己的python胚迫，畢竟自己的軟件用起來方便：

# 下載地址 
wget https://www.python.org/ftp/python/3.6.5/Python-3.6.5.tgz
# .tgz 就等同tar.gz
# 解壓后會看到配置文件configure，設(shè)置成自己的軟件環(huán)境變量即可
./configure --prefix=/YOUR PATH/
make
make install
# 可能結(jié)尾會有報錯擂橘，但是不影響使用晌区，出現(xiàn)了可執(zhí)行文件python就成功了
# 將python復(fù)制到你的環(huán)境變量中就能直接調(diào)取
# 最后使用which python檢測是否安在了自己的目錄下

自己的python安裝完成了，進入到解壓好的multiqc文件夾下, python setup.py build 編譯會運行成功通贞，然后python setup.py install

安裝成功

接著朗若，multiqc就會出現(xiàn)在了自己的軟件環(huán)境變量中了，輸入multiqc就會看到

multiqc

---

分界線2: 同一件事昌罩，換個玩法

---

目的：安裝multiqc

途徑：conda自動安裝

緣由：好久不用conda安裝軟件哭懈，一直堅持源代碼，因為好掌控茎用。但是有的軟件依賴的包很多遣总，自己又很難發(fā)現(xiàn)這些潛在的關(guān)系，所以想重新試試conda轨功，但并不是傻瓜式的使用旭斥，而是讓conda聽我的話～授之以魚，不如授之以漁

Here we go古涧！

安裝： 安裝過程需要注意

1. 可以自定義新文件夾垂券，默認(rèn)是在home/miniconda3/

2. ??不要將conda添加到PATH中，我們只需要把它視作一個保姆babysitter羡滑，而不要當(dāng)一個管家housekeeper

3. 安裝完conda菇爪，把miniconda/bin下的conda復(fù)制到你的/biotools/PATH（這個環(huán)境變量因人而異）

添加源： 清華源和中科大源都不錯，別忘了再添加一個bioconda源

新建conda專屬下載目錄： 你可以在你的biotools目錄下新建一個conda軟件存放目錄柒昏，例如conda_install凳宙。然后把這個文件夾添加到環(huán)境變量。以后你用conda安裝的所有軟件都存放在這里职祷，

conda專屬目錄

??和你自己安裝的軟件要區(qū)分開氏涩， 然后利用conda install -p /PATH/conda_install multiqc 就可以實現(xiàn)multiqc的安裝了

---

寫在后面

---

1. multiqc的使用很簡單届囚，然后結(jié)果還是交互式的，能導(dǎo)出很多格式

2. 這次實現(xiàn)了人工和自動二者的有機結(jié)合：復(fù)雜的削葱，多依賴性的軟件還是要靠conda解決奖亚。但是如果一個軟件用兩種方法都安裝了，怎么選擇析砸？這個利用bashrc中設(shè)置的順序：哪個目錄優(yōu)先，就先搜索哪個目錄

結(jié)果一目了然：trim掉低質(zhì)量序列和接頭后爆袍，新的數(shù)據(jù)質(zhì)量值都在q30以上

MultiQC 1

MultiQC 2

歡迎關(guān)注我們的公眾號～_～　　
我們是兩個農(nóng)轉(zhuǎn)生信的小碩首繁，打造生信星球，想讓它成為一個不拽術(shù)語陨囊、通俗易懂的生信知識平臺弦疮。需要幫助或提出意見請后臺留言或發(fā)送郵件到Bioplanet520@outlook.com

bioinfoplanet