以Targeting super-enhancer-associated oncogenes
in oesophageal squamous cell carcinoma為例，實現(xiàn)其Chip-seq分析。在Chip-seq過程中，我們需要尋找食管癌致癌基因的超級增強子區(qū)域以預(yù)測其在臨床中有可能的應(yīng)用。

步驟概覽

1. 下載Chip-seq原始數(shù)據(jù)
2. fastqc質(zhì)量檢測
3. 下載人類參考基因組并建立index
4. 使用bowtie比對
5. 使用MACS獲得Chip-seq富集區(qū)
6. 使用IGV工具可視化
7. 使用ROSE篩選super-enhancer

下載Chip-seq原始數(shù)據(jù)

文獻(xiàn)中提到，其原始數(shù)據(jù)上傳在NCBI的GEO Dataset數(shù)據(jù)庫中，編號GSE76861，在NCBI數(shù)據(jù)庫中搜索到結(jié)果

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE76861

在Samples欄目下，可以看到其上傳的所有數(shù)據(jù)（點擊More…來展開），其中

GSM2039110  TE7_H3K27Ac
GSM2039111  TE7_Input
GSM2039112  KYSE510_H3K27Ac
GSM2039113  KYSE510_Input

為Chip-seq數(shù)據(jù)，我們使用aspera工具來下載，首先在ebi中找到相應(yīng)的序列，獲得4個數(shù)據(jù)的下載鏈接，我們?nèi)コ懊娴挠蛎幚淼梦谋救缦?/p>

/vol1/fastq/SRR310/001/SRR3101251/SRR3101251.fastq.gz
/vol1/fastq/SRR310/002/SRR3101252/SRR3101252.fastq.gz
/vol1/fastq/SRR310/003/SRR3101253/SRR3101253.fastq.gz
/vol1/fastq/SRR310/004/SRR3101254/SRR3101254.fastq.gz

以方便批量下載。
然后運行aspera來下載

> ascp -QT -k1 -l 100M -i ~/asperaweb_id_dsa.openssh --mode recv --host fasp.sra.ebi.ac.uk --user era-fasp --file-list file-list  .

這里的-QT表示開啟斷點續(xù)傳，-l表示最大限速帶寬，-i參數(shù)輸入密匙文件路徑，--file-list參數(shù)即代表我們創(chuàng)建的下載鏈接文本。
下載完成后，進(jìn)行解壓

> gunzip SRR3101251.fastq.gz SRR3101252.fastq.gz SRR3101253.fastq.gz SRR3101254.fastq.gz

fastqc質(zhì)量檢測

運行命令

> fastqc -o ../fastqcresult/ -f fastq SRR3101251.fastq SRR3101252.fastq SRR3101253.fastq SRR3101254.fastq

來對序列質(zhì)量進(jìn)行檢測，-o參數(shù)代表結(jié)果輸出位置，-f參數(shù)表示輸入文件的格式，我們是fastq。

下載人類參考基因組并建立index

原則上來講，我們需要下載人類參考基因組，然后使用bowtie-build命令來建立index，但是bowtie的官網(wǎng)已經(jīng)給出了index文件并且提供下載，直接在bowtie官網(wǎng)下載

> wget ftp://ftp.ccb.jhu.edu/pub/data/bowtie_indexes/hg19.ebwt.zip
> gunzip hg19.ebwt.zip

這樣，我們就能得到人類的index了。值得注意的是，bowtie的官網(wǎng)提供了bowtie2和bowtie的兩種index，同時，又分為NCBI提供的基因組數(shù)據(jù)和UCSC提供的基因組數(shù)據(jù)，我們使用和文獻(xiàn)一致的UCSC提供基因組數(shù)據(jù)。

使用bowtie比對

bowtie是整個分析的開始，bowtie會將所有的序列和基因組比對，以給出其在基因組所在的位置，整個過程會比較長，我們需要將其放在后臺中運行，以免其中斷

> nohup bowtie index/hg19 -q data/SRR3101251.fastq -v 2 -m 1 -3 1 -S 2>bowtieresult/SRR3101251.out>bowtieresult/SRR3101251.sam &
> nohup bowtie index/hg19 -q data/SRR3101252.fastq -v 2 -m 1 -3 1 -S 2>bowtieresult/SRR3101252.out>bowtieresult/SRR3101252.sam &
> nohup bowtie index/hg19 -q data/SRR3101253.fastq -v 2 -m 1 -3 1 -S 2>bowtieresult/SRR3101253.out>bowtieresult/SRR3101253.sam &
> nohup bowtie index/hg19 -q data/SRR3101254.fastq -v 2 -m 1 -3 1 -S 2>bowtieresult/SRR3101254.out>bowtieresult/SRR3101254.sam &

其中，nohup能夠使命令成為無主進(jìn)程，脫離ssh的進(jìn)程，使得整個程序的運行能夠在斷開ssh連接之后仍然在進(jìn)行。同時將其標(biāo)準(zhǔn)輸出放入sam文件，錯誤輸出放入out文件。

