1）自己實(shí)驗(yàn)獲得的原始測(cè)序文件（.fastq.gz），測(cè)序文件格式為.fastq格式，而.gz格式是壓縮格式，為了節(jié)省存儲(chǔ)空間。

可以將測(cè)序文件修改為自己樣品的名字
2）從NCBI的GEO數(shù)據(jù)庫中下載sra文件，需要轉(zhuǎn)換為.fastq文件。

$ wget -c ‘下載地址’ #wget下載數(shù)據(jù)，-c參數(shù)為斷點(diǎn)續(xù)傳

從GEO數(shù)據(jù)庫中下載的.sra文件
#如果是單端測(cè)序文件就只需要fastq-dump命令將sra文件轉(zhuǎn)換為fastq測(cè)序文件。
$ fastq-dump SRR10485695.sra

#如果是雙端測(cè)序文件轉(zhuǎn)換，需要加入–split-e參數(shù)，將.fastq測(cè)序文件轉(zhuǎn)換為*_1.fastq和*_2.fastq。
$ fastq-dump --split-e SRR10485695.sra

#轉(zhuǎn)換為fastq文件后需要壓縮
$ gzip *.fastq
2.數(shù)據(jù)處理

1）trim_galore質(zhì)量過濾，這一個(gè)非常流行的用于「去接頭序列」的軟件，用于處理高通量測(cè)序得到的原始數(shù)據(jù)。通常我們從測(cè)序公司拿到數(shù)據(jù)后，第一步就是評(píng)估數(shù)據(jù)的質(zhì)量以及對(duì)raw data去接頭處理。公司拿來的數(shù)據(jù)通常附帶了clean data以及去接頭的說明文件。

具體用法可見：

FelixKrueger/TrimGalore
?github.com/FelixKrueger/TrimGalore

參數(shù)說明：-o 輸出到trim文件夾中；–cores 8線程數(shù)為8，默認(rèn)是20 ；–fastqc對(duì)去完接頭的數(shù)據(jù)進(jìn)行fastqc質(zhì)量控制。具體參數(shù)可以使用$trim_galore -help查看。

fastqc具體用法見：

FastQC A Quality Control tool for High Throughput Sequence Data
?www.bioinformatics.babraham.ac.uk/projects/fastqc/

$ trim_galore .fastq.gz -gzip -o trim --cores 8 --fastqc
2）bowtie2序列比對(duì) (比對(duì)軟件還有bwa和bowtie<和bowtie2不同>，這里使用的是bowtie2)。Bowtie 2是一個(gè)超快且內(nèi)存高效的工具，用于將序列讀取對(duì)齊到參考基因組上。參數(shù)說明：-x是使用bowtie2索引，索引位置會(huì)根據(jù)不同物種發(fā)生改變，索引文件在/home/reference/下。-p為線程數(shù)(需根據(jù)服務(wù)器的總線程數(shù)考慮)；-U為需要比對(duì)的文件；|為通道符，將前面的值賦予后面的命令；samtools sort是對(duì)bam文件進(jìn)行排序，一般情況按照序列在fasta文件中的順序和序列從左到右的位點(diǎn)排序；-@為設(shè)置線程數(shù)，同上-p；-O輸出文件夾，-o輸出文件名。

$ bowtie2 -x /home/reference/Dmelanogaster/UCSC/dm6/Sequence/Bowtie2Index/dm6 -p 30 -U _trimmed.fq.gz | samtools sort -@ 30 -O bam -o bam/_sort.bam
3）gatk去除重復(fù)。在制備文庫的過程中，由于PCR擴(kuò)增過程中會(huì)存在一些偏差，也就是說有的序列會(huì)被過量擴(kuò)增。這樣，在比對(duì)的時(shí)候，這些過量擴(kuò)增出來的完全相同的序列就會(huì)比對(duì)到基因組的相同位置。而這些過量擴(kuò)增的reads并不是基因組自身固有序列，不能作為變異檢測(cè)的證據(jù)，因此，要盡量去除這些由PCR擴(kuò)增所形成的duplicates。主要的重復(fù)來源于PCR的擴(kuò)增，還有部分的測(cè)序產(chǎn)生的重復(fù)。去重復(fù)的過程是給這些序列設(shè)置一個(gè)flag以標(biāo)志它們，方便GATK的識(shí)別。這里定義的重復(fù)序列是這樣的：如果兩條reads具有相同的長(zhǎng)度而且比對(duì)到了基因組的同一位置，那么就認(rèn)為這樣的reads是由PCR擴(kuò)增而來，就會(huì)被GATK標(biāo)記。參數(shù)說明：-I為輸入需要去除重復(fù)的樣本。–ADD_PG_TAG_TO_READS為是否添加PG tag； --REMOVE_SEQUENCING_DUPLICATES為是否去除重復(fù)序列；–CREATE_INDEX為建造索引。-O為輸出設(shè)定；-M為創(chuàng)造矩陣文本。

