B_两个NCBI的chip-seq一起下载为例子

如何判断是生物学重复或者技术重复？

·查看 SRA 或 GEO 页面：

o同一个 GSMxxxxxx 下列出的多个 SRR，通常属于同一个样本的不同测序 run（技术重复）。

o如果它们是生物学重复，一般会有独立的 GSM 编号（例如 GSM1_rep1、GSM1_rep2）。

·查看文章方法部分：会说明每个 ChIP 样本是否有多个测序 lane。

数据分析建议：做 ChIP-Seq 分析时，通常会把同一个 GSM 下的两个 SRR 数据合并，作为该对照组的完整测序数据进行后续分析（如质控、比对、peak calling）。

# 合并两个SRR的fastq文件（假设已经下载好） cat SRR36732730.fastq SRR36732731.fastq > DMSO_FOXA1_ChIP_combined.fastq

一、数据库准备

C．实际操作流程（以两个NCBI SRA 为例）：

#打开终端，用 echo 命令创建一个 SRR_Acc_List.txt 文件

echo -e "SRR36732730\nSRR36732731" > SRR_Acc_List.txt

#用 cat 命令检查一下文件内容是否正确

cat SRR_Acc_List.txt

#使用 prefetch 命令配合 --option-file 参数，可以一次性读取列表中的所有 SRR 号并开始下载

prefetch --option-file SRR_Acc_List.txt

#批量转换为 FASTQ（单端专用命令）, -e是线程数（建议不要超过逻辑线程）

for srr in $(cat SRR_Acc_List.txt); do

fasterq-dump$srr -e 16

done

#双端测序转换为FASTQ文件（--split-files：此参数会将双端测序数据自动分割为 _1.fastq 和 _2.fastq 两个文件）

for srr in $(cat SRR_Acc_List.txt); do

fasterq-dump$srr --split-files -e 16

done

二、质控

mkdir qcOutdir

fastqc SRR36732730.fastq SRR36732731.fastq -t 16 -o qcOutdir

cd qcOutdir

firefox SRR36732731_fastqc.html

#过滤

过滤主要有三步

#A. 去除adaptor序列

for sample in SRR36732730 SRR36732731; do

cutadapt-a AGATCGGAAGAGCACACGTCTGAACTCCAGTC \

${sample}.fastq\

-m30 -j 18 \

--output=${sample}.cutAd1.fq

done

#B. 去除N序列 (removeN是个脚本，读取一个FASTQ文件，计算每条read中 N或n（未确定碱基）占该read总长度的比例。

如果N比例 ≤ 5%（默认），这条read写入 *. unknowNul1.fq（好reads）

如果 N比例 > 5%，则写入 *.N.fq（N含量高的reads）)

输出的文件SRR36732730.unknowNul1.fq是合格的reads（N比例＜5%）

输出的文件SRR36732730. cutAd1.N.fq 是不合格reads （N比例＞5%）

for sample in SRR36732730 SRR36732731; do

./removeN${sample}.cutAd1.fq ${sample}.unknowNul1.fq

done

#C去除低质量reads

# basename：这是一个 Linux 命令，功能是从一个完整的文件路径中提取出文件名，并可以去掉指定的后缀。

# f：循环变量，依次代表列表中的每个原始文件名。第一次循环：f = SRR36732730.unknowNul1.fq

for f in SRR36732730.unknowNul1.fq SRR36732731.unknowNul1.fq; do

base=$(basename$f .unknowNul1.fq)

fastp-i $f -o ${base}.clean1.fq \

--qualified_quality_phred20 --unqualified_percent_limit 30 \

--thread16

done

#查看质控结果

for f in SRR36732730.clean1.fq SRR36732731.clean1.fq; do

fastqc$f -t 16 -o qcOutdir

done

三、对比

# 即mapping需要的软件：software: bowtie2 (bowtie2) + samtools (samtools)

conda install bowtie2 samtools

#这个命令试图用 Bowtie2 将过滤后的测序 reads 比对到人类参考基因组 hg19 上，生成一个 SAM 格式的比对文件。-p 16：使用 16个 CPU 线程。-U（单端）或 -1/-2（双端）明确指定。

for f in SRR36732730.clean1.fq SRR36732731.clean1.fq; do

base=$(basename$f .clean1.fq)# 提取 SRR36732730 或 SRR36732731

bowtie2 -p 16 -x /media/user/E18AF54B502996622/A_linux/D_chip_seq/A_hg19/hg19 -U $f -S ${base}.sam

done

