ChIP-seq/ATAC-seq peak calling and QC

最新推荐文章于 2024-09-30 17:19:22 发布
最新推荐文章于 2024-09-30 17:19:22 发布
阅读量444 收藏 22
点赞数 19
文章标签： 数据库 linux
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/telasi2004/article/details/142636074
版权

计算生物学实验4: ChIP-seq/ATAC-seq peak calling and QC

1. 实验目的

通过对ChIP-seq/ATAC-seq的数据进行peak calling、quality control（QC），熟悉常用的peak calling 工具MACS2、GEM。同时也要了解什么是blacklist region，并知道为什么、如何去除blacklist region。学习IGV browser的基本用法，学会用IGV查看bam、bdg(bedGraph)、narrowPeak文件，通过本实验，能够了解ChIP-seq/ATAC-seq的基本分析方法。

2. 实验准备

2.1 实验平台

Linux JSvr02 3.10.0-1160.80.1.el7.x86_64 #1 SMP Tue Nov 8 15:48:59 UTC 2022 x86_64 x86_64 x86_64 GNU/Linux

2.2 数据简述
数据数据来源
ATAC-seq ：SRR7696734NCBI数据库
ChIP-seq ：SRR5179211NCBI数据库
2.3 软件配置
Softwareversion
IGVIGV_2.6.3
FileZilla3.67.1_win64
MACS22.2.9.1
GEM3.4
bedtoolsbedtools2-2.29.2

3. 实验内容

分别用MACS2、GEM对CHIP/ATAC-seq数据进行peak calling，然后通过bedtools去除黑名单序列，通过IGV查看去除黑名单序列前后的区别，给出结论。

3.1 Prepare the working directory and visualize ChIP-seq and ATAC-seq data in IGV

下载数据、FIlezilla，IGV，将ChIP-seq/ATAC-seq对应的bam文件、bam.bai文件置于IGV browser中查看。
(a)Download files to your laptop
Use Filezilla to download:

  • ChIP-seq_SRR5179211_Rep1.bam
  • ChIP-seq_SRR5179211_Rep1.bam.bai
  • ATAC-seq_Rep1_SRR7696734.bam
  • ATAC-seq_Rep1_SRR7696734.bam.bai

bai文件是bam文件的索引文件。

(b)Load sacCer3(酵母菌) genome into IGV
选择参考基因组sacCer3。
(c )Open BAM files,view the following genomic coordinates:

chrV:175,594-181,774
在这里插入图片描述
chrXII:666,176-671,951
在这里插入图片描述
chrXIV:449,980-454,181
在这里插入图片描述

3.2 ChIP-seq peaks calling

(a) Call ChIP-seq peaks with MACS2

  1. 下载MACS2
conda install bioconda::macs2
  1. Call peak
macs2 callpeak -t ChIP-seq_SRR5179211_Rep1.bam -f BAM -g 1.2e7 -n CS_Rep1 -B -q \
0.05 --nomodel --extsize 100 --keep-dup all --call-summits
macs2 callpeak -t ChIP-seq_SRR5179212_Rep2.bam -f BAM -g 1.2e7 -n CS_Rep2 -B -q \
0.05 --nomodel --extsize 100 --keep-dup all --call-summits

