TCGA计划的ATAC-seq数据发布 / 四六文摘

写在前面

你现在看到的是文献俱乐部2019年笔记分享第一弹，我将会在春节7天连续分享，目录如下：

2019年1月份第1周（总第49周）单细胞转录组探索CAFs的功能和空间异质性
2019年1月份第2周（总第50周）异常CRC病人的突变时空异质性与免疫
2019年1月份第3周（总第51周）探索PDAC癌前病变
2019年1月份第4周（总第52周）TCGA计划的ATAC-seq数据发布
2019年2月份第1周（总第53周）胃癌类器官
2019年2月份第2周（总第54周）测173个成年人的大脑的102个基因
2019年2月份第3周（总第55周）2.5万汉族人的GWAS乳腺癌风险基因

因为学业需要，我阅读的大量文献都是NGS组学相关，所以会涉及到很多数据处理，而这些文献基于的生物信息学数据处理技巧，我都在过去的5年里以各种形式分享讲解过，也有系列视频，希望你可以在方便的时候再次学习一遍，查漏补缺。也欢迎推荐给有需要的朋友

学习笔记目录

1.3个学生的linux视频学习笔记

2.生信人应该这样学R语言系列视频学习心得笔记分享

3.一万人陪你学习GEO数据库挖掘知识（公益视频听课笔4.记分享）

5.公共数据库挖掘视频学习心得体会

6.生信小技巧系列第一季完结版视频教程学习笔记分享

……期待有你……

如果，你不仅仅是对NGS组学应用文献感兴趣，也欢迎加入我们文献阅读小组分享自己的主页领域文献。

逆向收费读文献社群（第二年通知）

今天是大年初三，给大家带来的是TCGA计划的ATAC-seq数据发布，希望你能学到知识。

今天的文献解读有点特殊，虽然该文章在我2019年的48篇精品解读文献列表中！

但是呢，我们技能树的云晓早在年前就出了一个非常棒的解读推文，这里就直接摘抄她的，借花献佛分享给大家！

文献概要

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

发表时间 26 October 2018

首先我们回顾了ATAC-seq方法的原理和优点，以及与其他研究染色质可及性方法的比较，然后介绍了这篇文章的主要结果和亮点，以及提供的数据资源。

导读

染色质的可及性（chromatin accessibility）通常理解为开放染色质（open chromatin），指致密的核小体结构被破坏后，启动子、增强子、绝缘子、沉默子等顺式调控元件和反式作用因子可以接近的区域，与真核生物的转录调控密切相关。
目前研究染色质的可及性的方法有DNase-Seq，MNase-Seq，FAIRE-seq和ATAC-seq。ATAC-seq是2013年由斯坦福大学William J. Greenleaf和Howard Y. Chang实验室开发的，通过Tn5转座酶切割暴露的DNA并同时连接上特异性的adapters，然后连接上adapters的DNA片段被分离出来用于二代测序。由于所需细胞量少，实验简单，可以在全基因组范围内检测染色质的开放状态，被广泛使用。

https://doi.org/10.1186/1756-8935-7-33；doi: 10.1038/nmeth.2688

各种研究方法原理及优缺点比较：

2018年10月26日ATAC-seq技术的两个主要开发者William J. Greenleaf和Howard Y. Chang作为共同通讯在Science上又发表一重大成果——人类原发性肿瘤染色质可及性图谱。

这篇文章的主要结果包括以下几方面：

绘制了来自410个肿瘤样本横跨23种癌症类型的796个全基因组染色质可及性图谱；
在这些癌症类型中共发现了562,709个DNA调控元件；
整合ATAC-seq与TCGA其他的多组学数据，鉴定肿瘤特异的DNA调控元件，如远端增强子具有更强的组织特异性，根据增强子元件聚类鉴定到新的肿瘤亚型；
通过TF足迹分析找到了关键的TF, 然后通过预测TF和DNA的相互作用模式以及基因的表达识别不同的TF活性；
基因表达和染色质可及性的关联分析预测到大量远端增强子与启动子间的相互作用，包括一些重要的致癌基因和肿瘤免疫治疗的靶点，如MYC,SRC, BCL2和PDL1，为免疫治疗提供了新的视觉；
另一个亮点是结合GWAS和WGS探索肿瘤变异的影响。结果表明在调控元件处的变异通过产生或干扰转录因子的结合位点，可能增强或抑制染色质的可及性。如位于12号染色体FGD4基因上游的单碱基突变，会产生NKX 转录因子结合的基序，增强了染色质的可及性，促进了FGD4基因的表达。

Cancer gene regulatory landscape
另一个重要的亮点是这篇文章提供了丰富的数据的资源和肿瘤研究的一个新的视觉，但是他们做的是pan-cancer 研究，只是做了初步探索，具体到每个癌症还有很多东西值得挖掘。

公布数据资源：

原始数据和三级数据都是开放的，三级数据：https://gdc.cancer.gov/about-data/publications/ATACseq-AWG，包括counts,peaks,和bigwig文件等其他中间文件。
ATAC-seq原始数据和bam文件： https://portal.gdc.cancer.gov/

