文章目錄

• CHIP-Seq數(shù)據(jù)分析流程
• • 相關(guān)軟件安裝
  • 數(shù)據(jù)下載
  • 將SRA轉(zhuǎn)化為fastq文件
  • 數(shù)據(jù)質(zhì)控，過濾
  • • 數(shù)據(jù)質(zhì)控
    • 然后進(jìn)行質(zhì)控，過濾低質(zhì)量reads，去接頭
  • 比對(duì)
  • • 下載小鼠參考基因組的索引和注釋文件
    • 比對(duì)
  • sam文件轉(zhuǎn)bam
  • 為bam文件建立索引
  • 載入IGV查看
  • 用MACS call peak
  • • 安裝MACS
  • 結(jié)果注釋與可視化
  • • 使用deeptools進(jìn)行可視化
    • peaks注釋

CHIP-Seq數(shù)據(jù)分析流程

相關(guān)軟件安裝

所需要的軟件有sratoolkit, fastqc,bowtie2,samtools,macs2,htseq-count,bedtools和前面RNA-seq數(shù)據(jù)分析很多軟件相同，所以基本都安裝過了。后面遇到?jīng)]有安裝過的在進(jìn)行安裝

###igv axel -n 20 https://data.broadinstitute.org/igv/projects/downloads/2.8/IGV_Linux_2.8.12_WithJava.zip unzip IGV_Linux_2.8.12_WithJava.zip cd IGV_Linux_2.8.12/ ./igv.sh

數(shù)據(jù)下載

還是利用sratoolkit的prefetch下載，首先構(gòu)建SRR_Acc_List.txt

for ((i=204;i<=209;i++)) do echo SRR620$i >> SRR_Acc_List.txt done ###文件內(nèi)容如下 SRR620204 SRR620205 SRR620206 SRR620207 SRR620208 SRR620209 ##接下來(lái)進(jìn)行下載 wkd=/home/meiling/baiduyundisk/Chip-seq #設(shè)置工作目錄 cat SRR_Acc_List.txt | while read id do nohup prefetch ${id} & done

將SRA轉(zhuǎn)化為fastq文件

wkd=/home/meiling/baiduyundisk/Chip-seq #設(shè)置工作目錄for i in $wkd/rawdata/*sra doecho $ifastq-dump --split-3 --skip-technical --clip --gzip $i ## 批量轉(zhuǎn)換 done

數(shù)據(jù)質(zhì)控，過濾

數(shù)據(jù)質(zhì)控

fastqc生成質(zhì)控報(bào)告，multiqc將各個(gè)樣本的質(zhì)控報(bào)告整合為一個(gè)

wkd=/home/meiling/baiduyundisk/Chip-seq #設(shè)置工作目錄ls $wkd/rawdata/*gz | xargs fastqc -t 5 multiqc ./ ##得到以下結(jié)果 ls -lh multiqc_data/ total 7.1M -rw-r--r-- 1 meiling meiling 312K Nov 10 14:10 multiqc_data.json -rw-r--r-- 1 meiling meiling 1.5K Nov 10 14:10 multiqc_fastqc.txt -rw-r--r-- 1 meiling meiling 621 Nov 10 14:10 multiqc_general_stats.txt -rw-r--r-- 1 meiling meiling 24K Nov 10 14:10 multiqc.log -rw-r--r-- 1 meiling meiling 1.2M Nov 10 14:10 multiqc_report.html -rw-r--r-- 1 meiling meiling 604 Nov 10 14:10 multiqc_sources.txt -rw-r--r-- 1 meiling meiling 576K Nov 10 14:08 SRR620204_fastqc.html -rw-r--r-- 1 meiling meiling 431K Nov 10 14:08 SRR620204_fastqc.zip

每個(gè)id_fastqc.html都是一個(gè)質(zhì)量報(bào)告，multiqc_report.html是所有樣本的整合報(bào)告

然后進(jìn)行質(zhì)控，過濾低質(zhì)量reads，去接頭

