ChIP-seq分析
創(chuàng)建小環(huán)境
conda create -n ChIP python=3
安裝軟件
conda install -y sra-tools
conda install -y trim-galore samtools
conda install -y deeptools homer meme
conda install -y macs2 bowtie bowtie2 bwa
數(shù)據(jù)下載
prefetch SRR8476707 -O ~
數(shù)據(jù)轉換，sra格式轉換為fastq格式
for i in *sra
do
echo $i
/path/sratoolkit.2.3.5-2-mac64/bin/fastq-dump --split-3 $i
done
|單個命令如下
fastq-dump --gzip --split-3 -O ./ ./SRR8476707.sra
過濾，用trim-galore軟件，直接fastqc
因為文件多，所以寫個腳本直接運行
vim trim_galore_batch.sh
輸入以下內(nèi)容
#!/bin/bash
# This is for trimming a batch data
for i in SRR8476709 SRR8476710 SRR8476717 SRR8476718 SRR8476719 SRR8476720 SRR8476721
do
trim_galore -q 25 --phred33 --length 25 --stringency 4 --gzip --paired ~/Data/H3K9me2-ChIP_seq_VC/SRR/fastq/${i}_1.fastq.gz ~/Data/H3K9me2-ChIP_seq_VC/SRR/fastq/${i}_2.fastq.gz -o ~/Data/H3K9me2-ChIP_seq_VC/clean_2 --fastqc
done
運行
bash trim_galore_batch.sh
建索引
bwa index -a bwtsw ~/genome/mouse/ensembl/Mus_musculus.GRCm39.dna.toplevel.fa
| 索引和gtf&gfa文件建立在一個文件夾下面，會生成5個文件，分別是.amb .ann .bwt .pac .sa
比對
比對一般用BWA或者Bowtie2
|也可以寫成文件，運行文件，當時一直報錯，就一個一個比對的
加上nohup ***** &，讓代碼在后臺運行
nohup bwa mem -v 3 -t 4 ~/genome/mouse/ensembl/Mus_musculus.GRCm39.dna.toplevel.fa ~/Data/H3K9me2-ChIP_seq_VC/fastq_p/SRR8476710_1_val_1.fq ~/Data/H3K9me2-ChIP_seq_VC/fastq_p/SRR8476710_2_val_2.fq -o ~/Data/H3K9me2-ChIP_seq_VC/bwa/SRR8476710_bwa.sam &
轉換為bam文件及排序，或者有時候需要合并bam文件（將多個重復合并到一個）
|因為bwa生成的是sam文件，所以要轉換成bam文件以及進行排序（如果不排序后面用macs2 peak calling的時候可能會報錯）
vim bam.sh
輸入以下內(nèi)容
for i in SRR8476707 SRR8476708 SRR8476709 SRR8476710 SRR8476711 SRR8476717 SRR8476718 SRR8476719 SRR8476720 SRR8476721
do
samtools view -bS -h ~/Data/H3K9me2-ChIP_seq_VC/bwa/${i}_bwa.sam -o ~/Data/H3K9me2-ChIP_seq_VC/bam/${i}_bwa.bam
查看比對成功率
samtools flagstat -@ 8 SRR8476707.sam
轉為bam文件
vim bam.sh
輸入以下內(nèi)容
do
samtools view -bS -h ~/Data/H3K9me2-ChIP_seq_VC/bwa_NCBI/${i}_bwa.sam -o ~/Data/H3K9me2-ChIP_seq_VC/bam_NCBI/${i}_bwa.bam
done
bash bam.sh
對bam文件進行排序
vim bam.sort.sh
for i in SRR8476707 SRR8476708 SRR8476709 SRR8476710 SRR8476711 SRR8476717 SRR8476718 SRR8476719 SRR8476720 SRR8476721
do
samtools sort -@ 5 -o ~/Data/H3K9me2-ChIP_seq_VC/bam_NCBI/${i}_bwa.sort.bam ~/Data/H3K9me2-ChIP_seq_VC/bam_NCBI/${i}_bwa.bam
done
轉換為bw模式，導入IGV查看
安裝deeptools
conda install -c bioconda/label/cf201901 deeptools
在sorted bam文件下進行轉換
ls *.bam |xargs -i samtools index {}
ls *.bam |while read id;do \
bamCoverage --normalizeUsing CPM -b $id -o ${id%%.*}.bw & \
done
peakcalling
vim macs2_2.sh
輸入以下內(nèi)容
for i in SRR8476717 SRR8476718 SRR8476719 SRR8476720 SRR8476721
do
macs2 callpeak -f BAM -t ${i}_bwa.sort.bam -n ${i}_bwa.sort.bam -g mm --outdir ~/Data/H3K9me2-ChIP_seq_VC/mcs2_2 --bdg -q 0.05
done
運行代碼
nohup bash macs2_2.sh>macs_2.log &
看下找到了多少個peak
wc -l *bed
有3個生物學重復的組蛋白修飾ChIP-Seq數(shù)據(jù)怎樣合并比較好呢？文獻中目前有2種方法：（1）直接合并bam文件，再call peak；（2）3個重復分別call peak，然后取至少在2個樣品中同時出現(xiàn)的 peak 進行后續(xù)的分析。這兩種方法各有哪些優(yōu)缺點？哪種比較科學呢？期待老師的解答，謝謝！
callpeak會得到如下結果文件：
NAME_summits.bed：Browser Extensible Data，記錄每個peak的peak summits，話句話說就是記錄極值點的位置。MACS建議用該文件尋找結合位點的motif。能夠直接載入UCSC browser，用其他軟件分析時需要去掉第一行。
NAME_peaks.xls：以表格形式存放peak信息，雖然后綴是xls，但其實能用文本編輯器打開，和bed格式類似，但是以1為基，而bed文件是以0為基.也就是說xls的坐標都要減一才是bed文件的坐標。
NAME_peaks.narrowPeak，NAME_peaks.broadPeak 類似。后面4列表示為， integer score for display， fold-change，-log10pvalue，-log10qvalue，relative summit position to peak start。內(nèi)容和NAME_peaks.xls基本一致，適合用于導入R進行分析。
NAME_model.r：能通過$ Rscript NAME_model.r作圖，得到是基于你提供數(shù)據(jù)的peak模型。
.bdg：能夠用 UCSC genome browser 轉換成更小的 bigWig 文件。
————————————————
版權聲明：本文為CSDN博主「fengzhibifang」的原創(chuàng)文章，遵循CC 4.0 BY-SA版權協(xié)議，轉載請附上原文出處鏈接及本聲明。
原文鏈接：https://blog.csdn.net/weixin_41745858/article/details/82628775
