寫在前面

你現在看到的是文獻俱樂部2019年筆記分享第一彈,，我將會在春節(jié)7天連續(xù)分享，目錄如下：

• 2019年1月份第1周（總第49周）單細胞轉錄組探索CAFs的功能和空間異質性

• 2019年1月份第2周（總第50周）異常CRC病人的突變時空異質性與免疫

• 2019年1月份第3周（總第51周）探索PDAC癌前病變

• 2019年1月份第4周（總第52周）TCGA計劃的ATAC-seq數據發(fā)布

• 2019年2月份第1周（總第53周）胃癌類器官

• 2019年2月份第2周（總第54周）測173個成年人的大腦的102個基因

• 2019年2月份第3周（總第55周）2.5萬漢族人的GWAS乳腺癌風險基因

因為學業(yè)需要,，我閱讀的大量文獻都是NGS組學相關,，所以會涉及到很多數據處理，而這些文獻基于的生物信息學數據處理技巧,，我都在過去的5年里以各種形式分享講解過,，也有系列視頻，希望你可以在方便的時候再次學習一遍,，查漏補缺。也歡迎推薦給有需要的朋友

如果,，你不僅僅是對NGS組學應用文獻感興趣,，也歡迎加入我們文獻閱讀小組分享自己的主頁領域文獻。

逆向收費讀文獻社群（第二年通知）

今天是大年初三,，給大家?guī)淼氖?strong>TCGA計劃的ATAC-seq數據發(fā)布,，希望你能學到知識。

今天的文獻解讀有點特殊,，雖然該文章在我2019年的48篇精品解讀文獻列表中,！

但是呢，我們技能樹的云曉早在年前就出了一個非常棒的解讀推文,，這里就直接摘抄她的,，借花獻佛分享給大家！

文獻概要

原文鏈接：The chromatin accessibility landscape of primary human cancers
http://science./content/362/6413/eaav1898
DOI: 10.1126/science.aav1898
期刊：Science

發(fā)表時間 26 October 2018

首先我們回顧了ATAC-seq方法的原理和優(yōu)點,，以及與其他研究染色質可及性方法的比較,，然后介紹了這篇文章的主要結果和亮點，以及提供的數據資源,。

導讀

染色質的可及性（chromatin accessibility）通常理解為開放染色質（open chromatin）,，指致密的核小體結構被破壞后，啟動子,、增強子,、絕緣子、沉默子等順式調控元件和反式作用因子可以接近的區(qū)域,，與真核生物的轉錄調控密切相關,。
目前研究染色質的可及性的方法有DNase-Seq，MNase-Seq,，FAIRE-seq和ATAC-seq,。ATAC-seq是2013年由斯坦福大學William J. Greenleaf和Howard Y. Chang實驗室開發(fā)的，通過Tn5轉座酶切割暴露的DNA并同時連接上特異性的adapters，然后連接上adapters的DNA片段被分離出來用于二代測序,。由于所需細胞量少,，實驗簡單，可以在全基因組范圍內檢測染色質的開放狀態(tài),，被廣泛使用,。

各種研究方法原理及優(yōu)缺點比較：

2018年10月26日ATAC-seq技術的兩個主要開發(fā)者William J. Greenleaf和Howard Y. Chang作為共同通訊在Science上又發(fā)表一重大成果——人類原發(fā)性腫瘤染色質可及性圖譜,。

這篇文章的主要結果包括以下幾方面：

• 繪制了來自410個腫瘤樣本橫跨23種癌癥類型的796個全基因組染色質可及性圖譜；

• 在這些癌癥類型中共發(fā)現了562,709個DNA調控元件,；

• 整合ATAC-seq與TCGA其他的多組學數據,，鑒定腫瘤特異的DNA調控元件，如遠端增強子具有更強的組織特異性,，根據增強子元件聚類鑒定到新的腫瘤亞型,；

• 通過TF足跡分析找到了關鍵的TF, 然后通過預測TF和DNA的相互作用模式以及基因的表達識別不同的TF活性；

• 基因表達和染色質可及性的關聯(lián)分析預測到大量遠端增強子與啟動子間的相互作用,，包括一些重要的致癌基因和腫瘤免疫治療的靶點,，如MYC,SRC, BCL2和PDL1，為免疫治療提供了新的視覺,；