wkd=/home/meiling/baiduyundisk/Chip-seq #設(shè)置工作目錄 mkdir $wkd/cleandata cd $wkd/cleandata ls $wkd/rawdata/*.fastq.gz >config ##開始質(zhì)控 wkd=/home/meiling/baiduyundisk/Chip-seq/cleandata #設(shè)置工作目錄 cd $wkd cat config |while read id doarr=(${id})fq1=${arr[0]} # fq2=${arr[1]} trim_galore -q 25 --phred33 --length 36 --stringency 3 -o $wkd $fq1 done

比對(duì)

下載小鼠參考基因組的索引和注釋文件

還是選擇在ensemble上下載小鼠參考基因組

選擇primary進(jìn)行下載

axel -n 20 ftp://ftp.ensembl.org/pub/release-101/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna.primary_assembly.fa.gz axel -n 20 ftp://ftp.ensembl.org/pub/release-101/gtf/mus_musculus/Mus_musculus.GRCm38.101.chr_patch_hapl_scaff.gtf.gz ##解壓 gzip -d Mus_musculus.GRCm38.dna.primary_assembly.fa.gz ls -lh total 2.7G -rw-r--r-- 1 meiling meiling 33M Nov 9 21:32 Mus_musculus.GRCm38.101.chr_patch_hapl_scaff.gtf.gz -rw-r--r-- 1 meiling meiling 2.6G Nov 9 21:31 Mus_musculus.GRCm38.dna.primary_assembly.fa #構(gòu)建索引 nohup bowtie2-build ../Mus_musculus.GRCm38.dna.primary_assembly.fa mouse &

比對(duì)

這里選擇bowtie,或者bwa
將得到的fastq文件用bowtie2比對(duì)小鼠參考基因組上

wkd=/home/meiling/baiduyundisk/Chip-seq #設(shè)置工作目錄 cd $wkd/cleandata ls *gz|cut -d"_" -f 1 |sort -u |while read id;do ls -lh ${id}_trimmed.fq.gz bowtie2 -q -p 10 -x $wkd/reference/mouse/bowtie_index/mouse -U ${id}_trimmed.fq.gz -S $wkd/align/${id}.sam done

比對(duì)結(jié)果輸出

bash align.sh >align.log 2>&1 #打開align.log -rw-r--r-- 1 meiling meiling 517M Nov 10 14:15 SRR620204_trimmed.fq.gz 12553187 reads; of these:12553187 (100.00%) were unpaired; of these:3148914 (25.08%) aligned 0 times7170824 (57.12%) aligned exactly 1 time2233449 (17.79%) aligned >1 times 74.92% overall alignment rate

sam文件轉(zhuǎn)bam

wkd=/home/meiling/baiduyundisk/Chip-seq #設(shè)置工作目錄 cd $wkd/align ls *.sam|cut -d"." -f 1 |while read id ;dosamtools view -@ 8 -S $id.sam -1b -o $id.bamsamtools sort -@ 8 -l 5 -o $id.sort.bam $id.bam done

為bam文件建立索引

ls *.sort.bam |xargs -i samtools index {}

載入IGV查看

從官網(wǎng)下載IGV, 解壓即可使用，linux下 igv.sh打開IGV界面
首先載入?yún)⒖蓟蚪M，可以載入自己下載好的參考基因組，也可選擇IGV中含有的參考基因組，ref_Genome 必須是fasta格式。
載入比對(duì)的文件，比對(duì)的文件必須先經(jīng)過sort 和 index, 才可加載。
比對(duì)可視化結(jié)果：其他組有的峰在對(duì)照組沒有，即為peak.

用MACS call peak

安裝MACS

