Chip-seq数据分析
发布时间
阅读量：
阅读量
1.数据质量评估
1.fastqc质控
# 进入到到自己的文件夹，使用FastQC软件对fastq文件进行质量评估，结果输出到qc/文件夹下
vim qc.sh
i
fastqc -t 6 -o ./ *.fq.gz
# 按esc
:wq
nohup bash qc.sh >qc.log &
# 报告整合
multiqc *zip
cat qc.sh
cat qc.log
# 质控的结果给每个fq文件生成一个zip文件和一个html网页报告
# zip文件里是一些报告图片等，主要需要看的是html网页报告，可以将其下载到本地通过浏览器打开
#另外一个跑qc的代码：
#放前台跑
ls ../raw/*gz |xargs fastqc -t 10 -o ./
2.使用Trim Galore修剪低质量reads
# 使用Trim Galore修剪低质量reads
trim_galore --paired --fastqc --length 50 input_R1.fastq input_R2.fastq
--paired 表示处理双端测序数据。
--fastqc 表示在修剪后生成FastQC报告以评估数据质量。
#1：单个数据处理示例代码如下
trim_galore --fastqc --gzip --paired ../AJV32_H2AK119ub_1.fq.gz ../AJV32_H2AK119ub_2.fq.gz
trim_galore --fastqc --gzip --paired ../JV21_H2AK119ub_1.fq.gz ../JV21_H2AK119ub_2.fq.gz
##尝试批量处理数据：
vim trim_galore.sh
analysis_dir=/Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/clean/
bin_trim_galore="trim_galore"
cd ../raw/
#1：批量处理双端测序数据
for fq1 in *_1.fq.gz; do
# 构造R2文件名
fq2=${fq1/_1/_2}
# 检查R1和R2文件是否存在:
if [ -f "$fq1" ] && [ -f "$fq2" ]; then
# 使用nohup运行trim_galore，处理R1和R2文件
$bin_trim_galore --fastqc --gzip --paired -o $analysis_dir $fq1 $fq2
else
echo "警告：文件 $fq1 或 $fq2 不存在，跳过处理"
fi
done
nohup bash trim_galore.sh > trim_galore.sh.log &
#仿照RNA-seq的trim_galore过滤
vim trim_galore.sh
rawdata=/Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/raw/
cleandata=/Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/clean2/
ID=/Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/clean2/ID
cat $ID | while read id
do
trim_galore -q 20 --length 35 --max_n 3 --stringency 3 --fastqc --gzip --paired -o ${cleandata} ${rawdata}/${id}_1.fq.gz ${rawdata}/${id}_2.fq.gz
done
nohup bash trim_galore.sh >trim_galore.log &
过滤完可以再查看QC一下：
3.序列比对
1、比对
使用比对软件（如BWA、Bowtie2等）将读段比对到参考基因组上。这个过程中，读段会根据其序列相似性被映射到基因组的特定位置。
#使用Bowtie2进行序列比对
#构建参考基因组的索引或者已经有现成的索引。
#如果没有，需要先使用 Bowtie2 构建索引。
#1.构建索引
# 进入参考基因组目录
mkdir bowtie2Index
vim bowtie2Index.sh
fa=Mus_musculus.GRCm39.dna.primary_assembly.fa
fa_baseName=GRCm39.dna
bowtie2-build -f ${fa} bowtie2Index/${fa_baseName}
#运行
nohup bash bowtie2Index.sh > bowtie2Index.log &
## ----2.比对
# 到bowtie2文件夹下先生成一个变量,为样本ID
cd /Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/align/
ls /Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/clean/*_1.fq.gz |awk -F '/' '{print $NF}' | cut -d '_' -f 1 > ID
#bowtie2批量处理多个样本：
cd /Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/clean2/
vim bowtie2.sh
# 设置输入数据目录
inputdir="/Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/clean2/"
# 设置ID文件路径
ID="/Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/clean2/ID"
# 设置输出结果目录
outdir="/Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/align/"
# 设置bowtie2索引目录
index=/Data/wu_lab/zyy/database/GRCm39/bowtie2Index/GRCm39.dna
# 确保输出目录存在
mkdir -p "$outdir"
# 读取ID文件并进行批量比对
cat $ID | while read id; do
# 执行bowtie2比对并将结果排序为BAM格式
bowtie2 -p 5 -x ${index} \
-1 ${inputdir}/${id}_1_val_1.fq.gz \
-2 ${inputdir}/${id}_2_val_2.fq.gz \
--phred33 --end-to-end --very-sensitive --no-mixed --no-discordant --no-unal | samtools sort -@ 16 -O bam -o ${outdir}/${id}.bam
echo "Alignment and sorting for $id completed."
done
nohup bash bowtie2.sh > bowtie2.sh.log &
4.Sambamba去重PCR重复
cd /Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/mergeBam/
samtools sort ../align/AJV12-Ring1b.bam -o AJV12-Ring1b_sorted.bam
samtools sort ../align/JV38-Ring1b.bam -o JV38-Ring1b_sorted.bam
samtools index -M *_sorted.bam
sambamba markdup -r AJV12-Ring1b_sorted.bam AJV12-Ring1b_sorted_rm.bam
sambamba markdup -r JV38-Ring1b_sorted.bam JV38-Ring1b_sorted_rm.bam
#跑出来的结果
(seq)wu_lab@bio:align$ ls -lh *.bam
-rw-rw-r-- 1 wu_lab wu_lab 1.3G 12月 23 15:25 AJV12-Ring1b.bam
-rw-rw-r-- 1 wu_lab wu_lab 889M 12月 23 16:28 JV38-Ring1b.ba
(seq) wu_lab@bio:align$ ls -lh *.bam | cut -d' ' -f 5-
1.3G 12月 23 15:25 AJV12-Ring1b.bam
889M 12月 23 16:28 JV38-Ring1b.bam
查看比对率：
cd /Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/align/
#samtools index 命令为每个BAM文件创建一个索引文件（.bai 文件）
ls *.bam |xargs -i samtools index {}
#为每个BAM文件生成flagstat统计
ls *.bam | while read id ;do (nohup samtools flagstat $id > $(basename $id ".bam").stat & );done
cat *stat| grep 'N/A' *stat |grep '%'
#输出结果
AJV12-Ring1b.stat:22697954 + 0 mapped (100.00% : N/A)
AJV12-Ring1b.stat:22697954 + 0 primary mapped (100.00% : N/A)
AJV12-Ring1b.stat:22697954 + 0 properly paired (100.00% : N/A)
AJV12-Ring1b.stat:0 + 0 singletons (0.00% : N/A)
JV38-Ring1b.stat:13730006 + 0 mapped (100.00% : N/A)
JV38-Ring1b.stat:13730006 + 0 primary mapped (100.00% : N/A)
JV38-Ring1b.stat:13730006 + 0 properly paired (100.00% : N/A)
JV38-Ring1b.stat:0 + 0 singletons (0.00% : N/A)
4.peak calling
Peak Calling 分析 通常用于识别染色质免疫共沉淀测序（ChIP-Seq）、开放染色质测序（ATAC-Seq）或其他类似实验中的基因组区域（称为峰或peaks），这些区域通常代表了蛋白质-基因组相互作用位点或活跃的染色质区域。
这些峰是指在基因组特定区域中测序读数显著富集的区域，通过Peak Calling 方法，可以从背景噪声中提取这些生物学信号。
cd /Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/macs2/
macs3 callpeak -f BAMPE -t ../mergeBam/AJV12-Ring1b_sorted_rm.bam -g mm --nomodel -q 0.05 -n AJV12
macs3 callpeak -f BAMPE -t ../mergeBam/JV38-Ring1b_sorted_rm.bam -g mm --nomodel -q 0.05 -n JV38
#macs2 callpeak参数：
#-t 输入文件
#–control 对照组文件。
#-g 基因组大小。输入具体的数字（如1.0e+9或1000000000）。
#对于人、小鼠、线虫或果蝇，可以分别用hs、mm、ce或dm。
#–outdir 输出文件的路径
#-n 文件名前缀
#-q Q值
bamCoverage
inputdir=/Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/mergeBam/
outdir=/Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/bamCoverage/
ID="/Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/clean2/ID"
cat $ID | while read id; do
bamCoverage --binSize 10 -p 15 --normalizeUsing RPKM -b ${inputdir}/${id}_sorted_rm.bam -o ${outdir}${id}.bw
done
nohup bash bamCoverage.sh > bamCoverage.sh.log 2>&1 &
–binSize 10：设置 bin 的大小为 10 bp，这是计算覆盖度时使用的数据区块大小。
-p 15：指定并行处理的线程数为 15，提高处理速度。
–normalizeUsing RPKM：一共有RPKM,CPM,BPM,RPGC,None这几种方式，选择使用 RPKM 方法标准化数据。
-b *_sort.bam ：指定输入的 BAM 文件，根据提取的 ID 动态指定。
-o *.bw：指定输出的 BigWig 文件路径和名称。
deeptools可视化
多样本分析
multiBamSummary：这个工具用于生成多个BAM文件的覆盖度矩阵。它可以为每个bin（基因组中的一个区域）提供覆盖度统计信息，如平均覆盖度、最大覆盖度、中位数覆盖度等。这对于比较不同样本或不同条件下的基因组覆盖模式非常有用。
plotCorrelation：这个工具用于生成两个样本集之间的相关性散点图。它可以显示两个样本集在基因组覆盖度上的相似性或差异性。这对于评估样本间的一致性或发现技术变异很有帮助。
plotPCA：这个工具执行主成分分析（PCA），这是一种降维技术，可以将高维数据（如基因组覆盖度数据）转换为低维数据，同时尽可能保留数据的变异性。PCA图可以帮助识别样本间的群体结构和异常值。
BAM到BW：如果BAM文件已经被转换为BW格式，那可以使用 multiBigWigSummary 代替 multiBamSummary 来处理这些文件。multiBigWigSummary 是 deepTools 中用于处理BW文件的等效工具。
path="/home/lm/Z_Projects/chipseq"
mkdir deeptools
# 统计reads在全基因组范围的情况
# 双端测序
#multiBamSummary bins -bs 1000 --bamfiles ${path}/intersect/*_sort.bam --extendReads 130 -out treat_results.npz
# 单端测序
nohup multiBamSummary bins -bs 1000 --bamfiles ${path}/intersect/*_sort.bam -out ./deeptools/treat_results.npz > multibam.log &
# 相关性散点图
plotCorrelation -in ./deeptools/treat_results.npz -o ./deeptools/treat_results.pdf --corMethod spearman -p scatterplot
# 热图
plotCorrelation -in ./deeptools/treat_results.npz -o ./deeptools/treat_results_heatmap.pdf --corMethod spearman -p heatmap
# 主成分分析
plotPCA -in ./deeptools/treat_results.npz -o ./deeptools/pca.pdf
从UCSC下载TSS的bed文件：
#先进入目录解压缩bed文件
cd /Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/tss/
gunzip -d mm39.tss.bed.gz
#批量处理-绘制TSS附近3k区域
path="/Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/"
bed="/Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/tss/mm39.tss.bed"
for id in ${path}/bamCoverage/*.bw ;
do
echo $id
file=$(basename $id)
sample=${file%.*}
echo $sample
computeMatrix reference-point --referencePoint TSS -p 15 \
-b 3000 -a 3000 \
-R $bed \
-S $id \
--skipZeros -o ${path}/tss/${sample}_TSS_3K.gz \
--outFileSortedRegions ${path}/tss/${sample}_TSS_3K.bed
## both plotHeatmap and plotProfile will use the output from computeMatrix
plotHeatmap -m ${path}/tss/${sample}_TSS_3K.gz -out ${path}/tss/${sample}_Heatmap_3K.png
plotHeatmap -m ${path}/tss/${sample}_TSS_3K.gz -out ${path}/tss/${sample}_Heatmap_3K.pdf --plotFileFormat pdf --dpi 720
plotProfile -m ${path}/tss/${sample}_TSS_3K.gz -out ${path}/tss/${sample}_Profile_3K.png
plotProfile -m ${path}/tss/${sample}_TSS_3K.gz -out ${path}/tss/${sample}_Profile_3K.pdf --plotFileFormat pdf --perGroup --dpi 720
done
nohup bash computeMatrix-3k.sh > 3k.log 2>&1 &
1.computeMatrix reference-point 命令
这个命令用于生成一个矩阵文件，该文件包含了以转录起始位点（TSS）为中心的检测的蛋白或者组蛋白修饰结合信号的分布情况。参数解释如下：
–referencePoint TSS: 指定参考点为转录起始位点。
-p 15: 使用 15 个处理器并行计算。
-b 10000: 指定参考点上游的区域长度为 10000 个碱基对。
-a 10000: 指定参考点下游的区域长度为 10000 个碱基对。
-R : 指定包含转录起始位点的 BED 文件。
-S : 指定包含蛋白结合信号的 bigWig 文件。
–skipZeros: 跳过信号值为零的区域。
-o *_TSS.gz: 指定输出矩阵文件的名称。
–outFileSortedRegions *_genes.bed: 指定输出排序后的区域列表文件的名称。
2. plotHeatmap 命令
这个命令用于生成热图，展示基因组区域的信号强度。参数解释如下：
-m *_TSS.gz: 指定输入的矩阵文件。
-out test.png: 指定输出的热图文件名，格式为 PNG。
–plotFileFormat pdf: 指定输出文件的格式为 PDF。
–dpi 720: 指定输出图像的分辨率为 720 DPI。
3.plotProfile 命令
这个命令用于生成信号的平均分布图，展示基因组区域的信号强度随位置的变化。参数解释如下：
-m *_TSS.gz: 指定输入的矩阵文件。
-out test.png: 指定输出的轮廓图文件名，格式为 PNG。
–plotFileFormat pdf: 指定输出文件的格式为 PDF。
–perGroup: 为每个组（如每个基因或每个样本）生成单独的轮廓图。
–dpi 720: 指定输出图像的分辨率为 720 DPI。
R语言注释peak
#用MACS2的bed文件,bed/peak文件第一列需要以chr开头，已经制作符合规格的bed格式文件
cd /Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/macs2/
#先查看bed文件
#给第一列都加上chr
find /Data/wu_lab/zyy/2022-jv-ajv-ring1b-chip-no-sp/macs2/ -type f -name "*.bed" | while read bed_file; do
new_file="${bed_file%.bed}_modified.bed"
awk -v OFS="\t" '{if ($1 ~ /^[0-9]+$/) {$1="chr"$1}; print}' "$bed_file" > "$new_file"; done
#下载制作好的bed文件，进入R进行下游可视化分析
##1.#加载R包,指定参考基因组
library(ChIPseeker)
library(ChIPpeakAnno)
library(clusterProfiler)
library("org.Hs.eg.db")
library("org.Mm.eg.db")
library(ggplot2)
library(plotrix)
library(TxDb.Mmusculus.UCSC.mm39.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm39.knownGene
txdb
##2.读入peak文件
AJV12 <- readPeakFile("AJV12_summits_modified.bed")
keepChr <- !grepl("_",seqlevels(AJV12)) #去除带有“_”的染色体
seqlevels(AJV12,pruning.mode="coarse") <- seqlevels(AJV12)[keepChr]
##3.以tss上游2000下游500的区域作为promoter区域进行区域注释
peakAnno <- annotatePeak(AJV12, tssRegion=c(-2000, 500),
TxDb=txdb, annoDb="org.Mm.eg.db")
##4.绘制区域位置饼图
plotAnnoPie(peakAnno)
plotAnnoBar(peakAnno)
vennpie(peakAnno)
upsetplot(peakAnno)
plotDistToTSS(peakAnno,title="Distribution of transcription factor-binding loci \n relative to TSS")
##5.查看注释表格
peak.df=as.data.frame(peakAnno)
##6.功能富集分析
提取peakAnnolist中的基因，结合clusterProfiler包对peaks内的邻近基因进行富集注释。
可视化
4.1 IGV
#示例代码
bamCoverage --bam * --outFileName * --blackListFileName --verbose -- normalizeUsing * #### 分类后的bam文件 输出文件名 显示进程 标化方法，为了防止测序深度不同，可根据help文档自选
参数的解释：
--bam *：
这个参数指定了要分析的BAM文件。* 是一个通配符，表示当前目录下的所有BAM文件
--outFileName *：
这个参数指定输出文件的名称。* 代表每个输入文件的基本名称，输出文件将以这个名称为基础进行命名。
--blackListFileName：
这个参数后面应该跟一个文件名，该文件包含要排除的区域（比如富含重复元素的区域）。这些区域在分析覆盖度时会被忽略。由于您没有提供具体的文件名，所以这部分需要您根据实际情况补充。
--verbose：
这个参数用于增加程序运行时的输出信息量，使得用户能够看到更多的处理细节和进度。
--normalizeUsing *：
这个参数用于指定标准化方法，以校正不同样本之间的测序深度差异。* 应该被替换为具体的标准化方法，例如 RPKM、RPGC 或 RAW 等。
deepTools 文档提供了不同标准化方法的详细说明和适用场景。
RAW：不进行标准化，直接使用原始的reads计数。
RPGC：标准化到每个基因组区域的reads数除以该区域的长度，再乘以基因组的大小。
RPKM：标准化到每千个碱基的reads数乘以百万个reads。
为了防止由于不同样本之间的测序深度不同而导致的比较偏差，选择合适的标准化方法是非常重要的。
#实例数据：
bamCoverage --bam JV38-Ring1b_sorted_rm.bam --outFileName JV38-Ring1b --verbose --normalizeUsing RPKM
bamCoverage --bam AJV12-Ring1b_sorted_rm.bam --outFileName AJV12-Ring1b --verbose --normalizeUsing RPKM
全部评论 (0)
还没有任何评论哟~