- 更新于2020.10.27
-
二代測(cè)序方式分為三種：
- single read 單端測(cè)序
- paired-end read 雙端測(cè)序
- mate-pair read 配對(duì)測(cè)序
雖然有三種方式，但是大多數(shù)數(shù)據(jù)（包括本課題）為雙端測(cè)序欣鳖，至于三者之間的差別察皇、優(yōu)缺點(diǎn)以及適用場(chǎng)合可以自行搜索。接下來講的內(nèi)容都是基于雙端測(cè)序的泽台。
-
測(cè)序儀原始下機(jī)的數(shù)據(jù)我們稱為raw data什荣，二代測(cè)序是將DNA片段打斷了再測(cè)的，每個(gè)測(cè)序片段我們稱為read怀酷，質(zhì)控完的數(shù)據(jù)稱為clean data稻爬。既然是雙端測(cè)序，那么文件就是成對(duì)出現(xiàn)的蜕依，分別記錄reads兩端的信息：一般的命名是*.fq1.gz桅锄、*.fq2.gz（' * ' 表示通配符），這是一個(gè)fastq文件样眠，通常以fq或fastq作為后綴友瘤，具體內(nèi)容下方會(huì)介紹。這里我給大家展示我們課題的一個(gè)真實(shí)數(shù)據(jù)檐束，給大家看看大小以及命名方式辫秧。
原始數(shù)據(jù)展示 我們可以看到，首先我把我的用戶名打碼了被丧，呵呵傅是，開個(gè)玩笑吵取。
首先這個(gè)樣本總共有五個(gè)文件，這就是公司給的原始數(shù)據(jù)，我原封不動(dòng)的上傳到了我的服務(wù)器上嫌术，也沒改名也沒怎樣渗饮，我在圖片上畫了五個(gè)藍(lán)色的框这橙，下面是分別代表的意思：
W2018001：是樣本的名字革为，其實(shí)我一開始提交給公司的名字是2018001，為什么多了個(gè)W我也不知道啊糕簿。
NZTD...：這一串是公司自動(dòng)生成的編號(hào)探入，他們內(nèi)部使用的，不知道也不需要知道啥意思懂诗，不同的公司蜂嗽，可能沒有這個(gè)
HCV...：flowcell_ID的信息，一般情況下一個(gè)樣本是一樣的
L1：flowcell_lane的信息殃恒，圖中有L1和L2植旧，這也是為什么一個(gè)數(shù)據(jù)他給我四個(gè)文件的原因辱揭，他在不同的lane上測(cè)的，這里還要注意的一點(diǎn)就是病附，這個(gè)lane的值可能是同一個(gè)问窃，就是即使碰到都是L1也是可能的
1，2：分別代表reads的兩端
- 我們可以看到完沪，這個(gè)30X的WGS的數(shù)據(jù)量還是很大的域庇，所以為什么他會(huì)是壓縮文件，不同公司給的命名可能不一樣覆积，不過大體不會(huì)相差太多听皿，具體命名的含義也可以問問公司的。（注：flowcell_ID宽档、flowcell_lane等信息一般在原始數(shù)據(jù)的內(nèi)也可以看到）
- 這里給到的是一個(gè)樣本2對(duì)文件尉姨，實(shí)際情況是，可能有多對(duì)或者只有一對(duì)文件吗冤，對(duì)于這樣的情況又厉，我們的處理方式一般是看作一個(gè)文件處理，就是merge一下椎瘟，有見過有些課題組是在這一步直接merge馋没，我的處理方式是后期bam文件再merge。
- 好了降传，拿到這么一個(gè)數(shù)據(jù)，我們做的第一件事情應(yīng)該是校驗(yàn)數(shù)據(jù)的完整性９磁Ｆ排拧！這點(diǎn)很重要笔链，即使一般情況下段只，傳輸都不會(huì)出錯(cuò)，但這一步還是要做的鉴扫，怎么去做呢赞枕，就看這個(gè)截圖里的第一個(gè)文件MD5.txt，其實(shí)這就是一個(gè)文本文件坪创，打開了看看是這樣的
- 就是這樣的四行信息炕婶，也就是一個(gè)文件名對(duì)應(yīng)一個(gè)校驗(yàn)碼，最簡(jiǎn)單的校驗(yàn)方式莱预，linux系統(tǒng)自帶的MD5校驗(yàn)命令：
$ md5sum -c MD5.txt
- 這個(gè)還是需要一點(diǎn)點(diǎn)時(shí)間的柠掂，大家要耐心等待。展示一下我這里的結(jié)果吧（其實(shí)在上傳完數(shù)據(jù)的時(shí)候我就校驗(yàn)過）
好了依沮，校驗(yàn)完了涯贞，沒問題枪狂，那么就給大家介紹一下文件的格式。
-
fastq文件是什么呢宋渔，其實(shí)就是文本文件州疾，可以直接查看，以下是示例：fastq格式
一條read一條記錄皇拣，一條記錄占四行严蓖，第一行是注釋，第二行是序列审磁，就是我們所說的ATCG堿基序列谈飒，第三行是‘+’，第四行是對(duì)應(yīng)的每個(gè)堿基的測(cè)序質(zhì)量态蒂，也就是fastq中的“q”杭措。每條記錄之間是沒有空格的。
貼一下關(guān)于第一行的解釋
| EAS139 | The unique instrument name |
|---|---|
| 136 | Run ID |
| FC706VJ | Flowcell ID |
| 2 | Flowcell lane |
| 2104 | Tile number within the flowcell lane |
| 15343 | 'x'-coordinate of the cluster within the tile |
| 197393 | 'y'-coordinate of the cluster within the tile |
| 1 | Member of a pair, 1 or 2 (paired-end or mate-pair reads only) |
| Y | Y if the read fails filter (read is bad), N otherwise |
| 18 | 0 when none of the control bits are on, otherwise it is an even number |
| ATCACG | Index sequence |
第一個(gè)是測(cè)序機(jī)器ID钾恢，獨(dú)一無二的手素，你可以去illumina官網(wǎng)查到的（我沒查過）
第二個(gè)是這個(gè)機(jī)器跑的次數(shù)，一般機(jī)器都是有使用壽命的瘩蚪，所以次數(shù)越多肯定結(jié)果越差泉懦，那么一般在什么區(qū)間合適呢，我的老師給的建議是200-9999疹瘦，為啥還有下限呢崩哩，那是因?yàn)椋瑱C(jī)器剛開始跑的那幾次言沐，需要磨合嘛邓嘹，不是很準(zhǔn)∠找龋看來這個(gè)展示數(shù)據(jù)不是很好啊汹押，嚇得我趕緊去看了下我的數(shù)據(jù)，還好起便，在212次（懶得貼圖了棚贾，簡(jiǎn)書弄個(gè)圖真麻煩），在范圍內(nèi)榆综。
倒數(shù)第三個(gè)妙痹，表示數(shù)據(jù)有沒有被過濾過，一般我們的原始數(shù)據(jù)肯定是N的鼻疮，要是給的原始數(shù)據(jù)是Y细诸，你就得好好問問公司原因了。
其余具體的我就不解釋了陋守，有問題可以評(píng)論交流震贵。
第二行要注意的是利赋，可能有N的出現(xiàn)，那是因?yàn)橛行A基沒被測(cè)出來猩系。
第四行是ASCII碼表示的堿基質(zhì)量媚送，這樣就能保證質(zhì)量是用一個(gè)字符表示的，和堿基一一對(duì)應(yīng)寇甸。公式如下：
????????????P=10^-[(n-33)/10]舉個(gè)栗子：比如“?”在ASCII碼表上對(duì)應(yīng)的編號(hào)是63塘偎，那么n就是63，減去33以后就是30拿霉，也就是我們說的Q30了吟秩，所以
Q = n - 33，P算出來就是0.001绽淘，這個(gè)就是錯(cuò)誤率涵防，反過來，準(zhǔn)確率就是99.9%沪铭，Q30就是準(zhǔn)確率99.9%壮池，同理Q20就是準(zhǔn)確率99%。
- ok杀怠，介紹完格式椰憋，介紹兩款質(zhì)控的軟件，fastqc和fastp
fastqc：
$ wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.7.zip
$ unzip fastqc_v0.11.7.zip
$ cd FastQC/
$ chmod 755 fastqc
fastp：
$ wget http://opengene.org/fastp/fastp
$ chmod 755 fastp
-
以上是我的安裝方式赔退，fastqc是一個(gè)很常用的軟件橙依，我想只要不是小白都用過。對(duì)于ubuntu用戶請(qǐng)用我提供的方式進(jìn)行安裝硕旗，用
apt-get install fastqc可能在使用的時(shí)候出現(xiàn)如下報(bào)錯(cuò)
fastqc錯(cuò)誤 至于fastp窗骑，是用來清洗質(zhì)量不好的reads的，其實(shí)類似的軟件很多卵渴，包括trimmomatic等，之所以選擇這個(gè)是因?yàn)橛眠^那么多軟件后鲤竹，發(fā)現(xiàn)fastp無論在安裝還是使用上都很順手浪读，僅此而已。（注：fastp的UMI操作可能導(dǎo)致結(jié)果不能復(fù)現(xiàn)辛藻，詳見issue）
好了碘橘，現(xiàn)在讓我們來使用fastqc對(duì)raw data的質(zhì)量進(jìn)行分析并查看
$ fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN
#這是fastqc的使用方法，其實(shí)很簡(jiǎn)單吱肌，對(duì)于不懂的軟件痘拆，我一般的推薦是，先看官方說明氮墨，再逛逛論壇纺蛆，這里不展開說這個(gè)軟件怎么使用了吐葵，下面是我的使用代碼
$ fastqc -o your/path/to/out -t 4 ${datadir}/*.fq.gz
#是不是超級(jí)簡(jiǎn)單，-o是報(bào)告輸出的目錄桥氏，這個(gè)目錄必須是存在的温峭，-t是線程數(shù)
- 結(jié)果如下，要放在早幾年字支，估計(jì)是很難見到這么干凈的原始數(shù)據(jù)了凤藏，為此我還專門請(qǐng)教了一些老師，這么干凈的數(shù)據(jù)的可信程度堕伪。得益于illumina公司對(duì)儀器的升級(jí)改善揖庄，數(shù)據(jù)的質(zhì)量也是越來越好了。
- 但是即便如此欠雌，我們還是要對(duì)raw data進(jìn)行清洗的蹄梢，數(shù)據(jù)的質(zhì)量直接決定了后面分析的準(zhǔn)確性，是一切分析的前提桨昙，我們一直都在說QC检号，QC也是基礎(chǔ)，但是它的重要性不容忽視蛙酪。下面是我要查看的質(zhì)控的內(nèi)容：
- read各個(gè)位置的堿基質(zhì)量值分布
- 堿基的總體質(zhì)量值分布
- read各個(gè)位置上堿基分布比例齐苛，目的是為了分析堿基的分離程度
- GC含量分布
- read各位置的N含量
- read是否還包含測(cè)序的接頭序列
- 下面是我的對(duì)于clean data的指標(biāo)：
- 將含有接頭的reads對(duì)去除，或者cut接頭
- 測(cè)序錯(cuò)誤率分布檢查桂塞，一般情況下凹蜂，每個(gè)堿基位置的測(cè)序錯(cuò)誤率都應(yīng)該低于1%
- GC含量，理論上A=T阁危，C=G玛痊，前幾個(gè)堿基可能會(huì)有波動(dòng)，GC含量整體在41.5%左右
- 平均質(zhì)量值分布：Q30平均比例≥80%狂打，平均錯(cuò)誤率≤0.1%
- 為了達(dá)到這樣的目標(biāo)擂煞，我的篩選是條件是：
- 去除含接頭（adapter）的一對(duì)reads；可以酌情考慮去除接頭保留reads
- 當(dāng)某一端read中的N含量超過該條read長(zhǎng)度比例10%趴乡，去除該對(duì)reads
- 當(dāng)某一端read中低質(zhì)量（Q≤20）堿基數(shù)超過該read長(zhǎng)度比例的50%時(shí)对省，去除該對(duì)reads
$ fastp -i 1.fq.gz -I 2.fq.gz \
--adapter_sequence ${Adapter_R1} \
--adapter_sequence_r2 ${Adapter_R2} \#不同的試劑盒、儀器晾捏，會(huì)有不同的接頭蒿涎，這個(gè)可以咨詢公司，也可以選擇去掉這兩個(gè)參數(shù)惦辛，影響不大
-o your/path/of/data/cleaned.1.fq.gz \
-O your/path/of/data/cleaned.2.fq.gz \#輸出的clean data的位置和名稱
-c -q 20 -u 50 -n 15 -5 20 -3 20 -w 16 \
-h your/path/of/report/clean.html -j your/path/of/report/clean.json #這倆參數(shù)分別是不同格式報(bào)告的輸出位置和文件名
- 別的不多說劳秋，解釋一下我用到的幾個(gè)參數(shù)，
- -c：對(duì)overlap區(qū)域進(jìn)行糾錯(cuò)，適用于PE
- -q：設(shè)置低質(zhì)量的標(biāo)準(zhǔn)玻淑，默認(rèn)是15嗽冒，多數(shù)公司這里設(shè)為5
- -u：低質(zhì)量堿基所占比例，默認(rèn)40代表40%岁忘，只要有一條read不滿足條件就成對(duì)丟掉
- -n：過濾N堿基過多的reads辛慰，15代表個(gè)數(shù)，因?yàn)橐话鉖E150的reads長(zhǎng)度是150
- -w：線程數(shù)干像，這里要注意帅腌！范圍是1-16，建議是8麻汰，親測(cè)速客，8和16差不多。
- -5 -3：根據(jù)質(zhì)量值來截取reads五鲫，分別對(duì)應(yīng)5‘端和3’端溺职，得到reads長(zhǎng)度可能不等
- 具體請(qǐng)參考官網(wǎng)說明。
- 這個(gè)時(shí)候再看clean data結(jié)果位喂，就無需fastqc了浪耘，因?yàn)閒astp也會(huì)生成一份報(bào)告。
- 看到結(jié)果塑崖，其實(shí)可以發(fā)現(xiàn)七冲，數(shù)據(jù)還是保留了大部分的，所以完全不用擔(dān)心cut太狠规婆，導(dǎo)致數(shù)據(jù)不夠澜躺。我的原則是，在數(shù)據(jù)量面前抒蚜，更應(yīng)該保證的是數(shù)據(jù)的準(zhǔn)確性掘鄙，即使會(huì)有一點(diǎn)損失，也是值得的嗡髓，畢竟research更注重的也是準(zhǔn)確性操漠，相對(duì)于降低假陰性，我更傾向于提高真陽性饿这。
- 那么得到clean data以后就可以進(jìn)行下一步mapping了浊伙，mapping內(nèi)容下次再與大家分享。
- 水平有限蛹稍，要是存在什么錯(cuò)誤請(qǐng)?jiān)u論指出吧黄！請(qǐng)大家多多批評(píng)指正部服，相互交流唆姐，共同成長(zhǎng)，謝謝＠恕７盥赵抢！