#下載安裝 axel -n 20 https://files.pythonhosted.org/packages/e2/61/85d30ecdd34525113e28cb0c5a9f393f93578165f8d848be5925c0208dfb/MACS2-2.2.7.1.tar.gz tar -xvf MACS2-2.2.7.1.tar.gz cd MACS2-2.2.7.1/ ##安裝 python setup.py install ##驗(yàn)證 macs2 -h ##出現(xiàn)以下結(jié)果 usage: macs2 [-h] [--version]{callpeak,bdgpeakcall,bdgbroadcall,bdgcmp,bdgopt,cmbreps,bdgdiff,filterdup,predictd,pileup,randsample,refinepeak}...macs2 -- Model-based Analysis for ChIP-Sequencingpositional arguments:{callpeak,bdgpeakcall,bdgbroadcall,bdgcmp,bdgopt,cmbreps,bdgdiff,filterdup,predictd,pileup,randsample,refinepeak} ##call peak wkd=/home/meiling/baiduyundisk/Chip-seq #設(shè)置工作目錄 cd $wkd/align ls *.sort.bam|cut -d"." -f 1 |sort -u |while read id do macs2 callpeak -c SRR620208.bam -t ${id}.bam -q 0.05 -f BAM -g mm -n ${id} 2> ${id}.log & done # （在當(dāng)前目錄下）統(tǒng)計(jì) *bed 的行數(shù)（peak數(shù)） wc -l *bed 7765 SRR620204_summits.bed 1963 SRR620205_summits.bed 6775 SRR620206_summits.bed 0 SRR620207_summits.bed 0 SRR620208_summits.bed##對(duì)照組 0 SRR620209_summits.bed##對(duì)照組 16503 total #RYBP SRR620207 無(wú)數(shù)據(jù)，估計(jì)是數(shù)據(jù)上傳錯(cuò)誤，可以下載作者上傳的peak數(shù)據(jù)。 ● 下載RYBP的peak數(shù)據(jù) wget ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE42nnn/GSE42466/suppl/GSE42466_RYBP_peaks_5.txt.gz gzip -d GSE42466_RYBP_peaks_5.txt.gz mv GSE42466_RYBP_peaks_5.txt SRR620207_summits.bed

callpeak會(huì)得到如下結(jié)果文件：

NAME_summits.bed：Browser Extensible Data，記錄每個(gè)peak的peak summits，話句話說(shuō)就是記錄極值點(diǎn)的位置。MACS建議用該文件尋找結(jié)合位點(diǎn)的motif。能夠直接載入U(xiǎn)CSC browser，用其他軟件分析時(shí)需要去掉第一行。

NAME_peaks.xls：以表格形式存放peak信息，雖然后綴是xls，但其實(shí)能用文本編輯器打開，和bed格式類似，但是以1為基，而bed文件是以0為基.也就是說(shuō)xls的坐標(biāo)都要減一才是bed文件的坐標(biāo)。

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基本一致，適合用于導(dǎo)入R進(jìn)行分析。

NAME_model.r：能通過$ Rscript NAME_model.r作圖，得到是基于你提供數(shù)據(jù)的peak模型。

.bdg：能夠用 UCSC genome browser 轉(zhuǎn)換成更小的 bigWig 文件。

結(jié)果注釋與可視化

結(jié)果的注釋用的是Y叔 的 Chipseeker包。

ChIPseeker的功能分為三類: ● 注釋：提取peak附近最近的基因， 注釋peak所在區(qū)域 ● 比較：估計(jì)ChIP peak數(shù)據(jù)集中重疊部分的顯著性；整合GEO數(shù)據(jù)集，以便于將當(dāng)前結(jié)果和已知結(jié)果比較 ● 可視化： peak的覆蓋情況；TSS區(qū)域結(jié)合的peak的平均表達(dá)譜和熱圖；基因組注釋；TSS距離；peak和基因的重疊。
后續(xù)在R語(yǔ)言中實(shí)現(xiàn)。

使用deeptools進(jìn)行可視化

deeptools提供bamCoverage和bamCompare進(jìn)行格式轉(zhuǎn)換，為了能夠比較不同的樣本，需要對(duì)先將基因組分成等寬分箱(bin)，統(tǒng)計(jì)每個(gè)分箱的read數(shù)，最后得到描述性統(tǒng)計(jì)值。對(duì)于兩個(gè)樣本，描述性統(tǒng)計(jì)值可以是兩個(gè)樣本的比率，或是比率的log2值，或者是差值。如果是單個(gè)樣本，可以用SES方法進(jìn)行標(biāo)準(zhǔn)化。

##deeptools安裝 git clone https://github.com/fidelram/deepTools cd deepTools python setup.py install

bamCoverage的基本用法

peaks注釋

總結(jié)

以上是生活随笔為你收集整理的CHIP-Seq数据分析流程的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。