专栏名称: 生信菜鸟团
生信菜鸟团荣誉归来,让所有想分析生物信息学数据的小伙伴找到归属,你值得拥有!
目录
相关文章推荐
生物制品圈  ·  A股IPO排队医药企业名单:仅17家等待上市 ... ·  昨天  
BioArt  ·  Sci ... ·  17 小时前  
生信人  ·  多组学联合如何关联临床,全国可约 ·  2 天前  
BioArt  ·  ​Nat ... ·  昨天  
BioArt  ·  Mol Cell | ... ·  2 天前  
51好读  ›  专栏  ›  生信菜鸟团

ChIP-Seq 分析流程-上游

生信菜鸟团  · 公众号  · 生物  · 2025-02-10 21:17

主要观点总结

本文介绍了使用HBC教学团队成员开发的文件资料进行的ChIP-seq分析流程,包括数据下载、质量控制、比对到参考基因组、peak calling等步骤。

关键观点总结

关键观点1: 数据下载和准备

文章介绍了从公开数据库下载原始测序数据的方法,包括原始FASTQ文件的下载和参考基因组的下载。

关键观点2: 质量控制

使用seqkit工具对原始数据进行统计,以检查数据质量。

关键观点3: 比对到参考基因组

使用Bowtie2工具进行基因序列比对,包括构建索引、比对、SAM转BAM、BAM排序和过滤唯一映射的读取等步骤。

关键观点4: Peak calling

使用macs2进行peak calling,找出转录因子结合的区域。

关键观点5: 后续分析

文章提到下游分析和每一步的详细讨论后续会继续更新。


正文

写在前面的话:

参考使用的文件资料是由哈佛陈生物信息学核心 (HBC) 教学团队成员开发的。另外也看了多个公众号文章和书籍。

Website: https://hbctraining.github.io/Intro-to-ChIPseq/

Github: https://github.com/hbctraining/Intro-to-ChIPseq

碎碎念:因为会有代码存在,代码框中既有命令,又有注释和输出信息,为了避免误解,这里统一一下(仅限代码框内)。

“#” 字符后代表注释信息

“$” 字符后代表输入的命令代码, 可主要关注这个

其余信息均为命令输出内容

全流程概述

分析工具

目前来说,分析工具肯定不止这些,每一步都有可替换的工具,甚至有些工具已经有所更新。但我感觉都差不多,感觉不行咱就换!

数据格式

实际上,生信分析只是分析数据。软件工具只是将数据从一种格式转换为另一种格式。新生成的数据通常会有不同的 优化 。它可能更加丰富,因为它可能结合了来自其他数据的信息。但它仍然是数据——生物学洞察力来自于对这些数据的理解和研究。

分析过程

创建目录

让我们直接开始吧!

这里使用的是学习资料给到的目录结构,但这个结构好坏因人而异,这里为了方便,我也就全盘接收了。

  1. 创建项目名称
$ mkdir chipseq
cd chipseq/
  1. 创建分析目录
$ mkdir raw_data reference_data scripts logs meta
$ mkdir -p results/fastqc results/bowtie2
  1. 查看目录结构
# 使用tree 查看目录结构
$ tree
.
|-- logs
|-- meta
|-- raw_data
|-- reference_data
|-- results
|   |-- bowtie2
|   `-- fastqc
`-- scripts
  • logs :跟踪运行的命令和使用的具体参数,同时还记录运行命令时生成的任何标准输出。
  • meta :用于描述您正在使用的样本的任何信息,我们将其称为元数据
  • raw_data :用于计算分析之前获得的任何 未修改的(原始)数据,例如来自测序中心的 FASTQ 文件。我们强烈建议在分析过程中不要修改此目录。
  • reference_data :用于分析的参考基因组相关的已知信息,例cd如基因组序列(FASTA)、与基因组相关的基因注释文件(GTF)。
  • results :用于您在工作流程中实施的不同工具的输出。在此文件夹中创建特定于工作流程的每个工具/步骤的子文件夹。
  • scripts :用于您编写并用于运行分析/工作流程的脚本。

下载原始数据

原始数据依旧是 FASTQ 格式。

此资料使用的数据包含六个样本。其中两个样本作为 control.

这里简单解释一下样本情况:这个测序的目的是找到 Nanog 和 Pou5f1 两个转录因子的结合位点。在实验部分,首先将细胞分为三份(不一定等分),分别不做抗体富集(input control)、Nanog 抗体富集和 Pou5f1 抗体富集,生成三组数据,再进行重复,生成六个测序样本。

原始数据链接

原教程的原始数据比较难找,就不详细说明过程了,这里直接给出下载 URL 好了。

  • Nanog
# rep1
https://www.encodeproject.org/files/ENCFF000OQA/@@download/ENCFF000OQA.fastq.gz
# rep2
https://www.encodeproject.org/files/ENCFF000OQM/@@download/ENCFF000OQM.fastq.gz
  • Pou5f1
# rep1
https://www.encodeproject.org/files/ENCFF000ORV/@@download/ENCFF000ORV.fastq.gz
# rep2
https://www.encodeproject.org/files/ENCFF000ORW/@@download/ENCFF000ORW.fastq.gz
  • Input