• 另一個亮點是結合GWAS和WGS探索腫瘤變異的影響,。結果表明在調控元件處的變異通過產生或干擾轉錄因子的結合位點，可能增強或抑制染色質的可及性,。如位于12號染色體FGD4基因上游的單堿基突變,，會產生NKX 轉錄因子結合的基序，增強了染色質的可及性,，促進了FGD4基因的表達,。
  

  

• 另一個重要的亮點是這篇文章提供了豐富的數據的資源和腫瘤研究的一個新的視覺，但是他們做的是pan-cancer 研究,，只是做了初步探索,，具體到每個癌癥還有很多東西值得挖掘。

公布數據資源：

• 原始數據和三級數據都是開放的,，三級數據：https://gdc./about-data/publications/ATACseq-AWG,，包括counts,peaks,和bigwig文件等其他中間文件。
  

  

• ATAC-seq原始數據和bam文件： https://portal.gdc./

文章納入的23種腫瘤類型

如下：

主要結果解讀

1. 數據質控

數據分析前的質控,，TSS富集,、片段長度分布特征、peaks的reads比例,。
其中常用于判斷實驗是否失敗的片段長度分布特征原理如下（摘抄于ATAC-seq交流群中來自菲沙基因小伙伴的總結）：

核小體完整的話,，纏繞核小體的DNA不會被切,，被切的是核小體之間的片段，纏繞核小體的146bp的DNA是完整的,，加上接頭,，就會形成特有的280-300bp左右的片段。而且以這個為基數,，兩倍就是兩個核小體,，三倍就是三個核小體，依次類推,，就形成了層次分明的特征條帶,。如果核小體解聚，那么纏繞核小體的146bp也會被切除,，自然就沒有特定大小的DNA富集,，形成的就是一片模糊的彌散條帶。

• 這篇文章中他們首先做的是確保數據與捐贈者對應,。通過ATAC-seq的基因型與其對應的捐贈者的SNP芯片的基因型的相關性分析確保ATAC-seq的數據與捐贈者對應,。

• 然后也有核小體片段特征圖。片段大小的分布特征具有明顯的核小體周期性分布,，數據質量符合要求。
  

  

下面就是數據分析了：

2. 鑒定到多少DNA調控元件

call peaks 后對peaks數目統(tǒng)計,，當然對peaks數目統(tǒng)計前也會做質控,，只保留重復性好的peaks。 通過對410個腫瘤樣本的23種腫瘤類型（其中386個樣本有技術重復）進行ATAC-seq分析,，共繪制出796個染色質可及性圖譜,，鑒定到 562,709個調控元件，即ATAC-seq數據分析中對peaks數目統(tǒng)計,，共有562,709個可重復,、轉座酶可接近的染色質可及性位點。其中562,709peaks是總的數目,，實際每種癌癥類型的peaks數目從 56,125 到215,978,，數目不等。

• ATAC-seq peaks與Roadmap DNase peaks的overlapy以及與chromHMM-defined regulatory states的overlap
  另外，將ATAC-seq得到的腫瘤特異的peaks與Roadmap Epigenetics project中DNase-seq測序得到peaks比較,，該研究中的peaks與以往發(fā)現的調控元件共有65%的overlap,。該結果一表明了此研究中觀測到的調控元件與以往研究中的一致性，二揭示了腫瘤樣本中新的染色質可及性敏感位點,。
  

3. ATAC-seq遠端元件,，promotor，RNA-seq得到的表達矩陣聚類

對遠端元件和啟動子，還有RNA-Seq的表達量做Pearsn相關性層次聚類,，遠端元件展現出更好的腫瘤特異性,。

4. 聚類,、腫瘤亞型分型,、鑒定cluster-specific peaks