使用MACS獲得Chip-seq富集區(qū)

由于bowtie僅給出了每個read的位置，為了統(tǒng)計其數(shù)量，我們需要使用，MACS進(jìn)行富集

> nohup macs14 -t SRR3101251.sam  -c SRR3101252.sam --format SAM --name "TE7" --keep-dup 1 --wig --single-profile --space=50 --diag &
> nohup macs14 -t SRR3101253.sam  -c SRR3101254.sam --format SAM --name "KYSE510" --keep-dup 1 --wig --single-profile --space=50 --diag &

與上面相同，MACS同樣需要運行比較長的時間，所以仍然需要其在后臺運行。其次，這里的參數(shù)-t表示實驗組，-c表示對照組，--format表示輸入文件的格式，--name輸出文件的附加前綴，--keep-dup對于重復(fù)序列的處理方式，1效果最好，--wig表示輸出wig文件，--single-profile表示輸出單文件，--space是文獻(xiàn)中的要求。

使用IGV工具可視化

首先在IGV官網(wǎng)下載并安裝軟件，采用其windows版本，為了能夠可視化.wig文件，我們需要將其轉(zhuǎn)換為.bw文件，首先需要使用fetchChromSizes

> wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64.v287/fetchChromSizes
> chmod 777 fetchChromSizes

使用chmod命令來獲得其執(zhí)行權(quán)限，然后我們獲取基因組長度信息

> fetchChromSizes hg19 >hg19.chrom.sizes

接下來下載wigToBigWig來轉(zhuǎn)換格式

> wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/wigToBigWig
> chmod 777 wigToBigWig

使用wigToBigWig轉(zhuǎn)換格式

> wigToBigWig TE7_control_afterfiting_all.wig   ../index/hg19.chrom.sizes TE7_control_afterfiting_all.bw
> wigToBigWig TE7_treat_afterfiting_all.wig   ../index/hg19.chrom.sizes TE7_treat_afterfiting_all.bw
> wigToBigWig  KYSE510_control_afterfiting_all.wig ../index/hg19.chrom.sizes KYSE510_control_afterfiting_all.bw
> wigToBigWig  KYSE510_treat_afterfiting_all.wig  ../index/hg19.chrom.sizes KYSE510_treat_afterfiting_all.bw

然后將得到的.bw文件導(dǎo)入IGV軟件即可看到效果圖

KYSE510效果圖

TE7效果圖

可以看到圖中用方框標(biāo)示的區(qū)域，在實驗組中，峰值都超過了域值，而對照組依然在域值內(nèi)，因此，我們認(rèn)為這些基因就是peeks。

使用ROSE篩選super-enhancer

為了找到哪些是super-enhancer，我們需要使用ROSE進(jìn)行篩選，首先需要將.sam文件轉(zhuǎn)換為.bam,我們使用samtools進(jìn)行轉(zhuǎn)換

> samtools view -b -u bt/SRR3101251.sam >SRR3101251.bam
> samtools view -b -u bt/SRR3101252.sam >SRR3101252.bam
> samtools view -b -u bt/SRR3101253.sam >SRR3101253.bam
> samtools view -b -u bt/SRR3101254.sam >SRR3101254.bam

然后對其進(jìn)行排序

> samtools sort SRR3101251.bam SRR3101251.sorted
> samtools sort SRR3101252.bam SRR3101252.sorted
> samtools sort SRR3101253.bam SRR3101253.sorted
> samtools sort SRR3101254.bam SRR3101254.sorted

再建立索引

> samtools index SRR3101251.sorted.bam
> samtools index SRR3101252.sorted.bam
> samtools index SRR3101253.sorted.bam
> samtools index SRR3101254.sorted.bam

安裝ROSE，我們直接從其托管在Bitbucket倉庫中克隆Python腳本

git clone https://bitbucket.org/young_computation/rose.git

最后進(jìn)行篩選

> nohup python ROSE_main.py -g HG19 -i /macsresult/TE7_peaks.bed -r /bowtieresult/SRR3101251.sorted.bam -c /bowtieresult/SRR3101252.sorted.bam -o /roseresult -s 12500 -t 2000 &
> nohup python ROSE_main.py -g HG19 -i /macsresult/KYSE510_peaks.bed -r /bowtieresult/SRR3101253.sorted.bam -c /bowtieresult/SRR3101254.sorted.bam -o /roseresult -s 12500 -t 2000 &

其中-g為基因組名，-i為輸入的文件，-r為實驗組數(shù)據(jù)，-c為對照組數(shù)據(jù)，-o為輸出文件夾，-s為相鄰的峰合并，-t除去起始子的距離。
最后我們得到結(jié)果如下圖

TE7組結(jié)果

KYSE510組結(jié)果

從圖中，我們也可以清晰地分辨super-enhancer的分界同時，所有的super-enhancer也列舉在*_peaks_SuperEnhancers.table.txt文件中，至此我們已經(jīng)完成了所有的步驟，獲得了結(jié)果。

?