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

可以將測序文件修改為自己樣品的名字
2）從NCBI的GEO數據庫中下載sra文件，需要轉換為.fastq文件。

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

從GEO數據庫中下載的.sra文件
#如果是單端測序文件就只需要fastq-dump命令將sra文件轉換為fastq測序文件。
$ fastq-dump SRR10485695.sra

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

#轉換為fastq文件后需要壓縮
$ gzip *.fastq
2.數據處理

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

具體用法可見：

FelixKrueger/TrimGalore
?github.com/FelixKrueger/TrimGalore

參數說明：-o 輸出到trim文件夾中；–cores 8線程數為8，默認是20 ；–fastqc對去完接頭的數據進行fastqc質量控制。具體參數可以使用$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序列比對 (比對軟件還有bwa和bowtie<和bowtie2不同>，這里使用的是bowtie2)。Bowtie 2是一個超快且內存高效的工具，用于將序列讀取對齊到參考基因組上。參數說明：-x是使用bowtie2索引，索引位置會根據不同物種發(fā)生改變，索引文件在/home/reference/下。-p為線程數(需根據服務器的總線程數考慮)；-U為需要比對的文件；|為通道符，將前面的值賦予后面的命令；samtools sort是對bam文件進行排序，一般情況按照序列在fasta文件中的順序和序列從左到右的位點排序；-@為設置線程數，同上-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去除重復。在制備文庫的過程中，由于PCR擴增過程中會存在一些偏差，也就是說有的序列會被過量擴增。這樣，在比對的時候，這些過量擴增出來的完全相同的序列就會比對到基因組的相同位置。而這些過量擴增的reads并不是基因組自身固有序列，不能作為變異檢測的證據，因此，要盡量去除這些由PCR擴增所形成的duplicates。主要的重復來源于PCR的擴增，還有部分的測序產生的重復。去重復的過程是給這些序列設置一個flag以標志它們，方便GATK的識別。這里定義的重復序列是這樣的：如果兩條reads具有相同的長度而且比對到了基因組的同一位置，那么就認為這樣的reads是由PCR擴增而來，就會被GATK標記。參數說明：-I為輸入需要去除重復的樣本。–ADD_PG_TAG_TO_READS為是否添加PG tag； --REMOVE_SEQUENCING_DUPLICATES為是否去除重復序列；–CREATE_INDEX為建造索引。-O為輸出設定；-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格式轉換。為了能夠比較不同的樣本，需要對先將基因組分成等寬分箱(bin)，統(tǒng)計每個分箱的read數，最后得到描述性統(tǒng)計值。

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

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

$ 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
參數說明：awk是一種處理文本文件的語言，是一個強大的文本分析工具。行匹配語句 awk ’ ’ 只能用單引號。首先將compareInput.bedgraph中第四列小于0的值轉換為0；sum+=$4表示sum=sum+$4，再輸出sum=totalreads；-v是賦值一個用戶定義變量。將第二步計算出的值賦予此命令中，將原來第四列的值比上totalreads再×1000000，最后輸出到第四列(改變第四列的值)。

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

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

7）將bedgraph轉換為bw文件。參數說明：根據dm6的染色體具體位置，將bedgraph轉為bw。Bedgraph文件中的染色體位置是相對位置，需要實際染色體的文件，轉換為絕對位置，即真實位置。

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

8）bamCoverage標準化處理。將bam文件轉為bw文件，同時進行標準化。參數說明：-p設置線程數；-b輸入bam文件；–skipNAs同上；–normalizeUsing標準化處理設置為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
參數說明：-t實驗組，IP的數據文件；-c對照組，input。腳本會將輸出的文件放在當前文件夾的peak文件夾中。

MACS2輸出文件解讀：

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

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

3.結果可視化

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

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

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

使用ChIPseeker進行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
參數說明：-a為peak文件；-b為注釋文件。輸出txt文件為注釋結果文件。(-a -b的表格盡量不要有tittle，不然容易出bug)

3）peak數據制作矩陣(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
參數說明：scale-regiuons mode簡單來說會將給定bed file范圍內的結合信號做一個統(tǒng)計（指的是一段長度），并將基因長度統(tǒng)一scale到設定regionBdoyLength的長度，加上統(tǒng)計基因上游和下游X bp的信號。reference-point mode則是給定一個bed file，以某個點為中心開始統(tǒng)計信號（TSS/TES/center）。但實際上我在嘗試的時候regionBdoyLength參數也還是可以用的，所以估計和scale-regions區(qū)別也不是太大，主要是作圖的一點區(qū)別。

4）作圖plotProfile & plotHeatmap

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

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

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。