MACS2全称Model-based Analysis of ChIP-Seq，是最常用的peak calling工具。

  • -t 实验组的比对结果；
  • -c 对照组（Control）的比对结果；
  • -f : 指定输入文件的格式。可以是“ELAND”, “BED”, “ELANDMULTI”, “ELANDEXPORT,“ELANDMULTIPET” (for pair-end tags), “SAM”, “BAM”, “BOWTIE”, “BAMPE” “BEDPE” 任意一个。如果不提供这项，则自动检测（对于Single-end数据），若是Paired-end的数据，则必须说明是“BAMPE”。如果使用了“BAMPE”，将跳过建立双峰，根据实际的插入大小来构建峰，同时参数–nomodel和–extsize失效；
  • -g：有效基因组大小。软件里给出默认的g 值：hs(人类): 2.7e9，mm(小鼠): 1.87e9，这里使用的是科学计数法，对于没有默认g值的需自己给出，比如本实验中的酵母菌，我们给出的是1.2e7；
  • -n：输出文件名的前缀；
  • -B/--bdg ：输出bedgraph格式的文件。
  • -q/--qvalue ：qvalue默认值是0.05，与pvalue不能同时使用，指定p值后MACS2就不会用q值了；
  • --nomodel : 这个参数说明不需要MACS去构建双峰模型，和extsize、shift是配套使用的，有这个参数才可以设置extsize和shift（MACS will bypass building the shifting model）；
  • --extsize : 当设置了nomodel时，MACS会用–extsize这个参数从5’->3’方向扩展reads修复fragments。比如说你的转录因子结合范围200bp，就设置这个参数是200。双端测序可以预测read具体长度，因此不适用此参数。
  • --keep-dup: 该选项控制在完全相同的位置的重复tags的行为。“auto”选项使MACS基于使用1e-5作为pvalue截止的二项式分布计算在完全相同的位置的maximum tags；“all”选项保留每个标签。如果给定一个整数，最多这个数量的标签将被保留在同一个位置，默认值为1。
  • --call-summits : MACS利用此参数重新分析信号谱，解析每个peak中包含的subpeak。

Output:
在这里插入图片描述

less CS_Rep1_peaks.narrowPeak

在这里插入图片描述
在narrowPeak文件中，各列的含义：

  • 1:染色体号
  • 2：peak起始位点
  • 3：结束位点
  • 4：peak name
  • 5：int(-10*log10qvalue)
  • 6 ：正负链
  • 7：fold change
  • 8：-log10pvalue
  • 9：-log10qvalue
  • 10: peak的中心，即summit距离peak起始位置的距离。

比如第二行的peak，第十列为70，起始位置（第二列）为12609，那么summit的位置大概在12680左右（对于这个例子，summit几乎位于peak的中心，但这只是一个特例）：

在这里插入图片描述
在IGV中查看，事实也确实如此。

  1. Filter peaks with blacklist file
    Remove peaks that overlap with our blacklist regions.(本实验中blacklist regions以bed文件的形式给出)
bedtools intersect -v -a CS_Rep1_peaks.narrowPeak -b sacCer3_blacklist.bed > \
CS_Rep1_peaks_filter.narrowPeak
bedtools intersect -v -a CS_Rep2_peaks.narrowPeak -b sacCer3_blacklist.bed > \
CS_Rep2_peaks_filter.narrowPeak

在这里插入图片描述

  • -v:输出在-a参数文件中没有overlap的区域（在此就是没有与blacklist region重叠的区域，即排除了所有与黑名单重叠的区域 ）；
  • -a, -b:输入文件（BAM/BED/GFF/VCF格式）。

比较去除blacklist前后的区别：

diff CS_Rep1_peaks.narrowPeak CS_Rep1_peaks_filter.narrowPeak

在这里插入图片描述

  • diff结果中出现的a=add，d=delete，c=change。

此处82,89d81的意思是第二个文件相比于第一个文件在第81行之后删除（delete）了第82-89行，202,219d193的意思是第二个文件相比于第一个文件在第193行后delete了第202-219行。（此处为什么数字不连贯了，因为193是相对于第二个文件，202-219相对于第一个文件，在第一次删除中删除了8行，因此此处不连贯，就比如说：假如第一次没有删除这8行，那就是从201开始，就连贯）。

  1. Visualize in IGV
    Use Filezilla to download:
    CS_Rep1_peaks_filter.narrowPeak
    CS_Rep2_peaks_filter.narrowPeak
    CS_Rep1_treat_pileup.bdg
    CS_Rep2_treat_pileup.bdg

在IGV中分别查看三个区域：
chrV:175,594-181,774
在这里插入图片描述
chrXII:666,176-671,951
在这里插入图片描述
chrXIV:449,980-454,181
在这里插入图片描述
注意：此处要右键选择“Autoscale”使其显示合理的尺寸。

(b) Call ChIP-seq peaks with GEM