取所有腫瘤類型中變化大的250,000 peaks的top 50主成分進行t-SNE聚類，鑒定到18個模塊,，且發(fā)現基于ATAC-seq的聚類與用TCGA中mRNA-seq, miRNA–seq, DNA 甲基化 reverse-phase protein array (RPPA)和DNA拷貝數數據做iCluster的結果高度一致,。聚類的結果表明一些癌癥被分成兩個模塊，如乳腺癌分為基底和非基底的,，食管癌分為鱗癌和腺癌,；來自相同組織的腫瘤樣本經常聚在一起，如腎透明細胞癌和腎乳頭狀細胞癌;有的是不同組織相同類型的腫瘤聚在一起,，如鱗癌,。

5. cluster-specific peaks的遠端元件、motif和甲基化分析

遠端調控元件和TFs都展現出組織和腫瘤特異性,。

6. Footprinting分析TF的活動

轉錄因子對染色質可及性和腫瘤的產生和轉移有什么影響,？首先做的是TF的足跡分析。
理解幾個概念：

• TF footprint：當一個或多個核小體的移位時在側翼序列中會產生高DNA可及性,，而TF與DNA的結合會保護蛋白質-DNA結合位點免于轉位,。這種現象認為是TF footprint。

• Flanking Accessbility: a measure of the accessibility of the DNA adjacent to a TF motif,；

• Footprint Depth: a measure of the relative protection of the motif site from transposition
  這篇文章用了兩種TF 足跡的分析方法：

• 1） quantifies the “flanking accessibility（doi: 10.1016/j.celrep.2017.05.003）,；

• 2） ChromVAR：（ doi：10.1038/nmeth.4401）
  兩種方法的得到TFs高度重合。
  

  

  
  結論：

• 一個TF如果能夠與DNA結合,，那么footprint depth 和 flanking accessibility 與其基因表達顯著相關

• flanking accessibility的增高和footprint depth的降低可能伴隨著甲基化水平的降低
  

  

7. DNA調控元件與基因的相互作用

為了識別鑒定ATAC-seq peaks和基因表達之間的假定因果聯(lián)系,，他們基于相關性的方法建立模型進行預測。具體方法如下圖所示：

8. 突變是如何影響染色質可及性和腫瘤的發(fā)生,？

這篇文章其中的一個亮點是揭示了調控元件處的突變影響著染色質的可及性,，如突變通過影響轉錄因子的結合調控染色質的可及性。他們通過整合WGS數據和ATAC-seq數據,，鑒定到上千個體細胞突變落在啟動子區(qū)域和調控元件區(qū)域,。如TERT啟動子區(qū)域的突變，這也表明ATAC-seq與外顯子測序相比,，可以鑒定調控元件處的突變,，而WES是不包括啟動子區(qū)域,。

除了TERT啟動子處的突變，還有增強子的突變,，如FGD4,。eFGD4的突變產生NKK TF的結合位點，導致在FGD4上游12-kb區(qū)域內的染色質可及性大幅增加,。FGD4的高表達與膀胱癌低生存率顯著相關,，進一步表明FGD4的突變可能對特定癌癥產生影響。

除了體細胞突變,，他們還整合了GWAS數據,，發(fā)現MYC位點附近有兩個腫瘤易感性SNP位點rs6983267和rs35252396。其中SNP rs6983267與結腸癌和前列腺癌的染色質可及性增強相關,，該結果與以往報道一致,，同時也與乳腺癌和鱗癌染色質可及性相關，后者是之前研究沒有發(fā)現的,。另一個SNP位點rs35252396與腎透明細胞癌（KIRC）,、乳腺癌、甲狀腺癌的染色質可及性有強烈的相關性,。這些結果都表明SNP可能在腫瘤中扮演者潛在的重要作用,。

9. 鑒定與腫瘤免疫治療相關的DNA調控元件

他們基于已知的免疫細胞特異性調節(jié)元件的可及性來估計免疫浸潤的水平。這些區(qū)域的可及性還與腫瘤純度成反比,，為去卷積體瘤數據提供了附加信息,。
同時他們還研究了與免疫治療的一個重要靶標PDL1相關的峰值，PDL1的表達受50kb之內的四種調控元件的影響,，并且不同癌癥類型的可及性不同,。了解PDL1和其他藥物靶點的調節(jié)元件的狀態(tài)或許可以為個性化治療提供一個路徑,。

數據處理方法詳解

1. ATAC-seq數據預處理和比對