文章纳入的23种肿瘤类型

如下：

ACC, adrenocortical carcinoma肾上腺皮质癌; BLCA, bladder urothelial carcinoma 膀胱上皮癌; BRCA, breast invasive carcinoma乳腺浸润性癌; CESC, cervical squamous cell carcinoma宫颈鳞癌; CHOL, cholangio carcinoma胆管癌; COAD, colon adenocarcinoma结肠癌; ESCA, esophageal carcinoma食道癌; GBM, glioblastoma multiforme胶质母细胞瘤; HNSC, head and neck squamous cell carcinoma头颈部鳞状细胞癌; KIRC, kidney renal clear cell carcinoma肾透明细胞癌; KIRP, kidney renal papillary cell carcinoma肾乳头状细胞癌; LGG, low grade glioma低级别胶质瘤; LIHC, liver hepatocellular carcinoma肝癌; LUAD, lung adenocarcinoma肺腺癌; LUSC, lung squamous cell carcinoma肺鳞状细胞癌; MESO, mesothelioma间皮细胞瘤; PCPG, pheochromocytoma and paraganglioma嗜铬细胞瘤和副神经节瘤; PRAD, prostate adenocarcinoma前列腺癌; SKCM, skin cutaneous melanoma皮肤黑色素瘤; STAD, stomach adenocarcinoma胃腺癌; TGCT, testicular germ cell tumors睾丸肿瘤; THCA, thyroid carcinoma; 甲状腺癌 UCEC, uterine corpus endometrial carcinoma子宫内膜癌.

主要结果解读

1. 数据质控

数据分析前的质控，TSS富集、片段长度分布特征、peaks的reads比例。
其中常用于判断实验是否失败的片段长度分布特征原理如下（摘抄于ATAC-seq交流群中来自菲沙基因小伙伴的总结）：

核小体完整的话，缠绕核小体的DNA不会被切，被切的是核小体之间的片段，缠绕核小体的146bp的DNA是完整的，加上接头，就会形成特有的280-300bp左右的片段。而且以这个为基数，两倍就是两个核小体，三倍就是三个核小体，依次类推，就形成了层次分明的特征条带。如果核小体解聚，那么缠绕核小体的146bp也会被切除，自然就没有特定大小的DNA富集，形成的就是一片模糊的弥散条带。

这篇文章中他们首先做的是确保数据与捐赠者对应。通过ATAC-seq的基因型与其对应的捐赠者的SNP芯片的基因型的相关性分析确保ATAC-seq的数据与捐赠者对应。
然后也有核小体片段特征图。片段大小的分布特征具有明显的核小体周期性分布，数据质量符合要求。

Pan-cancer ATAC-seq data from frozen tissues is high quality and internally robust

下面就是数据分析了：

2. 鉴定到多少DNA调控元件

call peaks 后对peaks数目统计，当然对peaks数目统计前也会做质控，只保留重复性好的peaks。通过对410个肿瘤样本的23种肿瘤类型（其中386个样本有技术重复）进行ATAC-seq分析，共绘制出796个染色质可及性图谱，鉴定到 562,709个调控元件，即ATAC-seq数据分析中对peaks数目统计，共有562,709个可重复、转座酶可接近的染色质可及性位点。其中562,709peaks是总的数目，实际每种癌症类型的peaks数目从 56,125 到215,978，数目不等。

肿瘤类型、样本数量和peaks的基因组分布特征和数量

ATAC-seq peaks与Roadmap DNase peaks的overlapy以及与chromHMM-defined regulatory states的overlap
另外，将ATAC-seq得到的肿瘤特异的peaks与Roadmap Epigenetics project中DNase-seq测序得到peaks比较，该研究中的peaks与以往发现的调控元件共有65%的overlap。该结果一表明了此研究中观测到的调控元件与以往研究中的一致性，二揭示了肿瘤样本中新的染色质可及性敏感位点。

3. ATAC-seq远端元件，promotor，RNA-seq得到的表达矩阵聚类

对远端元件和启动子，还有RNA-Seq的表达量做Pearsn相关性层次聚类，远端元件展现出更好的肿瘤特异性。

远端元件、启动子、RNA-seq相关性聚类比较

4. 聚类、肿瘤亚型分型、鉴定cluster-specific peaks

取所有肿瘤类型中变化大的250,000 peaks的top 50主成分进行t-SNE聚类，鉴定到18个模块，且发现基于ATAC-seq的聚类与用TCGA中mRNA-seq, miRNA–seq, DNA 甲基化 reverse-phase protein array (RPPA)和DNA拷贝数数据做iCluster的结果高度一致。聚类的结果表明一些癌症被分成两个模块，如乳腺癌分为基底和非基底的，食管癌分为鳞癌和腺癌；来自相同组织的肿瘤样本经常聚在一起，如肾透明细胞癌和肾乳头状细胞癌;有的是不同组织相同类型的肿瘤聚在一起，如鳞癌。

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

远端调控元件和TFs都展现出组织和肿瘤特异性。

