[Python知识库]Chip-seq分析笔记

前言

自己做笔记自己看的

一、软件安装

? 设备：mac m1 电脑

• 软件程序安装
• 参考教程Macbook Pro M1芯片Python开发环境配置_朝歌晚酒南栀雪的博客-CSDN博客
  • xcode
  • item2
  • git
  • homebrew

    brew install wget #安装 wget
    

  • anaconda-pkg

?? ? ? ? ? ? ? ? ? ???????Anaconda | Individual Edition??下载链接

? ? ? ? ? ? ? ? ? ?Installing on macOS — Anaconda documentation?帮助链接

? ? ? ? ? ? ? ?pycharm CE

二、创建环境及安装软件

1.创建环境 chipseq

conda? create -n chipseq? python=3
	conda activate chipseq #激活环境
	conda deactivate

2.chipseq 环境下安装软件

ps:一项一项安装，不然容易卡住

conda install -y bioconda parallel-fastq-dump ?#SRR→FASTQ 转换

conda install -y trim-galore #质控

conda install -y?bioconda bowtie2 # fastq比对到 hg19

brew tap homebrew/science?

brew install samtools?#sam 转化为 bam

conda install -y macs2 #峰值定量，差异分析

conda install -y conda-forge libgfortran #macs2 差异分析辅助用

conda install -y bioconda deeptools #可视化
   #包含 bamcoverage
    https://deeptools.readthedocs.io/en/latest/content/tools/bamCoverage.html

IGV 下载https://software.broadinstitute.org/software/igv/download #可视化

ROSE 安装

https://bitbucket.org/young_computation/rose/src/master/ ROSE 下载
#把所有要处理的文件都放到 rose 文件夹：sorted.bam, bed,bam.bai
#创建 Python2.7
	conda  create -n rose  python=2.7
	conda activate rose

Homer安装

• 公众号参考?HOMER | chipseq数据进行peaks的差异分析
• homer安装：?chipseq 环境下，conda install -c bioconda homer
  • 下载configureHomer.pl ?http://homer.ucsd.edu/homer/configureHomer.pl
  • 文件放在想要安装 homer?的文件夹中，/Users/xusiqi/chip/homer
    • 目录文件夹中安装?homer :?

perl configureHomer.pl -install

open -a TextEdit ~/.bash_profile

PATH=$PATH:/Users/chucknorris/homer/bin/??#路径加入到系统路径中

source ~/.bash_profile

• perl?configureHomer.pl -list #查看所有的?list

三、具体分析步骤

1.数据来源

使用数据Wang J, Zou JX, Xue X, Cai D et al.?ROR-γ drives androgen receptor expression and represents a therapeutic target in castration-resistant prostate cancer.?Nat Med?2016 May;22(5):488-96. PMID:?27019329

2.下载数据并重命名

wget https://sra-pub-run-odp.s3.amazonaws.com/sra/SRR8557353/SRR8557353
mv???? SRR8557353????? input.h3k.con
?
wget https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos3/sra-pub-run-20/SRR8557351/SRR8557351.1
mv? SRR8557351.1?? ip.h3k.con
?
wget https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos3/sra-pub-run-21/SRR8557354/SRR8557354.1
mv? SRR8557354.1?? input.h3k.xy
?
wget https://sra-pub-run-odp.s3.amazonaws.com/sra/SRR8557352/SRR8557352
mv?? SRR8557352???? ip.h3k.xy

3.fastq 文件转换（--sra-id?此步骤把 SRR 的名称改掉，加快运行速度，速度非常快）

parallel-fastq-dump —sra-id input.h3k.con —threads 35? —outdir out/ —split-files —gzip
parallel-fastq-dump —sra-id ip.h3k.con?? —threads 35? —outdir out/ —split-files —gzip
parallel-fastq-dump —sra-id? input.h3k.xy?? —threads 35? —outdir out/ —split-files —gzip
parallel-fastq-dump —sra-id?? ip.h3k.xy?? —threads 35? —outdir out/ —split-files —gzip

4.质控

4.1 trim_galore

• trim_galore（只能单线程）
  • Cutadapt:去接头，去除3端低质量碱基，去除长度太短的序列
    • ????report.txt
    • ????fq.gz(用来比对)