$ gatk MarkDuplicates -I _sort.bam --ADD_PG_TAG_TO_READS false --REMOVE_SEQUENCING_DUPLICATES true --CREATE_INDEX true -O bam/.rmdup.bam -M bam/.rmdup.matrix.txt
4）bamCompare格式轉(zhuǎn)換。為了能夠比較不同的樣本，需要對(duì)先將基因組分成等寬分箱(bin)，統(tǒng)計(jì)每個(gè)分箱的read數(shù)，最后得到描述性統(tǒng)計(jì)值。

$ bamCompare -p 10 -b1 .rmdup.bam -b2 ._IN.rmdup.bam --skipNAs --scaleFactorsMethod readCount --operation subtract --outFileFormat bedgraph -o …/bw_data/.compareInput.bedgraph --extendReads 200
參數(shù)說明：-p為線程數(shù)；-b1為ip；-b2為input；–skipNA是跳過bin中沒有覆蓋的部分； --scaleFactorsMethod是選擇縮放樣本的方法，可選參數(shù)為：readCount/SES/ None； --operation默認(rèn)情況下輸出兩個(gè)示例的log2比率，此腳本選用subtract為做差值；–outFileFormat設(shè)置輸出文件格式，可選bigwig及bedgraph；-o為輸出設(shè)置； --extendReads此參數(shù)允許將讀取擴(kuò)展為片段大小。通常不建議對(duì)拼接讀取的數(shù)據(jù)(比如RNA-seq)使用此特性，因?yàn)樗鼤?huì)在跳過的區(qū)域上擴(kuò)展讀取。默認(rèn)參數(shù)為200。

5）compareinput to move0 to rpm.bedgraph：上一步做完差值后，可能會(huì)存在負(fù)值，所以這一步需要將其矯正為0，為之后的統(tǒng)計(jì)做準(zhǔn)備。

$ awk ‘{if($4<0){$4=0};print}’ .compareInput.bedgraph > .move0.bedgraph
$ totalreads=awk '{sum+=$4}END{print sum}' <sample>.move0.bedgraph |echo 'reads left after remove duplication: '$totalreads
$ awk -v totalreads=$totalreads ‘{$4=$4/totalreads*1000000;print}’ .move0.bedgraph > .rpm.bedgraph
參數(shù)說明：awk是一種處理文本文件的語言，是一個(gè)強(qiáng)大的文本分析工具。行匹配語句 awk ’ ’ 只能用單引號(hào)。首先將compareInput.bedgraph中第四列小于0的值轉(zhuǎn)換為0；sum+=$4表示sum=sum+$4，再輸出sum=totalreads；-v是賦值一個(gè)用戶定義變量。將第二步計(jì)算出的值賦予此命令中，將原來第四列的值比上totalreads再×1000000，最后輸出到第四列(改變第四列的值)。

6）sort排序：sort將文件的每一行作為一個(gè)單位，相互比較，比較原則是從首字符向后，依次按ASCII碼值進(jìn)行比較，最后將他們按升序輸出。參數(shù)說明：-k是指定列數(shù)；-k1,1 -k2,2n是先按照第一列排序，如果相同就按照第二列，n表示數(shù)字類型。n后面跟著r的，表示reverse。

$ sort -k1,1 -k2,2n .rpm.bedgraph > .sort.bedgraph

7）將bedgraph轉(zhuǎn)換為bw文件。參數(shù)說明：根據(jù)dm6的染色體具體位置，將bedgraph轉(zhuǎn)為bw。Bedgraph文件中的染色體位置是相對(duì)位置，需要實(shí)際染色體的文件，轉(zhuǎn)換為絕對(duì)位置，即真實(shí)位置。

$ bedGraphToBigWig .sort.bedgraph /public/software/ucsc-tools/dm6.chrom.sizes .rpm.bw

8）bamCoverage標(biāo)準(zhǔn)化處理。將bam文件轉(zhuǎn)為bw文件，同時(shí)進(jìn)行標(biāo)準(zhǔn)化。參數(shù)說明：-p設(shè)置線程數(shù)；-b輸入bam文件；–skipNAs同上；–normalizeUsing標(biāo)準(zhǔn)化處理設(shè)置為CPM，-o輸出。

$ bamCoverage -p $nThread?bbam/${sample}.rmdup.bam --skipNAs --normalizeUsing CPM -o bw/${sample}.CPM.bw
9）call peak：支持bam/sam/bed等文件格式。

