第三章 研究方法
3.4 實驗流程
圖 3.2 實驗流程圖
實驗的流程如圖 3.2 所示，首先我們將從 ENCODE 的資料庫中，下載我們 欲分析的 ChIP-seq 資料，選擇了 2014 年後且為成年人的 ChIP-seq 資料如下圖 3.3 所示。
圖 3.3 ENCODE 資料庫下載資料頁面
ARID3A
ChIP-seq of GM12878 (Homo sapiens, adult)
ATF7
ChIP-seq of K562 (Homo sapiens, adult 53 year)
BACH1
ChIP-seq of A549 (Homo sapiens, adult 58 year)ChIP-seq of GM12878 (Homo sapiens, adult)
BHLHE40 ChIP-seq of GM12878 (Homo sapiens,
adult)
BMI1
ChIP-seq of K562 (Homo sapiens, adult 53 year)ChIP-seq of K562 (Homo sapiens, adult 53 year)
BRCA1
ChIP-seq of A549 (Homo sapiens, adult17
58 year)
CBX5
ChIP-seq of K562 (Homo sapiens, adult 53 year)CEBPB
ChIP-seq of GM12878 (Homo sapiens, adult)CHD1
ChIP-seq of A549 (Homo sapiens, adult 58 year)ChIP-seq of HeLa-S3 (Homo sapiens, adult 31 year)
CREB3L1 ChIP-seq of K562 (Homo sapiens, adult
53 year)
CREM
ChIP-seq of GM12878 (Homo sapiens, adult)ChIP-seq of K562 (Homo sapiens, adult 53 year)
表 3-1 轉錄因子與其對應的細胞株
結果如表 3-1 所示，我們可以看到不同細胞株探討的轉錄因子，統計後得到 2014 年 ENCODE 資料庫成人 ChIP-seq 資料總共有 68 個轉錄因子和 20 個細胞 株，且各個細胞株所探討的疾病分別為下表 3-2。
細胞株 相關疾病
ChIP-seq of A549 (Homo sapiens, adult 58 year)
肺癌ChIP-seq of GM12878
皰疹病毒ChIP-seq of T47D (Homo sapiens, adult 54 year)
乳癌ChIP-seq of NB4 (Homo sapiens, adult 23 year)
白血病ChIP-seq of MCF-7 (Homo sapiens, adult 69 year)
乳腺癌ChIP-seq of HeLa-S3 (Homo sapiens, adult 31 year) 子宮頸癌
ChIP-seq of Oci-Ly-3 (Homo sapiens, adult 52 year) 淋巴癌
ChIP-seq of GM12878 (Homo sapiens, adult)
皰疹病毒ChIP-seq of Ishikawa (Homo sapiens, adult)
子宮內膜腺癌ChIP-seq of LNCaP clone FGC (Homo sapiens,
前列腺癌adult 50 year)
ChIP-seq of DOHH2 (Homo sapiens, adult 60 year)
淋巴癌ChIP-seq of Caco-2 (Homo sapiens, adult 72 year)
大腸腺癌ChIP-seq of Oci-Ly-7 (Homo sapiens, adult 48 year) 淋巴癌 ChIP-seq of SUDHL6 (Homo sapiens, adult 43
year)
淋巴癌
ChIP-seq of H54 (Homo sapiens, adult 36 year)
多形性膠質母 細胞瘤ChIP-seq of Karpas-422 (Homo sapiens, adult 73
year)
淋巴癌
ChIP-seq of K562 (Homo sapiens, adult 53 year)
髓細胞性白血 病ChIP-seq of HL-60 (Homo sapiens, adult 36 year)
骨髓細胞白血 病ChIP-seq of U-87 MG (Homo sapiens, adult 44 year)
多形性膠質母 細胞瘤
ChIP-seq of Ishikawa (Homo sapiens, adult 39
year)
NR3C1(Glucocorticoid receptor)、SP(Transcription factor Sp1)、MAFK(Transcription factor MafK)，分別對於這些轉錄因子進行分析。19
轉錄因子 細胞株
CEBPB
A549 Ishikawa HeLa-S3 MCF-7 K562MAFK
A549 MCF-7NR3C1
A549 Ishikawa GM12878SP
A549 K562表 3-3 為本實驗使用的細胞株和轉錄因子的免疫沉澱資料
3.4.1 序列回貼
從 ENCODE 資料庫，下載上述所說的細胞株與其轉錄因子的資料，這些從 ENCODE 資料庫下載的資料檔案，都為 fasta 格式，如圖，我們可以看到 Chip-seq 資料，是很多的 reads 合成的檔案，接下來就把這些資料，回貼到人類參考基因
圖 3.4 為 MACS 輸出的 bed 檔格式
圖 3.4，為一 MACS 輸出檔案，可以看到輸出檔的資訊，其格式每一列都是一個 峰值，第一欄是這個峰值出現在參考基因組的第幾條染色體，二三欄分別是這個 峰值出現的起始位和終點，第四欄是這個峰值的總長度是多少，當然長度越長的 峰值比較有可能包含到特定蛋白的結合位，五六欄分別為峰值的高度和多少 reads 被貼到這個片段，後面就是一些分數的檢定，也因為 MACS 所釋出的峰值 數目很多，因此我們要在這些輸出檔裡做一些篩選，在此有興趣的是最後一行的 FDR 值，因為 FDR 值為找尋波峰時，對於這個波峰的真實性為何，FDR 實際就 是 false discovery rate 的縮寫，所以我們設定一些 FDR 值，並且觀察這些 FDR 值下的波峰數目，以及準確性的關係。
21
3.4.3 篩選 FDR
為了要選擇不同 FDR 值下的峰值的數目以及其準確率，我們拿了
TRANSFAC 資料庫的資訊，來幫助整個流程的進行，先以 CEBPB 這個轉錄因子 的分析為例，TRANSFAC 資料庫有的資訊是，已經被實驗驗證過的轉錄因子結 合位點位，和其轉錄因子為何，目前被收錄在 TRANSFAC 資料庫中的位置，總 共有 17441，所以我們把這些資訊，先把他們配對起來，把轉錄因子他自己的 ID 與名字，對到屬於他自己的位置上，如下圖 3.5。
圖 3.5 TRANSFAC 資料庫紀錄的位置與其對應的基因
我們可以比較容易看出 TRANSFAC 資料庫的資訊，第二行是代表這個位置 的資訊是出現在哪一條染色體上，二三行分別是該點位在其染色體出現的起點和 終點，後面是代表在這個位置上，有什麼基因出現在這個位置。
接下來我們要統計在 CEBPB 轉錄因子在 TRANSFAC 資料庫所收入的位置 總共有多少，整理完後得知 CEBPB 轉錄因子在 TRANSFAC 資料庫裡總共有 400 個位置，整理完的資訊如圖 3.6，
圖 3.6 CEBPB 在 TRANSFAC 資料庫被記載的位置
接下來我們把這 400 個位置與我們做出的 MACS 輸出檔進行比較，我們先 取了 5 個 FDR 值(%)，分別為 FDR(0%)、FDR(0.5%)、FDR(1%)、FDR(5%)和 FDR(無 篩選)，且在 ENCODE 得到的 ChIP-seq 資料中，擁有 CEBPB 這個轉錄因子的細 胞株有 ChIP-seq of A549 (Homo sapiens, adult 58 year)、ChIP-seq of K562 (Homo sapiens, adult 53 year)、ChIP-seq of HeLa-S3 (Homo sapiens, adult 31 year)、
ChIP-seq of MCF-7 (Homo sapiens, adult 69 year)和 ChIP-seq of Ishikawa (Homo sapiens, adult 39 year)這五個細胞株，所以這些 ChIP-seq 資料，分別都篩選了 FDR 值且與 TRANSFAC 資料庫做比較，我們要看這些 ChIP-seq 資料所得到的峰值，
是否出現在 TRANSFAC 所記錄的位置之中。
23
3.4.4 結合位特徵(Consensus, annotated motifs)
想要了解 ChIP-seq 峰值的可信程度，我們把 TRANSFAC 資料庫所記載的轉 錄因子的 consensus 拿出來，CEBPB 在 TRANSFAC 資料庫中紀載的 consensus ID 分別是 V$CEBPB_01、V$CEBPB_02、V$CEBPB_03、V$CEBPB_04、V$CEBPB_06 和 V$CEBPB_07，這些 consensus 都是由對應到的位置，透過這些位置的資訊建 立的 PWM 所得到。
圖 3.7，是 TRANSFAC 資料庫收入的 CEBPB 的 consensus，以 V$CEBPB_01 這個 ID 為例，他的 consensus 為 RNRTKDNGMAAKNN，意思是轉錄因子 CEBPB 結合到 DNA 上的序列為[AG][ATCG][AG]T[GT][AGT][ATCG]G[AC]AA[GT][ATCG][ATCG]
，這是經由統計結位置建立的 PFM 得到的 consensus，在這個 V$CEBPB_01ID 的 consensus 大小為 14 鹼基對(bp)，所以我們就把這樣的 consensus 從 MACS 輸 出的峰值中，對回去找這些峰值是否存在著這樣的 consensus，進一步來確定 ChIP-seq 資料的可信度，因此把 MACS 輸出的 bed 檔找到的峰值，寫程式去把 整個峰值的序列把它抽取出來。
圖 3.7 為 TRANSFAC 資料庫中紀錄的的 CEBPB 結合位特徵
圖 3.8 峰值被抽取出來的序列
這樣抽取出來的序列，是根據峰值在哪一條染色體的位置，從參考基因體 hg19 取出，所以圖 3.8，代表在參考基因體 hg19 的第一條染色體(chr1)的
11203(start)到 11620(end)的位置，然後我們將 TRANSFAC 資料庫記載的 CEBPB 轉錄因子的 consensus，寫程式從這些峰值的序列中，去尋找有沒有符合的峰值 和 pattern。
25
圖 3.9 結合位特徵找到的峰值與其找到的特徵
從圖 3.9，我們可以看到，在這些峰值的序列中，確實有發現符合 TRANSFAC 資料庫所記載的 consensus，我們找到的這些 motif 在正股和反股端都有出現，在 反股端出現的 motif 我們加了 rc 作為註解，因此我們把所有細胞株的 MACS 輸 出檔，和所篩選的 FDR 五個值都做了以上的步驟，然後使用 R 語言[16]繪製出 圖表，將在本論文的第四章進行討論。
3.4.5 TOP500 模序探勘 (De novo Motif Discovery)
這邊我們使用了 eTFBS[14]來做 De novo motif discovery，我們取了所有峰值 的前 500 名和後 500 名(用 P-value 排列)，來當正相關(positive)序列和負相關 (negative)序列，我們主要是想要在前 500 名的序列中找出共同存在的模序特徵，
後 500 名是來當作濾掉一些可能被找到不是正確的模序特徵，所以求出了 CEBPB 中五個細胞株(A549、HeLa-S3、Ishikawa、MCF-7、K562)的模序特徵和位置頻 率矩陣，所以我們將找到前三名的模序特徵對回輸出的總數峰值，重複之前的動 作來比較。