trim_galore input.h3k.con_1.fastq.gz? -q 25 —phred33 —length 25 -e 0.1 —stringency 4 
trim_galore input.h3k.xy_1.fastq.gz? -q 25 —phred33 —length 25 -e 0.1 —stringency 4 
trim_galore ip.h3k.con_1.fastq.gz? -q 25 —phred33 —length 25 -e 0.1 —stringency 4 
trim_galore ip.h3k.xy_1.fastq.gz? -q 25 —phred33 —length 25 -e 0.1 —stringency 4 

4.2 FastQC质控报告生成（旧的数据）

• fastqc -o?/Users/xusiqi/CHIP/out?-t 6?/Users/xusiqi/CHIP/out/SRR2242680_1_trimmed.fq.gz
• fastqc -o?/Users/xusiqi/CHIP/out?-t 6?/Users/xusiqi/CHIP/out/SRR2242690_1_trimmed.fq.gz
• fastqc -o?输出绝对路径?-t 线程数 输入文件绝对路径
• -o?后面为文件输出绝对路径，-t 6(为线程数)，?最后为输入文件绝对路径；大约 5min/file

5.bowtie2比对(15min/file)

bowtie2 -p 35 -x /Users/xusiqi/CHIP/index/hg19/hg19 -U? /Users/xusiqi/CHIP/lesson23/out/input.h3k.con_1_trimmed.fq.gz? -S input.h3k.con.sam
bowtie2 -p 35 -x /Users/xusiqi/CHIP/index/hg19/hg19 -U? /Users/xusiqi/CHIP/lesson23/out/input.h3k.xy_1_trimmed.fq.gz? -S input.h3k.xy.sam
bowtie2 -p 35 -x /Users/xusiqi/CHIP/index/hg19/hg19 -U? /Users/xusiqi/CHIP/lesson23/out/ip.h3k.con_1_trimmed.fq.gz? -S ip.h3k.con.sam
bowtie2 -p 35 -x /Users/xusiqi/CHIP/index/hg19/hg19 -U? /Users/xusiqi/CHIP/lesson23/out/ip.h3k.xy_1_trimmed.fq.gz? -S ip.h3k.xy.sam

6.samtools sam?转换?bam

6.1view &sort

samtools view -S -b input.h3k.con.sam > input.h3k.con.bam
samtools sort input.h3k.con.bam -o input.h3k.con.sorted.bam
?
samtools view -S -b input.h3k.xy.sam > input.h3k.xy.bam
samtools sort input.h3k.xy.bam -o input.h3k.xy.sorted.bam
?
samtools view -S -b ip.h3k.con.sam >ip.h3k.con.bam
samtools sort ip.h3k.con.bam -o ip.h3k.con.sorted.bam
?
samtools view -S -b ip.h3k.xy.sam > ip.h3k.xy.bam
samtools sort ip.h3k.xy.bam -o ip.h3k.xy.sorted.bam

6.2 index

samtools index input.h3k.con.sorted.bam
samtools index input.h3k.xy.sorted.bam
samtools index ip.h3k.con.sorted.bam
samtools index? ip.h3k.xy.sorted.bam

7.bamCoverage?BW转换 bam 转换成 bw

bamCoverage -e 170 -bs 10 -p 35 -b input.h3k.con.sorted.bam -o input.h3k.con.sorted.bw
bamCoverage -e 170 -bs 10 -p 35 -b input.h3k.xy.sorted.bam -o input.h3k.xy.sorted.bw
bamCoverage -e 170 -bs 10 -p 35 -b ip.h3k.con.sorted.bam -o ip.h3k.con.sorted.bw
bamCoverage -e 170 -bs 10 -p 35 -b ip.h3k.xy.sorted.bam -o ip.h3k.xy.sorted.bw

8.IGV查看 bw文件 peak

• IGV?下载?Downloads | Integrative Genomics Viewer
  • IGV?打开?bw?文件

9.MACS2 callpeak

• macs2 callpeak -t?实验?sorted.bam?-c???对照inputsorted.bam?-g hs -B -f BAM -n?输出名称（无后缀，简单命名）?-q 0.05

  macs2 callpeak -t ip.h3k.con.sorted.bam -c? ?input.h3k.con.sorted.bam -g hs -B -f BAM -n con? -q 0.05
  macs2 callpeak -t ip.h3k.xy.sorted.bam -c? ?input.h3k.xy.sorted.bam -g hs -B -f BAM -n xy -q 0.05
  

10.MACS2 差异分析（2021-10-15 待完善）

参考教程：

Call differential binding events · macs3-project/MACS Wiki · GitHub

11.Homer 分析 motif