ATAC-seq預處理和比對使用的是PEPATAC pipeline(http://code./PEPATAC/),。
PEPATAC pipeline是一個打包的ATAC-seq數據預處理流程，包括對原始數據的去接頭,、比對,、call peak、創(chuàng)建bigwig,、TSS富集文件等其他一些統(tǒng)計結果文件,。
輸出的圖如：

具體包括：

• Bowtie2 比對，移除比對到chrM和重復序列

• 排序去除重復
  使用Picard 的MarkDuplicates去除重復,。

2. call peaks(MACS2)

這里他們選用固定寬度（fixed-width）的peaks,優(yōu)點有：1）對大量的peaks進行counts和motif分析時可以減小誤差,；2）對于大量數據集的可以合并峰得到一致性的peaks; 
使用的是macs2 call peaks,參數如下：

同時根據hg38 blacklist過濾，并除去染色體兩端以外的峰,。
一個樣本的overlaps他們是通過迭代移除的方法,，首先保留最顯著的peak,然后任何與最顯著peak有直接overlap的peaks都被移除,；接著對另一個最顯著性的peak進行相同的操作，最終保留所有更顯著的peaks,，移除與其有直接overlaps的peaks,。

3. ATAC-seq數據分析—— 構建counts矩陣并標準化

為了獲得每個峰中獨立的Tn5插入的數量，首先用RRsamtools “scanbam”對BAM文件矯正Tn5偏移量（“+” stranded +4 bp, “-” stranded -5 bp）并存入Genomic Ranges對象,。然后用“countOverlaps”對矯正后的插入位點計數,，最終得到 562,709 x 796 counts 矩陣。
counts矩陣用edgeR “cpm(matrix , log = TRUE,prior.count = 5)”標準化,，然后用R中的preprocessCore’s “normalize.quantiles”做分位數標準化,。

4. ATAC-seq data analysis – Transcription factor footprinting

TF足跡的分析：
一是參考了文章doi: 10.1016/j.celrep.2017.05.003：

• 首先確定peaks內的TF motif的位置，用pan-cancer peak set 結合CIS-BP motifs計算motif的位置,，motifmatchr “matchMotifs(positions = “out”)

• 然后計算flanking accessibility 和 footprint depth

• 最后確定哪個TF的足跡與基因的表達是顯著相關
  通過將flanking accessibility or footprint depth與250個隨機的TFs的關聯(lián)分析生成零均值和標準偏差,。

5. ATAC-seq data analysis – chromVAR for transcription factor activity

除了足跡分析，他們還用chromVAR包評估TF的活動,，首先用chromVAR deviations函數計算GC矯正偏差,，然后將矯正偏差與motif相關的TFs關聯(lián)，最后5000個轉錄因子基序和非相關轉錄因子基因的RNA-seq基因表達之間的隨機相關性,，以計算每個相關性的FDR,。具體參考：Week4— chromVAR:預測染色質可及性相關的轉錄因子

6. ATAC-seq data analysis – chromVAR for GWAS enrichment

• 首先從GWAS catalog（https://www./gwas/docs/file-downloads）下載SNPs位點，過濾和16種癌癥類型相關的SNPs位點,。

• 加上連鎖不平衡（Linkage Disequilibrium ,，LD) 信息（ r 2 > 0.8）
  LD信息從haploreg 網站下載 http://archive./mammals/haploreg/data/

• 移走位于exons或UTR區(qū)域的SNPs位點，得到最后的SNP列表

• 將最后的SNP列表與遠端 binarization peak 集overlap,，得到一個二元匹配矩陣,。每列代表不同癌癥癌癥類型的GWAS SNP，每行代表一個peak,，這個peak來自遠端 binarization peak 集,。

• 用chromVAR deviations函數計算GC矯正偏差

• 用PNAMER將“偏差分數”轉換為p值，并使用Bejimi-HocHBG程序調整

后記

如果你完全沒有看懂文章說了些什么,，卻仍然堅持到了最后,，說明你有可能對生物信息學感興趣，你缺乏的是一個入門的契機,！

也歡迎推薦你身邊的朋友參與我們的線下培訓,，如果有緣的話