$ macs2 callpeak -t ._sort.bam -c ._sort.bam
參數(shù)說明：-t實(shí)驗(yàn)組，IP的數(shù)據(jù)文件；-c對(duì)照組，input。腳本會(huì)將輸出的文件放在當(dāng)前文件夾的peak文件夾中。

MACS2輸出文件解讀：

① NAME_peaks.xls——包含peak信息的tab分割的文件，前幾行會(huì)顯示callpeak時(shí)的命令。輸出信息包含：染色體號(hào)、peak起始位點(diǎn)、peak結(jié)束位點(diǎn)、peak區(qū)域長(zhǎng)度、peak的峰值位點(diǎn)(summit position)、peak 峰值的高度(pileup height at peak summit, -log10(pvalue) for the peak summit)、peak的富集倍數(shù)(相對(duì)于random Poisson distribution with local lambda)。<XLS里的坐標(biāo)和bed格式的坐標(biāo)還不一樣，起始坐標(biāo)需要減1才與narrowPeak的起始坐標(biāo)一樣>

② NAME_peaks.narrowPeak——1：染色體號(hào)；2：peak起始位點(diǎn)；3：結(jié)束位點(diǎn)；4：peak name；5：int(-10*log10qvalue)；6：正負(fù)鏈；7：fold change；8：-log10pvalue；9：-log10qvalue；10：relative summit position to peak start。

3.結(jié)果可視化

1）利用bedtools intersect取overlap （可選步驟！！ 這個(gè)是求兩個(gè)peak的交集）

$ bedtools intersect -a .bed -b .bed -wa -wb > overlap_.txt
$ bedtools window -w n #利用設(shè)置一定的大小的window進(jìn)行overlap，-w為設(shè)置window大小為n
參數(shù)說明：bedtools intersect為兩個(gè)取交集；bedtools window是取一定大小的窗口，進(jìn)行overlap，即可以存在最多窗口大小的差距再取overlap。-a -b是指定兩個(gè)文件；-wa -wb是將兩個(gè)文件中overlap的列出來；>輸出到指定文件中。

2）利用bedtools intersect做peak注釋。也可以用R進(jìn)行peak注釋，可以參考：

使用ChIPseeker進(jìn)行peak注釋_廬州月光的博客-CSDN博客
?blog.csdn.net/weixin_43569478/article/details/108079450

$ bedtools intersect -a stage6_peak.bed -b UCSC.D.melanogaster.dm6.repeats.gtf -wa -wb > overlap_repeat_anno.txt
參數(shù)說明：-a為peak文件；-b為注釋文件。輸出txt文件為注釋結(jié)果文件。(-a -b的表格盡量不要有tittle，不然容易出bug)

3）peak數(shù)據(jù)制作矩陣(scale-regions/reference-point)

$ computeMatrix scale-regions -S .bw input..bw -R .bed -a 1000 -b 1000 -o ip_input.computeMatrix
$ computeMatrix reference-point --referencePoint center -b .bed -S .bw .bw --skipZeros -o .gz --outFileNameMatrix
參數(shù)說明：scale-regiuons mode簡(jiǎn)單來說會(huì)將給定bed file范圍內(nèi)的結(jié)合信號(hào)做一個(gè)統(tǒng)計(jì)（指的是一段長(zhǎng)度），并將基因長(zhǎng)度統(tǒng)一scale到設(shè)定regionBdoyLength的長(zhǎng)度，加上統(tǒng)計(jì)基因上游和下游X bp的信號(hào)。reference-point mode則是給定一個(gè)bed file，以某個(gè)點(diǎn)為中心開始統(tǒng)計(jì)信號(hào)（TSS/TES/center）。但實(shí)際上我在嘗試的時(shí)候regionBdoyLength參數(shù)也還是可以用的，所以估計(jì)和scale-regions區(qū)別也不是太大，主要是作圖的一點(diǎn)區(qū)別。

4）作圖plotProfile & plotHeatmap

$ plotHeatmap -m .gz -o heatmap_.pdf -zMin 0 --zMax 8 --colorMap coolwarm --missingDataColor 1
參數(shù)說明：-m輸入矩陣；-o輸出pdf文件；-zMin 0 --zMax 8是熱圖的最大值和最小值；–colorMap是選擇熱圖顏色類型。做出的pdf熱圖可以直接在Ai里編輯。

ps. 溫馨提示，具體操作都需要自己實(shí)踐，實(shí)踐的路上總會(huì)出現(xiàn)bug。以上內(nèi)容如果有什么建議和意見，可以私信，歡迎討論。

總結(jié)

以上是生活随笔為你收集整理的ChIP-seq分析流程(基于linux系统)的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

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