ChIP-seq分析全记录

一、背景知识

CHIP-seq染色体免疫共沉淀技术。是研究==细胞内DNA与蛋白质互作==就是将chip实验得到的dna拿来测序,再进行下游分析。关键在于Chip实验。

ChIP被用来研究细胞内DNA与蛋白质相互作用,具体来说就是确定特定蛋白(如转录因子)是否结合特定基因组区域(如启动子或其它DNA结合位点)——可能定义顺反组。ChIP还被用来确定基因组上与组蛋白修饰相关的特定位点

用于在全基因组范围内研究DNA结合蛋白(相互反应)、组蛋白修饰(表观遗传标记)与 核小体 的技术。有助于了解基因之间的相互调控以及染色体的功能结构。

二、CHIP 实验

  • CHIP 实验流程分为:

    1. 蛋白交联到DNA上,保证蛋白和DNA能结合,找到互作位点。
    2. 超声波剪切DNA链
    3. 加上 带有抗体的磁珠用于 免疫沉淀靶蛋白。抗体很重要。
    4. 接触蛋白交联,纯化DNA
    5. 送测序,后续分析。
  • 实验设计的关键:

    1. 抗体质量、
    2. 样本量(10-50ng DNA,PCR扩增越少,偏差越少)
    3. 空白对照:(排除假阳性的存在,开放染色体区域、重复序列) input DNA对照???mock IP DNA(免疫共沉淀不含有抗体的DNA)

三、CHIP 分析流程:

  • CHIP-SEQ 分析流程,分析分为4步

    • 质量控制,用的是Fastqc等
    • 序列比对,Bowtie2或BWA
    • peak calling, MACS
    • peak注释, ChIPseeker
  • CHIP数据分析所特有的步骤:peak calling

    1. 染色体上信号波形的定义;
    2. 建立背景矫正模型;
    3. 建立搜索peaks的准则,即建立判断怎样可以是一个peak(一般会用到背景矫正模型中的背景值);
    4. 矫正模型,过滤掉假阳性的peaks;
    5. 给找到的peaks排序,给出显著度。

read只是跟随着TF一起沉淀下来的DNA fragment的末端,read的位置并不是真实的TF结合的位置。所以在peak-calling之前,延伸read是必须的。不同TF大小不一样,对read延伸的长度也理应不同。我们知道,测得的read最终其实会近似地平均分配到正负链上,这样,对于一个TF结合热点而言,read在附近正负链上会近似地形成“双峰”。MACS会以某个window size扫描基因组,统计每个window里面read的富集程度,然后抽取(比如1000个)合适的(read富集程度适中,过少,无法建立模型,过大,可能反映的只是某种偏好性)window作样本,建立“双峰模型”。最后,两个峰之间的距离就被认为是TF的长度D,每个read将延伸D/2的长度。

四、具体代码记录

目标:获取感兴趣的基因如ARF3在已发表的chip-seq数据AP3/PI/SEP3中是否存在相同的峰(位置,高度)

1. 文章数据下载

  1. SEP3: Target Genes of the MADS Transcription
    Factor SEPALLATA3: Integration of
    Developmental and Hormonal Pathways
    in the Arabidopsis Flower GSE:GSE14600
  2. AP3-PI:Molecular basis for the specification of floral organs by
    APETALA3 and PISTILLATA GSE:GSE38358
  3. AG : Control of Reproductive Floral Organ Identity Specification in
    Arabidopsis by the C Function Regulator AGAMOUS GSE :GSE45938
  4. AP1: Orchestration of Floral Initiation
    by APETALA1
  5. AP1/SEP3 GSE46986
  • 根据 GSE号 以及 PRJNA 号下载
### 以AG数据为例
tail -n +2 AG-SraRunTable.txt |cut -f 7 |xargs -i nohup prefetch {} \&
  • sra 数据转换fastq.gz,注意也可以不用--gzip参数
### 测序数据皆为单端测序,故不需要--split-3

for i in `ls *.sra`
do
nohup fastq-dump.2 --gzip -O 1_fastq/ $i &
done

2. 数据的质控

2.1 fastqc 处理

### 直接fastqc显示测序质量
ls *.fq.gz |parallel fastqc -o ./ --nogroup {} &

### MultiQC批量显示测序的质量

##在已经用fastqc得到测序的html文件和zip文件夹
multiqc *fastqc.zip --pdf

2.2 fastp处理

###fastp快速质控
fastp -i SRR502857.fastq.gz -o SRR502857_qc.fq.gz
### 出现adapter sequences的自动检测出现错误的问题,加上-A参数


###根据fastqc报告查看接头序列指定接头的序列
cat ~/miniconda3/envs/bioinfo/opt/fastqc-0.11.5/Configuration/adapter_list.txt
echo ">nextera" > nextera.fa
echo "CTGTCTCTTATACACATCTCCGAGCCCACGAGAC" >> nextera.fa


### 2018/7/10 -W sliding window -M 20

for i in `ls *.gz`; do i=${i%.fastq.gz}; 
fastp -3 -W 4 -M 20 -a AGATCGGAAGAG -i $i.fastq.gz -o ../clean_data/fastp_out_2.0/$i.qc.fq.gz & done
### 结果发现测序质量十分不错。