cluster-specific peaks相关的远端元件和TF

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和基因表达之间的假定因果联系，他们基于相关性的方法建立模型进行预测。具体方法如下图所示：

8. 突变是如何影响染色质可及性和肿瘤的发生？

这篇文章其中的一个亮点是揭示了调控元件处的突变影响着染色质的可及性，如突变通过影响转录因子的结合调控染色质的可及性。他们通过整合WGS数据和ATAC-seq数据，鉴定到上千个体细胞突变落在启动子区域和调控元件区域。如TERT启动子区域的突变，这也表明ATAC-seq与外显子测序相比，可以鉴定调控元件处的突变，而WES是不包括启动子区域。

WGS与ATAC的联合分析

除了TERT启动子处的突变，还有增强子的突变，如FGD4。eFGD4的突变产生NKK TF的结合位点，导致在FGD4上游12-kb区域内的染色质可及性大幅增加。FGD4的高表达与膀胱癌低生存率显著相关，进一步表明FGD4的突变可能对特定癌症产生影响。

WGS与ATAC的联合分析

除了体细胞突变，他们还整合了GWAS数据，发现MYC位点附近有两个肿瘤易感性SNP位点rs6983267和rs35252396。其中SNP rs6983267与结肠癌和前列腺癌的染色质可及性增强相关，该结果与以往报道一致，同时也与乳腺癌和鳞癌染色质可及性相关，后者是之前研究没有发现的。另一个SNP位点rs35252396与肾透明细胞癌（KIRC）、乳腺癌、甲状腺癌的染色质可及性有强烈的相关性。这些结果都表明SNP可能在肿瘤中扮演者潜在的重要作用。

GWAS和ATAC-seq联合分析

9. 鉴定与肿瘤免疫治疗相关的DNA调控元件

他们基于已知的免疫细胞特异性调节元件的可及性来估计免疫浸润的水平。这些区域的可及性还与肿瘤纯度成反比，为去卷积体瘤数据提供了附加信息。
同时他们还研究了与免疫治疗的一个重要靶标PDL1相关的峰值，PDL1的表达受50kb之内的四种调控元件的影响，并且不同癌症类型的可及性不同。了解PDL1和其他药物靶点的调节元件的状态或许可以为个性化治疗提供一个路径。

数据处理方法详解

1. ATAC-seq数据预处理和比对

ATAC-seq预处理和比对使用的是PEPATAC pipeline(http://code.databio.org/PEPATAC/)。
PEPATAC pipeline是一个打包的ATAC-seq数据预处理流程，包括对原始数据的去接头、比对、call peak、创建bigwig、TSS富集文件等其他一些统计结果文件。
输出的图如：

具体包括：

Bowtie2 比对，移除比对到chrM和重复序列

-k 1 -D 20 -R 3 -N 1 -L 20 -i S,1,0.50 -X 2000 –rg-id # remove repeats的参数 --very-sensitive -X 2000 --rg-id # bowtie2参数

排序去除重复
使用Picard 的MarkDuplicates去除重复。

-f 2 -q 10 -b -@ 20 # 排序参数 VALIDATION_STRINGENCY =LENIENT REMOVE_DUPLICATES = true #去重参数

2. call peaks(MACS2)

这里他们选用固定宽度（fixed-width）的peaks,优点有：1）对大量的peaks进行counts和motif分析时可以减小误差；2）对于大量数据集的可以合并峰得到一致性的peaks;
使用的是macs2 call peaks,参数如下：

--shift -75 --extsize 150 --nomodel --call-summits --nolambda --keep-dup all -p 0.01

同时根据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的关联分析生成零均值和标准偏差。

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

除了足迹分析，他们还用chromVAR包评估TF的活动，首先用chromVAR deviations函数计算GC矫正偏差，然后将矫正偏差与motif相关的TFs关联，最后5000个转录因子基序和非相关转录因子基因的RNA-seq基因表达之间的随机相关性，以计算每个相关性的FDR。具体参考：Week4— chromVAR:预测染色质可及性相关的转录因子

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

首先从GWAS catalog（https://www.ebi.ac.uk/gwas/docs/file-downloads）下载SNPs位点，过滤和16种癌症类型相关的SNPs位点。
加上连锁不平衡（Linkage Disequilibrium ，LD) 信息（ r 2 > 0.8）
LD信息从haploreg 网站下载 http://archive.broadinstitute.org/mammals/haploreg/data/
移走位于exons或UTR区域的SNPs位点，得到最后的SNP列表
将最后的SNP列表与远端 binarization peak 集overlap，得到一个二元匹配矩阵。每列代表不同癌症癌症类型的GWAS SNP，每行代表一个peak，这个peak来自远端 binarization peak 集。
用chromVAR deviations函数计算GC矫正偏差
用PNAMER将“偏差分数”转换为p值，并使用Bejimi-HocHBG程序调整

后记

如果你完全没有看懂文章说了些什么，却仍然坚持到了最后，说明你有可能对生物信息学感兴趣，你缺乏的是一个入门的契机！

也欢迎推荐你身边的朋友参与我们的线下培训，如果有缘的话

TCGA计划的ATAC-seq数据发布