一蜒蕾、背景知識(shí)
CHIP-seq染色體免疫共沉淀技術(shù)。是研究==細(xì)胞內(nèi)DNA與蛋白質(zhì)互作==就是將chip實(shí)驗(yàn)得到的dna拿來測(cè)序焕阿，再進(jìn)行下游分析咪啡。關(guān)鍵在于Chip實(shí)驗(yàn)。
ChIP被用來研究細(xì)胞內(nèi)DNA與蛋白質(zhì)相互作用暮屡，具體來說就是確定特定蛋白（如轉(zhuǎn)錄因子）是否結(jié)合特定基因組區(qū)域（如啟動(dòng)子或其它DNA結(jié)合位點(diǎn)）——可能定義順反組撤摸。ChIP還被用來確定基因組上與組蛋白修飾相關(guān)的特定位點(diǎn)
用于在全基因組范圍內(nèi)研究DNA結(jié)合蛋白（相互反應(yīng)）、組蛋白修飾（表觀遺傳標(biāo)記）與 核小體 的技術(shù)。有助于了解基因之間的相互調(diào)控以及染色體的功能結(jié)構(gòu)准夷。
二钥飞、CHIP 實(shí)驗(yàn)
-
CHIP 實(shí)驗(yàn)流程分為：
- 蛋白交聯(lián)到DNA上，保證蛋白和DNA能結(jié)合楔绞，找到互作位點(diǎn)结闸。
- 超聲波剪切DNA鏈
- 加上 帶有抗體的磁珠用于 免疫沉淀靶蛋白【贫洌抗體很重要桦锄。
- 接觸蛋白交聯(lián)，純化DNA
- 送測(cè)序蔫耽，后續(xù)分析结耀。
-
實(shí)驗(yàn)設(shè)計(jì)的關(guān)鍵：
- 抗體質(zhì)量、
- 樣本量（10-50ng DNA匙铡，PCR擴(kuò)增越少图甜，偏差越少）
- 空白對(duì)照：（排除假陽性的存在，開放染色體區(qū)域慰枕、重復(fù)序列） input DNA對(duì)照具帮？？低斋？蜂厅、mock IP DNA（免疫共沉淀不含有抗體的DNA）、
三膊畴、CHIP 分析流程：
-
CHIP-SEQ 分析流程，分析分為4步
- 質(zhì)量控制，用的是Fastqc等
- 序列比對(duì)买猖，Bowtie2或BWA
- peak calling改橘， MACS
- peak注釋， ChIPseeker
-
CHIP數(shù)據(jù)分析所特有的步驟：peak calling：
- 染色體上信號(hào)波形的定義玉控；
- 建立背景矯正模型飞主；
- 建立搜索peaks的準(zhǔn)則，即建立判斷怎樣可以是一個(gè)peak（一般會(huì)用到背景矯正模型中的背景值）；
- 矯正模型碌识，過濾掉假陽性的peaks碾篡；
- 給找到的peaks排序，給出顯著度筏餐。
read只是跟隨著TF一起沉淀下來的DNA fragment的末端开泽，read的位置并不是真實(shí)的TF結(jié)合的位置。所以在peak-calling之前胖烛，延伸read是必須的眼姐。不同TF大小不一樣，對(duì)read延伸的長度也理應(yīng)不同佩番。我們知道众旗，測(cè)得的read最終其實(shí)會(huì)近似地平均分配到正負(fù)鏈上，這樣趟畏，對(duì)于一個(gè)TF結(jié)合熱點(diǎn)而言贡歧，read在附近正負(fù)鏈上會(huì)近似地形成“雙峰”。MACS會(huì)以某個(gè)window size掃描基因組赋秀，統(tǒng)計(jì)每個(gè)window里面read的富集程度利朵，然后抽取（比如1000個(gè)）合適的（read富集程度適中猎莲，過少绍弟，無法建立模型，過大著洼，可能反映的只是某種偏好性）window作樣本樟遣，建立“雙峰模型”。最后身笤，兩個(gè)峰之間的距離就被認(rèn)為是TF的長度D豹悬，每個(gè)read將延伸D/2的長度。
四液荸、具體代碼記錄
目標(biāo)：獲取感興趣的基因如ARF3在已發(fā)表的chip-seq數(shù)據(jù)AP3/PI/SEP3中是否存在相同的峰（位置瞻佛，高度）
1. 文章數(shù)據(jù)下載
- SEP3: Target Genes of the MADS Transcription
Factor SEPALLATA3: Integration of
Developmental and Hormonal Pathways
in the Arabidopsis Flower GSE:GSE14600- AP3-PI:Molecular basis for the specification of floral organs by
APETALA3 and PISTILLATA GSE:GSE38358- AG : Control of Reproductive Floral Organ Identity Specification in
Arabidopsis by the C Function Regulator AGAMOUS GSE :GSE45938- AP1: Orchestration of Floral Initiation
by APETALA1- AP1/SEP3 GSE46986
- 根據(jù) GSE號(hào) 以及 PRJNA 號(hào)下載
### 以AG數(shù)據(jù)為例
tail -n +2 AG-SraRunTable.txt |cut -f 7 |xargs -i nohup prefetch {} \&
- sra 數(shù)據(jù)轉(zhuǎn)換fastq.gz，注意也可以不用--gzip參數(shù)
### 測(cè)序數(shù)據(jù)皆為單端測(cè)序娇钱，故不需要--split-3
for i in `ls *.sra`
do
nohup fastq-dump.2 --gzip -O 1_fastq/ $i &
done
2. 數(shù)據(jù)的質(zhì)控
2.1 fastqc 處理
### 直接fastqc顯示測(cè)序質(zhì)量
ls *.fq.gz |parallel fastqc -o ./ --nogroup {} &
### MultiQC批量顯示測(cè)序的質(zhì)量
##在已經(jīng)用fastqc得到測(cè)序的html文件和zip文件夾
multiqc *fastqc.zip --pdf
2.2 fastp處理
###fastp快速質(zhì)控
fastp -i SRR502857.fastq.gz -o SRR502857_qc.fq.gz
### 出現(xiàn)adapter sequences的自動(dòng)檢測(cè)出現(xiàn)錯(cuò)誤的問題伤柄，加上-A參數(shù)
###根據(jù)fastqc報(bào)告查看接頭序列指定接頭的序列
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
### 結(jié)果發(fā)現(xiàn)測(cè)序質(zhì)量十分不錯(cuò)。
3. 序列的比對(duì)
### 構(gòu)建索引 基本命令就是bowtie2-build 基因組序列.fa 索引名字
bowtie2-build Arabidopsis_thaliana.TAIR10.dna.toplevel.fa ath_index_bt1
### 比對(duì)
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對(duì)于bowtie1 需要參數(shù)-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
- samtools faidx Arabidopsis_thaliana.TAIR10.dna.toplevel.fa 1:14331064-14331662
4. Samtools轉(zhuǎn)換
samtools view -bS SRR851694.sam | samtools sort -@ 4 - -T $i -o SRR851694.sorted.bam &
samtools index SRR851694.sorted.bam
5. BAM 轉(zhuǎn) bigwig文件 后續(xù)在IGV里查看各個(gè)基因的峰
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環(huán)境适刀。
### macs 較為嚴(yán)格的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參數(shù)改進(jìn)
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)
###轉(zhuǎn)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)
其他
- 從染色體指定區(qū)段獲取序列，并進(jìn)行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 轉(zhuǎn) bigwig文件
bamCoverage -b AG_IP.sorted.bam --normalizeUsing RPKM --binSize 30 --smoothLength 300 -p 10 --extendReads 200 -o AG_IP.bw