# rep1
https://www.encodeproject.org/files/ENCFF000OSS/@@download/ENCFF000OSS.fastq.gz
# rep2
https://www.encodeproject.org/files/ENCFF000OST/@@download/ENCFF000OST.fastq.gz

批量下载

现在创建一个批量下载脚本

$ cat down_fq.sh 
#!/bin/bash
URL=~/project/chipseq/raw_data/url.txt
dir=~/project/chipseq/raw_data

cat ${URL} | parallel --dry-run aria2c -x 5 -d ${dir} {}
cat ${URL} | parallel aria2c -x 5 -d ${dir} {}

下载后解压

$ gunzip  *.gz

查看原始数据情况,我这里推荐使用 seqkit 工具去统计一下我们的原始数据情况,可以比对下看是否有缺失。

$ ls *fastq
ENCFF000OQA.fastq  ENCFF000OQM.fastq  ENCFF000ORV.fastq  ENCFF000ORW.fastq  ENCFF000OSS.fastq  ENCFF000OST.fastq
$ seqkit stats *fastq
processed files:  6 / 6 [======================================] ETA: 0s. done
file               format  type    num_seqs        sum_len  min_len  avg_len  max_len
ENCFF000OQA.fastq  FASTQ   DNA   15,329,142    551,849,112       36       36       36
ENCFF000OQM.fastq  FASTQ   DNA   63,707,371  2,293,465,356       36       36       36
ENCFF000ORV.fastq  FASTQ   DNA   18,153,153    653,513,508       36       36       36
ENCFF000ORW.fastq  FASTQ   DNA   35,253,613  1,269,130,068       36       36       36
ENCFF000OSS.fastq  FASTQ   DNA   55,951,481  2,014,253,316       36       36       36
ENCFF000OST.fastq  FASTQ   DNA   27,548,706  1,377,435,300       50       50       50

为了便于查看,我们可以手动改个名字,我根据需要将六个样本改成了这样。

$ ls *fastq 
H1hesc_Input_Rep1.fastq  H1hesc_Input_Rep2.fastq  H1hesc_Nanog_Rep1.fastq  H1hesc_Nanog_Rep2.fastq  H1hesc_Pou5f1_Rep1.fastq  H1hesc_Pou5f1_Rep2.fastq

下载参考基因组

cd reference_data/

# 下载fa文件
$ aria2c -x 5 -d ./ https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.40_GRCh38.p14/GCF_000001405.40_GRCh38.p14_genomic.fna.gz

# 下载gtf文件
$ aria2c -x 5 -d ./ 
https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.40_GRCh38.p14/GCF_000001405.40_GRCh38.p14_genomic.gtf.gz

解压

$ gunzip  *.gz

质量评估

cd raw_data
$ fastqc *.fastq

所有数据下载完后就可以正式进行数据分析了。简单进行数据质量评估。

质评结果就不详细说明了,因为现在测序技术越来越成熟,很多数据本来质量就很好,就不需要过滤了,如果查看评估结果后觉得需要过滤,就使用一些常见的过滤软件默认参数过滤一下就好了,基本上问题不大。

Bowtie 比对

Bowtie2是一种快速准确的对齐工具,它使用基于 Burrows-Wheeler 变换方法的 FM 索引对基因组进行索引,以保持对齐过程的内存要求较低。Bowtie2支持间隙、局部和双端对齐模式,最适合至少 50 bp 的读取(较短的读取长度应使用 Bowtie1)。默认情况下, Bowtie2 将执行全局端到端读取对齐,这最适合质量修剪的读取。但是,它还有一种局部对齐模式,它将执行软剪辑以从未修剪的读取中去除质量较差的碱基或适配器。我们将使用此选项,因为我们没有修剪我们的 reads。

构建索引

如果没有构建索引,第一步需要构建索引。

cd reference_data
$ bowtie2-build GCF_000001405.40_GRCh38.p14_genomic.fna GCF_000001405.40_GRCh38.p14_genomic

比对

我们先对单个原始 fastq 文件进行处理

$ mkdir results/bowtie2

$ bowtie2 -p 4 -q --local \
-x reference_data/GCF_000001405.40_GRCh38.p14_genomic \
-U raw_data/H1hesc_Input_Rep1.fastq \
-S results/bowtie2/H1hesc_Input_Rep1_aln_unsorted.sam

过滤reads

ChIP-seq数据中包含多重映射的读取(reads映射到参考基因组的多个位置)。允许多重映射的读取可以增加可用读取的数量和峰值检测的灵敏度;然而,假阳性的数量也可能增加。因此,我们需要过滤我们的对齐文件,只包含唯一映射的读取,以提高位点发现的置信度和改善可重复性。由于Bowtie2中没有参数可以只保留唯一映射的读取,我们需要执行以下步骤来生成只包含唯一映射读取的对齐文件:

  1. 将对齐文件格式从SAM转换为BAM
  2. 按读取坐标位置对BAM文件进行排序






请到「今天看啥」查看全文