java -Xmx8G -jar gem.jar --t 4 --q 0.05 --d Read_Distribution_default.txt --g \
sacCer3.chrom.sizes --genome . --exptCS ChIP-seq_SRR5179211_Rep1.bam --exptCS \
ChIP-seq_SRR5179212_Rep2.bam --f BAM --out PeakOutput --k_min 6 --k_max 13
  • -Xmx8G ：表示Java堆内存（Xmx）的最大容量为8G；
  • -jar：用于运行JAR（Java Archive）文件；
  • --t:Number of threads to run GEM in paralell.(default: physical CPU number)
  • --q:significance level for q-value, specified as -log10(q-value). For example, to enforce a q-value threshold of 0.001, set this value to 3.
  • --d：The path to the read distribution model file
  • --g：The path to a genome information file
  • --genome：The path to the genome sequence directory, which contains fasta files by chromosomes
  • --exptCS：The path to the aligned reads file for experiment (IP). X is condition name.
  • --f：Read file format: BED/SAM/BOWTIE/ELAND/NOVO. The SAM option allows SAM or BAM file input. (default = BED)
  • --out：Output folder and file name prefix
  • --k：the width of the k-mers
  • --k_min/max：minimum and maximum value of k

关于GEM的更多信息，可查看：GEM
在这里插入图片描述

Open ‘PeakOutput.results.htm’ in a browser of your choice. The ChIP-seq dataset is from the
Reb1 yeast protein with sequence specificity for the ‘TTACCCK’ motif. Did GEM correctly identify
the motif at the peak?

从结果HTML文件中，可以看到筛选出了’TTACCCG’motif，GEM identify the motif at the peak correctly.

3.3 ATAC-seq peak calling
  1. call peaks using MACS2
macs2 callpeak -t ATAC-seq_Rep1_SRR7696734.bam -f BAM -g 1.2e7 -n ATAC_Rep1 -B -q \
0.05 --shift -75 --extsize 150 --nomodel --SPMR --keep-dup all --call-summits
macs2 callpeak -t ATAC-seq_Rep2_SRR7696734.bam -f BAM -g 1.2e7 -n ATAC_Rep2 -B -q \
0.05 --shift -75 --extsize 150 --nomodel --SPMR --keep-dup all --call-summits
  • --shift:绝对的偏移值，会先于–extsize前对read进行整体移动。。正数表示从5’往3’偏移延长到片段中心，如果是负数则是3’往5’偏移延长到片段中心。ATAC-seq数据中，此值设置为-1/2*extsize。
  • --extsize:当设置了nomodel时，MACS会用–extsize这个参数从5’->3’方向扩展reads修复fragments。比如说你的转录因子结合范围200bp，就设置这个参数是200。双端测序可以预测read具体长度，因此不适用此参数。
  • --SPMR:进影响bedGraph文件的输出，将bedGraph文件中的读数转化为每百万reads的比例。
  1. Filter peaks with blacklist file：
bedtools intersect -v -a ATAC_Rep1_peaks.narrowPeak -b sacCer3_blacklist.bed > \
ATAC_Rep1_peaks_filter.narrowPeak
bedtools intersect -v -a ATAC_Rep2_peaks.narrowPeak -b sacCer3_blacklist.bed > \
ATAC_Rep2_peaks_filter.narrowPeak
  1. Visualize in IGV ：
    Use Filezilla to download:
    ATAC_Rep1_peaks_filter.narrowPeak
    ATAC_Rep2_peaks_filter.narrowPeak
    ATAC_Rep1_treat_pileup.bdg
    ATAC_Rep2_treat_pileup.bdg

chrV:175,594-181,774
在这里插入图片描述
chrXII:666,176-671,951
在这里插入图片描述
chrXIV:449,980-454,181
在这里插入图片描述

4. 实验总结

4.1 实验结论

ChIp-seq/ATAC-seq是表观基因组学(epigenomics)常用的两种测序手段，MACS2/GEM是peak calling中常用的工具，实验中，要注意去除黑名单序列，避免产生错误。