3. 序列的比对

### 构建索引 基本命令就是bowtie2-build 基因组序列.fa 索引名字
bowtie2-build Arabidopsis_thaliana.TAIR10.dna.toplevel.fa ath_index_bt1 
### 比对
for i in `ls *.fq.gz`
do
i=${i%.fq.gz}
bowtie -p 10 --best --chunkmbs 320 Database/Ath/ath_index_bt1 -q $i.fq.gz -S $i.sam
done


### 2018/7/10 unique mapping对于bowtie1 需要参数-m 1 
for i in `ls *.gz`; 
do 
i=${i%.qc.fq.gz}; 
nohup bowtie -p 5 -m 1 Database/plant/arabidopsis_th/ath_index_bt1 -q $i.qc.fq.gz -S ../3_mapping/unique_mapping/$i.sam >$i.log2 & 
done

新浪博客——关于map当中的unique mapped reads问题

  • samtools faidx Arabidopsis_thaliana.TAIR10.dna.toplevel.fa 1:14331064-14331662

4. Samtools转换

samtools view -bS SRR851694.sam | samtools sort -@ 4 - -T $i -o SRR851694.sorted.bam &
samtools index SRR851694.sorted.bam

5. BAM 转 bigwig文件 后续在IGV里查看各个基因的峰

bamCoverage -b AG_IP.sorted.bam --normalizeUsing RPKM --binSize 30 --smoothLength 300 -p 10 --extendReads 200 -o AG_IP.bw

6. call peaks 利用macs2软件,需conda下的pytho2.7环境。

### macs 较为严格的call peaks
macs -t SRR502859.sorted.bam -c SRR502860.sorted.bam -n AP3 -g 1.2e8 -p 1e-5 -O ../peaks_call/



### macs2 较为宽松的call peaks. 
macs2 callpeak -t SRR016811.sorted.bam -c SRR016812.sorted.bam -B --nomodel -g dm --keep-dup 1 -n SEP3-rep2 --trackline --SPMR --call-summits

### macs2 2014年AP1/SEP3的Genome biology文章方法里所提供的macs2参数改进
macs2 callpeak -t ap1_rep1.sorted.bam -c control.sorted.bam -g dm --mfold 2 20 -q 0.001 -n ap1_rep1 --outdir macs2_new/ &


7. 不同peaks的交集并集处理

### bedtools处理
bedtools intersect -a ap1_peaks.narrowpeaks -b ag_peaks.narrowpeaks -wa |wc -l

### bedops 处理
bedops --intersect ap1_peaks ap3_peaks pi_peaks sep3_peaks > petal_peaks.bed


8. 注释Rstudio中用ChIPseeker

library(ChIPseeker)
library(clusterprofiler)
library(TxDb.Athaliana.BioMart.plantsmart28)
library(org.At.tair.db)
stamen_peak <- readPeakFile("stamen_intersect.bed")

## 注释
stamen_peak_anno <- annotatePeak(peak = stamen_peak,tssRegion = c(-3000,3000),TxDb = TxDb.Athaliana.BioMart.plantsmart28,level="gene")
stamen_gene <- as.data.frame(stamen_peak_anno)

###利用lappy集中注释
peak<- list(sepal_peak_RE,petal_peak_RE,stamen_peak_RE,carpel_peak_RE)
peak_anno <- lapply(peak,annotatePeak,tssRegion = c(-3000,3000),TxDb = TxDb.Athaliana.BioMart.plantsmart28,level = "gene")
sepal_peak_anno_RE <- as.data.frame(peak_anno[[1]])
plotDistToTSS(peak_anno)


###转id
keytypes(org.At.tair.db)
sepal_gene_test<- select(org.At.tair.db,keys = sepal_gene$geneId,columns = c("SYMBOL","GENENAME"),keytype = "TAIR")

### clusterprofiler功能注释
stamen_mf <- enrichGO(gene = stamen_specific$stamen_specific_PEvsST,OrgDb = org.At.tair.db,ont = "MF",keyType = "TAIR",pvalueCutoff = 0.05)
dotplot(stamen_cc)



其他

  • 从染色体指定区段获取序列,并进行CArG-box的定位
samtools faidx /Database/plant/a_th/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa 1:25686061-25687376|seqkit locate -i -d -p CHWWWWWWDG|tail -n +2|cut -f 1,5,6,7|perl -lne '@F=split/[\:\-\t]/;print $F[5],"\t" ,$F[1]+$F[3]'
  • BAM 转 bigwig文件
bamCoverage -b AG_IP.sorted.bam --normalizeUsing RPKM --binSize 30 --smoothLength 300 -p 10 --extendReads 200 -o AG_IP.bw

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 156,630评论 4 359
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 66,405评论 1 289
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 106,382评论 0 237
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 43,548评论 0 203
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 51,853评论 3 285
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 40,276评论 1 209
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,638评论 2 309
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 30,350评论 0 195
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 34,015评论 1 238
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,309评论 2 240
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,836评论 1 256
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 28,206评论 2 251
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,767评论 3 231
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 25,972评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,715评论 0 192
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 35,319评论 2 269
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 35,222评论 2 258

推荐阅读更多精彩内容