公众号参考?HOMER | chipseq数据进行peaks的差异分析

11.1查看所有的?list:

perl configureHomer.pl -list

11.2下载注释信息

perl configureHomer.pl -install hg19   
#（安装hg19的GENOMES,会自动下载human数据，1.38G，多试几次，网速可达 6-7Mb)

11.3 Create a tag directory

• homer?下新建一个?out?文件夹，在此路径下运行以下代码

  makeTagDirectory?文件夹名（tag directory）??-genome?hg19(带绝对路径)??-checkGC??输入文件 sam（绝对路径）

• 输入文件名：input.h3k.con.sam?input.h3k.xy.sam??ip.h3k.con.sam????ip.h3k.xy.sam

• makeTagDirectory? input.h3k.con ?-genome /Users/xusiqi/chip/homer/data/genomes/hg19 -checkGC /Users/xusiqi/chip/lesson23/out/input.h3k.con.sam
  makeTagDirectory input.h3k.xy ?-genome /Users/xusiqi/chip/homer/data/genomes/hg19 -checkGC /Users/xusiqi/chip/lesson23/out/input.h3k.xy.sam
  makeTagDirectory? ip.h3k.con ?-genome /Users/xusiqi/chip/homer/data/genomes/hg19 -checkGC /Users/xusiqi/chip/lesson23/out/ip.h3k.con.sam
  makeTagDirectory? ip.h3k.xy ?-genome /Users/xusiqi/chip/homer/data/genomes/hg19 -checkGC /Users/xusiqi/chip/lesson23/out/ip.h3k.xy.sam
  

11.4 Differentially Bound Peaks（找差异peak）

getDifferentialPeaks <peak file> <target tag directory> <background tag directory> [options]

输入文件名：

con_peaks.narrowPeak ?xy_peaks.narrowPeak
input.h3k.con ? input.h3k.xy ? ?ip.h3k.con ? ? ip.h3k.xy

getDifferentialPeaks /Users/xusiqi/chip/lesson23/out/con_peaks.narrowPeak ip.h3k.con/ input.h3k.con/ > con.diffpeaks.csv
getDifferentialPeaks /Users/xusiqi/chip/lesson23/out/xy_peaks.narrowPeak ip.h3k.xy/ input.h3k.xy/ > xy.diffpeaks.csv

11.5?Annotate peaks（Peak注释）

Usage:?annotatePeaks.pl?<peak?file?|?tss>?<genome?version>??[additional?options...]

输入文件名：con.diffpeaks.csv ? ?xy.diffpeaks.csv

annotatePeaks.pl con.diffpeaks.csv /Users/xusiqi/chip/homer/data/genomes/hg19 > con.diffpeaks.anno.txt
annotatePeaks.pl xy.diffpeaks.csv /Users/xusiqi/chip/homer/data/genomes/hg19 > xy.diffpeaks.anno.txt

用 R 把差异peak 和 Peak 注释合并：fread,merge 函数

11.6 其他功能

• annotatePeaks.pl程序还可以用于创建显示相对于给定基因组特征(包括转录起始位点(TSS)或用户想要定义的任何其他区域)的相对读富集的直方图。由于TSS经常用于此目的，所以HOMER为TSS提供了一个内置注释(基于RefSeq转录本)。创建直方图的关键参数是“-hist #”和“-size #”选项，它们控制直方图的装箱大小和总长度。另一个重要的选项是“-d”，它指定要为哪些实验编译直方图。
  • annotatePeaks.pl tss hg19 -size 8000 -hist 10 -d h3k.con/ h3k.bay/ > output.txt
  • 用电子表格程序/Excel打开output.txt文件。将注意到第一列给出了到TSS的距离偏移量，然后是对应于每个实验的“覆盖率”、“+标记”和“-标记”的列。尝试用第一列作为X-Y线形图来查看模式。

12. ROSE 寻找超级增强子 SE（super enhancer）

12.1 参考教程：?ROSE | 超级增强子鉴定

12.2 ROSE 安装见上

12.3 ROSE 寻找 SE

?
python ROSE_main.py -g HG19 -i xy_summits.bed -r ip.h3k.xy.sorted.bam -c input.h3k.xy.sorted.bam -o /Users/xusiqi/chip/rose/xy
python ROSE_main.py -g HG19 -i con_summits.bed -r ip.h3k.con.sorted.bam -c input.h3k.con.sorted.bam -o /Users/xusiqi/chip/rose/con

（未完待续）