4.2 实验收获
  1. 学会了如何使用MACS2/GEM进行peak calling，对常用的各种参数有所了解；
  2. 学会了用IGV查看生物信息学中常见的bam、bgd、narrowpeak文件；
  3. 学会了用bedtools去除blacklist regions；
  4. 对ChIP-seq和ATAC-seq有最基本的认识。
4.3 其它心得

生物信息学中常用的一些工具/软件参数较多，可以通过官方网站学习各个参数的含义，再多加练习，以至于熟能生巧。

晚渔声
关注
ChIP-seq 数据分析1 ChIP-Seq技术2 ChIP-Seq数据分析
wx1871428的博客
07-08 3344
1 ChIP-Seq技术 1.1 概念 1.2 ChIP-seq技术原理 2 ChIP-Seq数据分析 2.1 数据下载 2.2 质量控制（data_assess） 2.3 比对到参考基因组（mapping_analysis） 2.4 搜峰（Peak_calling） MACS2 2.4.1 MACS2 核心: callpeak 用法 2.4.2 callpeak 结果文件说明 2.4.3 bdg file → wig file 2.5 峰注释（Peak_anno） ChIPseek.
CHIP-seq流程学习笔记（13）-ATAC_seq 数据加工处理
垚垚爸的博客
05-10 1251
今天第一次尝试处理ATAC_seq数据，希望能尽快做完吧。 先放个找好的参考文章：ATAC-seq/ChIP-seq分析方法 1.建立相应目录 对新数据建立对应实验人员（zhaoyingying）、测序类型（ATAC_seq）和日期（2021_05_03）的目录。 # 建立后如下： (base) zexing@DNA:~/projects/zhaoyingying/ATAC_seq/2021_05_03$ # 新建对应的目录 mkdir raw_data clean_data bam bam_bw bam
LanceOtron: a deep learning peak caller for ATAC-seq, ChIP-seq, and DNase-seq
dawnyi_yang的博客
07-10 535
ATAC-seq、ChIP-seq 和 DNase-seq 是基因组识别中重要的DNA编码元素，这些元素在分析覆盖轨迹的模拟信号中表现为峰值。本文提出了一个基于深度学习的峰值调⽤框架 LanceOtron，使用深度学习图像识别的方法来识别峰值形状，来进行富集测量。 ...
ATAC-seq学习记录
GeekFocus
10-07 3434
致密核小体结构被破坏后，启动子、增强子、绝缘子、沉默子等顺式调控元件和反式作用因子可以接近的特性，叫染色质的可及性，也叫染色质开放性(chromatin accessibility ），这段区域叫开放染色质（open chromatin） 。 2013年由斯坦福大学William J. Greenleaf和Howard Y. Chang实验室开发的ATAC-seq（Assay for Transposase-Accessible Chromatin with high throughput sequenc.
ATAC-seq技术原理简介
where is redhat
04-10 7093
ATAC-seq ATAC-seq (Assay for Transposase-Accessible Chromatin using sequencing) is a technique used in molecular biology to assess genome-wide chromatin accessibility.[1] In 2013, the technique was first described as an alternative advanced method for MNas
ChIP-seq 数据分析
热门推荐
fengzhibifang的博客
09-11 3万+
1 ChIP-Seq技术 1.1 概念 1.2 ChIP-seq技术原理 2 ChIP-Seq数据分析 2.1 数据下载 2.2 质量控制（data_assess） 2.2.1 质量评估 2.2.2 清理 reads 2.3 比对到参考基因组（mapping_analysis） fastq→sam→bam（H3K9me3 在全基因组的分布位置） 2.4 搜峰（Peak_callin...
在Windows下完成ChIP-Seq分析
Nuvolar
10-27 2382
纯粹就是记录下如何在Window下完成ChIP-Seq,并没有详细讲解ChIP-Seq的原理和每一步背后的含义。想深入了解可以翻生信技能树，生信宝典，组学大讲堂等等公众号。本文里的命令都是用(powershell or cmd, 部分需要cygwin64) 和 R脚本。 本文参考文章 Genome‑wide discovery of OsHOX24‑binding sites and regulation of desiccation stress response in rice 工具准备 conda （
ATAC-seq【Harvard FAS Informatics】
GeekFocus
10-14 1768
ATAC-seq数据质量评估注意 ENCODE的ATACseq数据标准。 必须有两次或更多次生物学重复（稀有样本也必须做两次技术重复） 每次重复要有25million非冗余，非线粒体，能够回帖的fragment（单端25 million reads，双端50 million reads=25 million fragment） 回帖率>95%, >80%可接受。用IDR计算重复一致性，rescue和self consisty ratios 都>2 用以下指标控制PCR扩增对文库复杂性的影
ATAC-seq分析：Peak Calling（8）
冷冻工厂
01-24 1868
1. 寻找开发区域 ATACseq 的一个共同目标是识别转录因子结合和/或转录机制活跃的无核小体区域。该核小体游离信号对应于小于一个核小体的片段（如 Greenleaf 论文中定义 < 100bp）。 然而，为了识别开放的染色质，我们可以简单地使用在测序中正确配对的所有读数（< 2000bp）。 2. 无核小体区域 有许多方法可用于从 ATACseq 数据中调用无核小体区域，其中许多方法借鉴自 ChIPseq 分析。一种非常流行且标准的 ATACseq 峰值调用程序是 MAC2。 2.1. 单端
CHIP-seq流程学习笔记（4）-call peak 软件macs2
垚垚爸的博客
09-22 4568
参考文章： MACS2 Call Peak 参数详细学习 1. 软件说明： macs2 callpeak --help usage: macs2 callpeak -t TFILE # treat组 [-c [CFILE]] # control 或 mock(非特异性抗体，如IgG)组 # control:input DNA，没有经过免疫共沉淀处理；
16/10/2019 一步步学会分析ATAC-seq
uqyuan0515的博客
10-29 2797
https://zhuanlan.zhihu.com/p/20702684 关于illumina 测序的原理 so, 我的测序结果是paired-end sequencing 并且是在不同的lanes上测序的。原始名称为： larvae1_S1_L001_R1_001.fastq.gz, larvae1_S1_L001_R2_001.fastq.gz, larvae1_S1_L002...
Postgresql特殊关键字
czjnoe的博客
09-25 327
表字段名称如果出现这些特殊关键字，则需要通过单引号访问字段，如offset、group。PostgreSQL特殊关键字。
MySQL之复合查询与内外连接
zdlynj的博客
09-30 626
前面我们讲解的mysql表的查询都是对一张表进行查询，即数据的查询都是在某一时刻对一个表进行操作的。而在实际开发中，我们往往还需要对多个表同时进行查询。我们这里使用的测试表，为雇员信息表(来自Oracle 9i的经典测试表)：EMP员工表，DEPT部门表，SALGRADE工资等级表。
k8s搭建一主三从的mysql8集群---无坑
oopxiajun博客专栏
09-27 902
对于构建基于MySQL的大规模、高性能应用来讲，需要使用水平扩展(集群)的数据库架构方式。在MySQL内建的复制功能可以实现，通过为服务器配置一个或多个备库的方式来进行数据同步。同时复制功能不仅有利于构建高性能的应用，也是高可用性、可扩展性、容灾、备份以及数据仓库等工作的基础。复制的基本原理是让一台服务器的数据与其他服务器保持同步。一台主库的数据可以同步到多台备库上，备库本身也可以被配置成另外一台服务器的主库。主库和备库之间可以有多种不同的组合方式。
JPA+Thymeleaf增删改查
y050505的博客
09-25 499
在使用JPA（Java Persistence API）和Thymeleaf作为模板引擎进行Web开发时，实现增删改查（CRUD）操作是一个常见的需求。下面，我将简要介绍如何使用Spring Boot框架结合JPA和Thymeleaf来实现这一功能。
[EBPF] 实时捕获DM数据库是否存在SQL阻塞
飞奔的屎壳郎
09-27 892
名称版本linux版本（centos8.5）DM88.1.4.27python利用ebpf技术在不登录数据库的情况下，通过对数据库的函数添加钩子，实时的对数据库是否存在阻塞判断现有的监控逻辑都是新建数据库用户执行SQL语句来判断阻塞情况，告警的实时性取决于监控的周期频率，而使用ebpf技术能够解决这个痛点ebpf有多种实现方式，这里便于理解测试的话 使用python语言内置c语言的形式进行了说明。涉及到生产环境应该是使用go语言或者c语言去做能够有效避免源码泄露等问题。
hive/impala/mysql几种数据库的sql常用写法和函数说明
biqidaoer的专栏
09-27 273
做大数据开发的时候，会在几种库中来回跳，同一个需求，不同库函数和写法会有出入，在此做汇总沉淀。
Redis事务
最新发布
zhzjn的博客
09-30 248
Redis事务可以通过MULTI、EXEC、WATCH和DISCARD指令来实现。1. MULTI指令：表示事务的开始，之后的所有指令都会被逐一记录，但不会立即执行。2. EXEC指令：表示事务的执行，会按照记录的指令逐一执行。3. WATCH指令：可以监视一个或多个键的变化，如果在事务执行期间被其他客户端修改，事务会被回滚。4. DISCARD指令：表示事务的撤销，可以放弃之前记录的指令。
检查并解释代码：#!/bin/bash #bwa ID=$(cat ID_list.txt) sample=$(cat sample_list.txt) for /date/lizichen/software/bwa-0.7.17/bwa mem -t 4 -R "@RG\tID:${ID}\tPL:ILLUMINA\tSM:${sample}" /date/lizichen/RAD-Seq/bwaout.file/GDDH13_1-1_formatted.fasta /date/lizichen/RAD-Seq/${sample}-I.1.fq.gz /date/lizichen/RAD-Seq/${sample}-I.2.fq.gz > /date/lizichen/RAD-Seq/bwaout.file/sam/${sample}.sam in $ID $sample do echo $ /date/lizichen/software/bwa-0.7.17/bwa mem -t 4 -R "@RG\tID:${ID}\tPL:ILLUMINA\tSM:${sample}" /date/lizichen/RAD-Seq/bwaout.file/GDDH13_1-1_formatted.fasta /date/lizichen/RAD-Seq/${sample}-I.1.fq.gz /date/lizichen/RAD-Seq/${sample}-I.2.fq.gz > /date/lizichen/RAD-Seq/bwaout.file/sam/${sample}.sam done
05-28
这是一个 Bash 脚本，用于运行 BWA 软件进行测序数据比对。下面是代码的解释和修正： ``` #!/bin/bash # 读取 ID 和样品名列表 ID=$(cat ID_list.txt) sample=$(cat sample_list.txt) # 循环运行 BWA 进行比对 for ID_sample in "${ID[@]} ${sample[@]}"; do # 运行 BWA，将比对结果存储为 SAM 文件 /date/lizichen/software/bwa-0.7.17/bwa mem -t 4 -R "@RG\tID:${ID_sample}\tPL:ILLUMINA\tSM:${sample}" /date/lizichen/RAD-Seq/bwaout.file/GDDH13_1-1_formatted.fasta /date/lizichen/RAD-Seq/${sample}-I.1.fq.gz /date/lizichen/RAD-Seq/${sample}-I.2.fq.gz > /date/lizichen/RAD-Seq/bwaout.file/sam/${sample}.sam # 输出 BWA 命令 echo "/date/lizichen/software/bwa-0.7.17/bwa mem -t 4 -R \"@RG\tID:${ID_sample}\tPL:ILLUMINA\tSM:${sample}\" /date/lizichen/RAD-Seq/bwaout.file/GDDH13_1-1_formatted.fasta /date/lizichen/RAD-Seq/${sample}-I.1.fq.gz /date/lizichen/RAD-Seq/${sample}-I.2.fq.gz > /date/lizichen/RAD-Seq/bwaout.file/sam/${sample}.sam" done ``` 修正后的代码中，`ID` 和 `sample` 分别改为数组类型，循环变量 `ID_sample` 用于遍历两个数组。在循环中，运行 BWA 进行比对，将比对结果存储为 SAM 文件。同时，输出 BWA 命令，方便调试和排错。
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值