#只留下唯一比对，去除多余比对

#先查看sam文件

less SRR36732730.sam

#去多余比对

·-b：输出 BAM 格式（二进制压缩格式），而不是默认的 SAM 文本格式。不加 -b 时，samtools view 输出 SAM 文本；加上 -b 则输出 BAM。

·-S：指明输入是 SAM 格式（而不是 BAM）。@后面是线程数

for sample in SRR36732730 SRR36732731; do

samtoolsview -h -@ 16 ${sample}.sam | grep -v 'XS:i:' | samtools view -bS -@ 16 - > ${sample}.unique.bam

done

#打开这个bam文件，查看头文件

samtools view -h SRR36732730.unique.bam | grep -c 'XS:i:'

#然后查看一下还有没有XS，这个代码出结果，最后面应该是0

四、call peak

#需要使用的软件macs3

conda install macs3

#call peak (peak finding) : -t 实验组；-c对照组；-f输入文件类型；-g 染色体选择的物种或者染色体长度（非模式物种）；-n 输出一大堆结果的头文件；-p p值的过滤条件；-B 输出bed文件

#用下面这个

macs3 callpeak -t Treatment_combined.unique.bam -c DMSO_combined.unique.bam -f BAM -g hs -n sample -p 1e-5 -B

#查看bed文件，按q键退出

less -SN sample_summits.bed

#主要是要查看xls文件

le sample_peaks.xls

#如果peaks还是非常多，那么后面也可以根据FDR值进行过滤

#查看summits.bed文件，会告诉我们每个peak的最高峰

# ，是实验组在某个区域比较低，对照组在这个区域反而是比较高的

# sample_peaks.narrowPeak 是 MACS2 / MACS3 输出的核心结果文件

#下面是peak的画图

#根据MACS3自己生成的R文件先画图，完事后会生成，sample_model.pdf这个文件

Rscript sample_model.r

#然后查看这个生成的pdf图

evince sample_model.pdf #双峰交汇处是正确的结合位点

#查看行号，为了是hist.R文件查看的时候去除多余的行，然后可以改峰的长度，更贴合xls表的数据

less -N sample_peaks.xls

#主要是需要根据sample_peaks.xls修改这个R脚本的行号。生成的是peakHist.pdf文件，看peak的长度

Rscript hist.R

evince peakHist.pdf

#进行peakview，peak的可视化

#得到的两个bdg文件，分别是treat和control。

#要修改peakView.R的具体参数

Rscript peakView.R

#Homer软件进行peak注释

conda install homer

#查看需要下载的注释基因组类型，如果没有，也可以自己设定

perl /home/user/miniconda3/envs/chip-seq_py3.10/share/homer/configureHomer.pl -list

#同时还要安装基因组注释软件，bowtie2用哪个，现在的注释软件就用哪个#同时还有四个依赖软件

conda install ghostscript weblogo ucsc-blat samtools

#bed文件的输出的sample_peaks.narrowPeak转为homer的输入

awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' sample_peaks.narrowPeak > sample_homer.bed

less -N sample_homer.bed

#安装基因组

perl /home/user/miniconda3/envs/chip-seq_py3.10/share/homer/.//configureHomer.pl -install hg19

#找motif。按照8,10,12的长度去寻找。motifDir是新建的文件夹存放文件的。

findMotifsGenome.pl sample_homer.bed hg19 motifDir -len 8,10,12

#已知基序富集：knownResults.html 相关文件，帮你回答「我的 ChIP 峰里，哪些已知转录因子的结合基序显著富集？」

#从头基序预测：homerResults.html 相关文件，帮你回答「我的 ChIP 峰里，有没有新的、未被报道的转录因子结合基序？」

cd motifDir

cd ..

#对peak进行注释。1是注释结果，2是错误结果。目前没有好的方法多文件并行

annotatePeaks.pl sample_homer.bed hg19 1>peakAnn.xls 2>annLog.txt

#也可以对peak进行kegg和go注释(在GODir文件中)

annotatePeaks.pl sample_homer.bed hg19 -go GODir 1>peakAnn.xls 2>annLog.txt

less -N peakAnn.xls

#统计peakAnn.xls中，第8列的类型

awk -F "\t" 'NR>1{ split($8, arr, "("); sub(" ", "",arr[1]); print arr[1]}' peakAnn.xls | sort | uniq -c | awk '{print $2"\t"$1}' > peakRegionStat.xls

less -N peakRegionStat.xls

Rscript pie.R

evince peakRegionStat.pdf#查看peak存在哪些区域