Chip-seq上游分析实战：从数据下载到质控全流程解析

1. 环境准备与软件安装别在第一步就卡住大家好我是老张在生信分析这个坑里摸爬滚打十来年了今天咱们来聊聊Chip-seq上游分析这个事儿。很多刚入门的朋友尤其是学生物的同学一看到命令行就头疼觉得这玩意儿比实验还难。其实不然只要你把环境搭好软件装对后面的流程就像搭积木按部就班来就行。咱们今天的目标就是让你能亲手把GSE274995这个数据集从头到尾跑一遍拿到像样的质控报告和比对结果。首先你得有个地方干活我们管这个叫“工作环境”。现在最省心的方式就是用Conda。你可以把它理解成一个“软件管家”帮你把不同分析流程需要的软件和它们各自的依赖包分别装在不同的“小房间”环境里互相不打架。比如你跑Chip-seq需要一个Python 3.9的环境跑别的项目可能需要Python 3.11用Conda就能轻松搞定。原始文章里提到了一个关键问题芯片架构。如果你用的是苹果的M1或M2芯片电脑ARM架构直接安装某些生信软件可能会遇到兼容性问题。我当初用M1 MacBook Pro的时候就踩过这个坑像bowtie2的一些版本就是装不上报错信息看得人云里雾里。解决办法其实很简单就是在创建Conda环境时强制指定使用传统的x86_64架构的软件包。具体命令如下# 针对苹果M1/M2芯片的用户 CONDA_SUBDIRosx-64 conda create -n chipseq_x86_64 python3.9 conda activate chipseq_x86_64这条命令里的CONDA_SUBDIRosx-64就是关键它告诉Conda“别给我ARM架构的包我要老式的Intel芯片能用的包”。这样创建出来的环境绝大多数生信软件都能顺利安装。当然如果你用的是Windows通过WSL2、Linux服务器或者Intel芯片的Mac那就省事了直接创建环境就行# 适用于大多数Linux服务器及Intel Mac用户 conda create -n chipseq python3.9 conda activate chipseq环境建好了接下来就是装软件。原始文章里列了一串sra-tools, trim-galore, samtools, deeptools, homer, meme, macs2, bowtie, bowtie2。如果你直接用conda install来装有时候会因为解析软件依赖关系特别慢等得人心急。我这里分享一个提速小技巧先安装mamba。mamba是conda的C重写版速度更快用它来安装软件体验好很多。# 先安装mamba conda install -c conda-forge mamba # 然后用mamba一次性安装所有需要的软件 mamba install -y sra-tools trim-galore samtools deeptools homer meme macs2 bowtie bowtie2这里解释一下几个核心软件是干嘛的让你心里有数sra-tools 里面包含prefetch和fastq-dump用来从NCBI的SRA数据库下载数据并转换成fastq格式。这是数据入口。trim-galore 它其实是Cutadapt和FastQC的包装用来做数据质控和修剪接头、低质量碱基一步到位。samtools 处理SAM/BAM格式比对文件的“瑞士军刀”排序、索引、格式转换都靠它。bowtie2 将测序reads比对到参考基因组上的工具速度快内存占用相对友好是Chip-seq比对的常用选择。macs2 下游找峰Peak Calling的核心工具不过咱们上游分析暂时用不到先装上备用。deeptools 用于生成各种可视化图比如覆盖度曲线、热图等质控和结果展示利器。homer meme 主要用于下游的motif分析和上游关系不大但既然装了环境一并装上也无妨。软件安装这一步最怕的就是网络问题。如果你觉得下载速度慢可以给Conda配置国内镜像源比如清华源或者中科大源速度会快很多。具体配置方法网上教程很多这里就不展开了。总之这一步顺利完成你的“厨房”就算备好料了接下来可以开始“烹饪”数据了。2. 数据下载与格式转换搞定原始数据环境准备好了我们得把“食材”——测序数据弄到手。这次我们实战用的数据集是GSE274995研究的是头颈部鳞癌细胞系Cal27。在NCBI的GEO数据库找到这个数据集后里面通常会有一个“SRA Run Selector”的链接点进去就能看到具体的样本列表比如SRR1234567这种编号。我们把所有需要的样本编号保存到一个文本文件里比如叫SRR_Acc_List.txt每行一个编号。下载原始数据我们用prefetch命令。原始文章里给出了一条复合命令用了nohup和放到后台运行这对于在服务器上操作、担心网络中断的同学非常实用。我们来拆解一下这条命令nohup cat SRR_Acc_List.txt | while read id; do prefetch $id; done download.log 21 这条命令看起来有点复杂我们一点点看nohup 保证命令在后台持续运行即使你关闭了终端窗口下载也不会停止。cat SRR_Acc_List.txt | while read id; do ... donecat读取列表文件通过管道|传给while循环。read id把每一行的内容读出来赋值给变量id。这个id可以是任何名字比如line、srr都行。prefetch $id 在循环体内对每一个id即SRR编号执行prefetch命令进行下载。 download.log 21 这是重定向组合拳。将命令的标准输出正常信息重定向到download.log文件。21的意思是“将标准错误2也重定向到标准输出1当前所在的地方”也就是也写入download.log。这样所有运行信息包括报错都记录在一个文件里方便查看。最后的 让整个命令在后台执行不占用当前终端。执行后你可以用tail -f download.log来实时查看下载日志。下载下来的文件是.sra格式这是NCBI的存储格式我们分析需要的是.fastq格式。转换就用fastq-dump。这里有个细节要注意先搞清楚数据是单端Single-end还是双端Paired-end测序。在SRA样本页面查看“Layout”属性。我们用的GSE274995数据是“SINGLE”即单端测序。如果是双端fastq-dump需要加上--split-files参数来生成两个文件。我习惯写一个简单的Shell脚本来批量转换这样更清晰也避免每次手动输入。创建一个文件比如叫sra_to_fastq.sh#!/bin/bash # 定义路径请根据你自己的实际情况修改 work_dir/home/your_project/chipseq sra_dir${work_dir}/sra # 存放.sra文件的目录 fastq_dir${work_dir}/fastq # 输出.fastq.gz的目录 list_file${work_dir}/SRR_Acc_List.txt # 创建输出目录 mkdir -p ${fastq_dir} # 循环处理 cat ${list_file} | while read srr_id do sra_file${sra_dir}/${srr_id}.sra if [ -f ${sra_file} ]; then echo 正在转换 ${srr_id} ... # 关键命令--gzip直接输出压缩格式-O指定输出目录 fastq-dump --gzip -O ${fastq_dir} ${sra_file} echo ${srr_id} 转换完成。 else echo 警告文件 ${sra_file} 不存在跳过。 fi done echo 所有转换任务已完成。给脚本加上执行权限 (chmod x sra_to_fastq.sh)然后运行 (./sra_to_fastq.sh) 即可。用--gzip参数直接输出压缩的.fastq.gz文件能节省大量的磁盘空间。转换完成后建议用ls -lh ${fastq_dir}看一眼文件大小再用zcat ${fastq_dir}/SRRxxx.fastq.gz | head -n 4快速瞄一眼前几条read确认文件没问题。通常一条read占4行序列ID、碱基序列、分隔符、质量值。3. 原始数据质控用FastQC和MultiQC给数据“体检”数据拿到手了千万别急着往下分析。这就好比做实验你拿到一批新试剂总得先测个OD值、跑个胶看看纯度吧测序数据也是一个道理原始数据里可能藏着小问题比如测序质量随循环数下降、接头污染、碱基组成异常等。质控QC这一步的目的就是发现这些问题并决定是否需要、以及如何进行修剪过滤。我们用的主力工具是FastQC。它不修改你的数据只是生成一份详细的“体检报告”。原始文章里用了fastqc -t 6 -o ./ ./SRR*.fastq.gz这个命令。我来解释一下-t 6 使用6个CPU线程并行运行加快速度。这个数字根据你电脑或服务器的核心数来调整。-o ./ 指定报告的输出目录为当前目录。./SRR*.fastq.gz 通配符*表示处理当前目录下所有以SRR开头、以.fastq.gz结尾的文件。运行后每个样本都会生成一个.zip压缩包和一个.html网页报告。打开.html文件你会看到十几个模块我挑几个最重要的说说怎么看Basic Statistics 看一眼总序列数、序列长度、GC含量等基本信息确认数据量符合预期。Per base sequence quality这是重中之重它展示每个测序循环位置的碱基质量。横坐标是位置纵坐标是质量分数Q值。Q20错误率1%、Q30错误率0.1%是常用阈值。理想情况是整条线都在绿色区域高质量区或只在末尾稍有下降。如果开头或中间就掉到黄色甚至红色说明测序起始质量不好或有其他问题。Per sequence quality scores 看所有序列的平均质量分布。通常应该是一个单一的高峰。Per base sequence content 检查每个位置A/T/C/G四种碱基的组成比例。在随机文库中每个位置的比例应该接近~25%并且四条线平行。如果开头几条线剧烈波动可能是接头污染或随机引物偏倚。Adapter Content 直接告诉你检测到了多少比例的接头序列。如果比例很高比如5%在后续的过滤步骤就必须把它切掉。当你样本不多时一个个点开.html报告看还行。但如果有几十上百个样本呢眼睛都得看花。这时候MultiQC就派上大用场了。它能自动搜集所有FastQC生成的.zip或.html文件整合成一份统一的、交互式的报告。命令简单到令人发指# 假设所有FastQC的.zip结果都在当前目录 multiqc . -o ./multiqc_report运行后在./multiqc_report目录下会生成multiqc_report.html。打开它你可以轻松地比较所有样本在各个QC指标上的表现一眼就能看出哪个样本质量差问题出在哪里。这才是高效的批量质控方式。原始文章也提到了这一点但我想强调先跑fastqc生成单个报告再用multiqc汇总这个顺序不能错。4. 数据过滤与修剪用Trim Galore做“数据清洗”看完FastQC报告如果发现了接头残留、或者3‘端质量普遍偏低等问题我们就需要对数据进行“清洗”也就是过滤和修剪。这一步的目标是去除低质量序列、接头序列以及过短的序列保留下高质量的数据用于后续比对。原始文章提到了fastp和trim-galore我这里强烈推荐新手使用Trim Galore。为什么呢因为它本质上是一个智能脚本自动调用了两个非常优秀的工具Cutadapt用于切除接头和FastQC用于质控。它简化了参数设置对常见的Illumina接头序列能自动检测而且能在修剪后自动再跑一次FastQC让你直观看到修剪前后的对比效果非常贴心。一个最基础的Trim Galore命令长这样trim_galore --quality 20 --length 20 --fastqc --gzip -o ./trimmed_reads SRR1234567.fastq.gz我们来拆解一下关键参数--quality 20 设定质量阈值。它会从3‘端开始扫描如果滑动窗口默认4个碱基的平均质量低于20就把后面部分全部切掉。Q20是常用标准要求高可以设为--quality 30。--length 20 设定保留序列的最短长度。修剪后如果序列长度小于20个碱基这条序列就会被直接丢弃。太短的序列比对时特异性差容易造成假阳性。--fastqc 修剪完成后自动对修剪后的文件运行FastQC生成新的质控报告。这是Trim Galore的一大亮点。--gzip 输出压缩的.fastq.gz文件。-o ./trimmed_reads 指定输出目录。最后跟上你的输入文件。对于双端测序数据需要同时指定两个文件并加上--paired参数这样它能保证一对reads被同时保留或丢弃避免后续比对出错。在实际项目中我们同样需要批量处理。可以写一个循环脚本#!/bin/bash input_dir./fastq output_dir./trimmed_reads mkdir -p ${output_dir} for fq in ${input_dir}/*.fastq.gz do base$(basename ${fq} .fastq.gz) echo 正在处理 ${base} ... trim_galore --quality 20 \ --length 20 \ --fastqc \ --gzip \ -o ${output_dir} \ ${fq} done echo 数据修剪全部完成。修剪后报告位于: ${output_dir}运行完检查输出目录。你会看到每个样本产生了几个新文件修剪后的.fastq.gz文件、一个修剪报告的.txt文件、以及一个新的FastQC报告。一定要打开那个.txt报告看看里面统计了原始读长、修剪掉的碱基数、丢弃的序列数等信息。然后再对比一下修剪前后的FastQC报告特别是“Per base sequence quality”和“Adapter Content”这两个图应该能看到明显的改善。数据清洗干净了我们才能放心地进行下一步——把reads“贴”到参考基因组上。5. 序列比对与SAM/BAM文件处理找到reads的“家”清洗干净的数据现在要找到它们在基因组上的位置这个过程就叫比对。Chip-seq分析里Bowtie2是最常用的比对工具之一它在速度和灵敏度之间取得了很好的平衡。比对的“地图”就是参考基因组我们需要先下载并构建索引。原始文章提到数据使用了hg19但现在强烈推荐使用更新的GRCh38或hg38因为其组装质量和注释都更完善。首先下载参考基因组序列FASTA文件和基因注释文件GTF文件。我推荐从ENSEMBL或UCSC官网下载。假设我们下载了Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz这个文件。接下来用bowtie2-build为这个基因组构建索引。这是一个非常耗时的步骤但只需要做一次。构建好的索引是一组以.bt2结尾的文件。# 解压基因组文件 gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz # 构建Bowtie2索引 (假设文件在当前目录) bowtie2-build Homo_sapiens.GRCh38.dna.primary_assembly.fa grch38_index # 这会生成 grch38_index.1.bt2, grch38_index.2.bt2 等文件索引准备好后就可以进行比对了。对于单端数据如我们的GSE274995命令如下bowtie2 -x ./grch38_index \ -U ./trimmed_reads/SRR1234567_trimmed.fq.gz \ -S ./alignment/SRR1234567.sam \ --threads 8 \ --very-sensitive-local 2 ./alignment/SRR1234567.bowtie2.log参数解释-x 指定索引的路径和前缀这里是./grch38_index。-U 指定单端测序文件Unpaired reads。-S 指定输出的SAM格式文件。--threads 8 使用8个线程加速。--very-sensitive-local 这是预设的比对模式之一在局部比对模式下提供更高的灵敏度适合Chip-seq允许部分比对捕捉结合位点信号。你也可以根据数据情况选择其他模式如--sensitive-local。2 ...log 将Bowtie2运行的统计信息标准错误重定向到一个日志文件方便查看总体比对率。比对完成后会生成SAM文件。这是一种文本格式可读但体积庞大。我们需要用samtools将它转换为二进制的BAM格式体积小并排序、建立索引以便后续工具快速访问。# 1. SAM转BAM并按坐标排序 samtools view - 4 -bS SRR1234567.sam | samtools sort - 4 -o SRR1234567.sorted.bam # 2. 为排序后的BAM文件建立索引生成 .bai 文件 samtools index SRR1234567.sorted.bam这里用了管道|将两个命令连起来避免生成巨大的中间文件。- 4指定使用4个线程。现在你得到了SRR1234567.sorted.bam和SRR1234567.sorted.bam.bai两个文件。你可以用samtools flagstat SRR1234567.sorted.bam快速查看比对的统计摘要比如总共多少条reads有多少比对上有多少是唯一比对等等。这个比对率是衡量实验和数据质量的一个重要指标通常希望至少在70%以上。6. 比对后质控与深度检查确保比对结果可靠你以为比对完生成BAM文件就结束了吗还早呢。比对的质量和深度直接决定了后续找峰Peak Calling的可靠性。这一步的质控主要关注几个方面比对率、插入片段长度分布、PCR重复率、以及测序深度。首先samtools flagstat给出的总体比对率是一个宏观指标。我们还需要更细致的检查。一个常见问题是PCR重复。由于PCR扩增来自同一个模板分子的多个拷贝会被测序产生完全相同的reads。这些重复会干扰信号高重复率可能意味着实验起始量不足或PCR循环数过多。我们可以用samtools markdup需要先按坐标排序或者picard MarkDuplicates来标记并移除重复。# 使用samtools标记重复 (需要已排序的BAM) samtools markdup -r -s SRR1234567.sorted.bam SRR1234567.dedup.bam # -r: 直接移除重复而不是仅标记。 -s: 打印统计信息到屏幕。其次对于ChIP-seq数据尤其是转录因子插入片段长度是一个关键参数。它代表了蛋白质结合DNA后被交联和片段化的大小。估计这个长度对于后续分析如MACS2建模很重要。我们可以用samtools配合awk简单估算或者用更专业的工具如picard CollectInsertSizeMetrics。一个简单的估算方法是提取正确配对的reads对于单端数据这一步需要一些假设双端数据更准确然后计算模板长度TLEN字段的分布。这里不展开复杂脚本但你需要知道这个概念。最后也是非常重要的一点测序深度。深度不够信号太弱找不到真正的峰深度过高又浪费资源还可能引入更多噪音。我们可以用samtools depth来计算基因组每个位点的覆盖深度然后看全局的平均深度。# 计算平均深度 samtools depth SRR1234567.dedup.bam | awk {sum$3} END {print 平均深度 , sum/NR}对于转录因子ChIP-seq有效深度去除重复后唯一比对的reads数通常在1千万到2千万条reads就足够了。对于组蛋白修饰可能需要更高的深度。你可以用samtools view -c来统计去重后的reads数。所有这些质控指标最好能整合进一份报告。除了之前提到的MultiQC它可以整合Bowtie2的日志、samtools flagstat的结果等还可以用qualimap或deepTools来生成更丰富的比对质量评估图。例如plotFingerprint来自deepTools可以用来评估样本的富集效果是否理想。完成这一步并确认数据质量合格后你的上游分析才算真正完成。干净、可靠的BAM文件是下游进行Peak Calling、可视化如用bamCoverage生成BigWig文件和生物学解读的坚实基础。整个流程从数据下载到质控完成虽然步骤不少但每一步都有明确的目的和工具。多跑几次把流程写成脚本固化下来以后遇到新的ChIP-seq数据你就能从容应对了。