在實(shí)踐之前苟呐，我先查找了一下有關(guān)組蛋白修飾的知識點(diǎn)：
（參考文章：http://www.reibang.com/p/8aca72809c5c）

組蛋白修飾能預(yù)測染色質(zhì)的類型（異染色質(zhì)或常染色質(zhì)）怎静、區(qū)分基因組功能元件（啟動子焊虏、增強(qiáng)子、基因主體）以及檢測決定這些元件處于活性狀態(tài)或是抑制狀態(tài)。例如H3K4me2和H3K4me3修飾大多數(shù)富集在轉(zhuǎn)錄起始位點(diǎn)附近的啟動子上激活基因表達(dá)隧饼，而H3K27me2和H3K27me3與基因抑制相關(guān)进统。因此可通過CHIP-seq分析組蛋白修飾的分布尋找基因的啟動子區(qū)和增強(qiáng)子區(qū)域及其是激活或抑制基因表達(dá)助币。

H3K4me1可作為增強(qiáng)子的標(biāo)志，H3K4me3作為啟動子標(biāo)志螟碎。研究表明眉菱，H3K4me1和H3K4me3與基因激活相關(guān)，H3K4me3主要富集在轉(zhuǎn)錄起始位點(diǎn)附近的啟動子區(qū)域掉分，而大多數(shù)H3K4me1修飾富集在增強(qiáng)子區(qū)域俭缓；H3K27ac與基因激活相關(guān)克伊，主要富集在增強(qiáng)子和啟動子區(qū)域，當(dāng)增強(qiáng)子區(qū)只有H3K4me1修飾富集時华坦，該增強(qiáng)子處于平衡狀態(tài)愿吹，而當(dāng)增強(qiáng)子區(qū)域同時富集H3K4me1和H3K27ac修飾時，該增強(qiáng)子就處于激活狀態(tài)促進(jìn)基因表達(dá)惜姐；H3K27的甲基化是可逆的過程犁跪，H3K27me1顯示出對轉(zhuǎn)錄具有正向影響,啟動子區(qū)域的H3K27me3甲基化修飾時抑制基因的轉(zhuǎn)錄,而H3K27me2廣泛分布并且在沉默非細(xì)胞類型特異性增強(qiáng)子中起作用。

下表為常見的組蛋白修飾的主要分布及功能：

image

異染色質(zhì)是染色質(zhì)的濃縮歹袁，轉(zhuǎn)錄無活性狀態(tài)坷衍，H3K9甲基化是異染色質(zhì)的標(biāo)志。H3K27me1和H3K9me3存在于著絲粒異染色質(zhì)區(qū)域条舔，而H3K27me3和H3K9me2共同存在于抑制的常染色質(zhì)區(qū)域中枫耳。H3K9ac也與H3K14ac和H3K4me3高度共存共同作為活性基因啟動子的標(biāo)志。

本次實(shí)踐文章：
Targeting super-enhancer-associated oncogenes in oesophageal squamous cell carcinoma
用不同濃度的OSCC抗癌物THZ1處理兩種OSCC細(xì)胞系：KYSE510和TE7逞刷。
尋找藥物處理后的超級增強(qiáng)子嘉涌。

1.下載SRA原始數(shù)據(jù)：

#!/bin/bash
for ((i=251;i<=254;i++))
do wget https://sra-download.ncbi.nlm.nih.gov/traces/sra28/SRR/003028/SRR3101$i
done

2.提取fastq文件：

#!/bin/bash
for i in SRR310125*
do
        echo $i
        fastq-dump --gzip --split-files $i
done

每運(yùn)行一個文件，都會在家目錄下的ncbi文件里生成一個臨時文件夸浅，非常占內(nèi)存仑最。由于我的存儲空間沒有了，所以我就直接下載了fastq文件帆喇。下載fastq文件的網(wǎng)址是：https://www.ebi.ac.uk/ena/data/view/SRX1530372警医。再復(fù)制下fastq.gz文件的地址。
網(wǎng)上有的文章提到：盡量不要用wget或curl去下載sra文件坯钦，某些情況下會導(dǎo)致下載的文件不完整预皇。所以這次我也嘗試了用ascp下載的方法。
這里說一下用ascp下載的方法婉刀。

安裝ASCP:

#下載
$ wget http://download.asperasoft.com/download/sw/connect/3.7.4/aspera-connect-3.7.4.147727-linux-64.tar.gz
$ tar zxvf aspera-connect-3.7.4.147727-linux-64.tar.gz
#安裝
$ bash aspera-connect-3.7.4.147727-linux-64.sh
#要去home目錄里查看是否有.aspera文件夾
$ cd /home
$ ls -a 
#如果看到.aspera文件夾，代表安裝成功
#永久添加環(huán)境變量
$ echo 'export PATH=~/.aspera/connect/bin:$PATH' >> ~/.bashrc
$ source ~/.bashrc
#查看是否可以調(diào)用
$ ascp --help

ascp的用法：ascp [參數(shù)] 目標(biāo)文件 目標(biāo)地址突颊，在線文檔
ascp命令的常用參數(shù):
-v verbose mode 嘮叨模式，能讓你實(shí)時知道程序在干啥爬橡，方便查錯。
-T 取消加密棒动，否則有時候數(shù)據(jù)下載不了。
-i 提供私鑰文件的地址柜裸。
-l 設(shè)置最大傳輸速度缕陕，一般200m到500m，如果不設(shè)置疙挺，反而速度會比較低，可能有個較低的默認(rèn)值。
-k 斷點(diǎn)續(xù)傳锦爵，一般設(shè)置為值1
-Q 一般加上它
-P 提供SSH port奥裸，一般是33001。

在EBI網(wǎng)站copy下載地址一般是ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR310/004/SRR3101254/SRR3101254.fastq.gz這樣的格式樟氢。需要把ftp://ftp.sra.ebi.ac.uk換成era-fasp@fasp.sra.ebi.ac.uk:侠鳄。NOTE: 這篇文章的ChIP-seq是單端測序。
我的代碼：

$ ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR310/001/SRR3101251/SRR3101251.fastq.gz .
$ ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR310/002/SRR3101252/SRR3101252.fastq.gz .
$ ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR310/003/SRR3101253/SRR3101253.fastq.gz .
$ ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR310/004/SRR3101254/SRR3101254.fastq.gz .

另外要注意的是：era-fasp@fasp.sra.ebi.ac.uk后面是:號碴开，不是路徑
最后gz結(jié)尾有一個空格潦牛，然后是"."挡育。如果沒有空格和點(diǎn)，則會報錯橡淆。無法下載蒿叠。

3.使用fastqc質(zhì)量檢查
4.bowtie2比對

下面這些代碼是按照上一次實(shí)踐的方法寫的明垢，直接生成bam文件市咽。因?yàn)檫@篇文獻(xiàn)的測序結(jié)果質(zhì)量很好，所以沒有去掉3'端的幾個堿基
#!/bin/bash
bowtie2 -p 6 -x /media/yanfang/FYWD/CHIPSEQ/fqfile/bowtie2hg19/hg19 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR3101251.fastq.gz | samtools sort -O bam -o /media/yanfang/FYWD/CHIPSEQ/sambam/TE7_H3K27Ac.bam
bowtie2 -p 6 -x /media/yanfang/FYWD/CHIPSEQ/fqfile/bowtie2hg19/hg19 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR3101252.fastq.gz | samtools sort -O bam -o /media/yanfang/FYWD/CHIPSEQ/sambam/TE7_input.bam
bowtie2 -p 6 -x /media/yanfang/FYWD/CHIPSEQ/fqfile/bowtie2hg19/hg19 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR3101253.fastq.gz | samtools sort -O bam -o /media/yanfang/FYWD/CHIPSEQ/sambam/KYSE510_H3K27Ac.bam
bowtie2 -p 6 -x /media/yanfang/FYWD/CHIPSEQ/fqfile/bowtie2hg19/hg19 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR3101254.fastq.gz | samtools sort -O bam -o /media/yanfang/FYWD/CHIPSEQ/sambam/KYSE510_input.bam

下面這些代碼是bowtie2比對后生成sam文件的方法溯革。有些文章里在進(jìn)行比對的時候是用的bowtie進(jìn)行比對的，并且設(shè)置了只保留比對一次的reads致稀。我這里用的是bowtie2抖单，所以沒有設(shè)置這個參數(shù)。
#!/bin/bash
bowtie2 -p 6 -x /media/yanfang/FYWD/CHIPSEQ/fqfile/bowtie2hg19/hg19 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR3101251.fastq.gz -S TE7_H3K27Ac.sam
bowtie2 -p 6 -x /media/yanfang/FYWD/CHIPSEQ/fqfile/bowtie2hg19/hg19 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR3101252.fastq.gz -S TE7_input.sam
bowtie2 -p 6 -x /media/yanfang/FYWD/CHIPSEQ/fqfile/bowtie2hg19/hg19 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR3101253.fastq.gz -S KYSE510_H3K27Ac.sam
bowtie2 -p 6 -x /media/yanfang/FYWD/CHIPSEQ/fqfile/bowtie2hg19/hg19 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR3101254.fastq.gz -S KYSE510_input.sam

#比對率如下：
64273075 reads; of these:
  64273075 (100.00%) were unpaired; of these:
    1285733 (2.00%) aligned 0 times
    46680241 (72.63%) aligned exactly 1 time
    16307101 (25.37%) aligned >1 times
98.00% overall alignment rate
[bam_sort_core] merging from 16 files and 1 in-memory blocks...
64402758 reads; of these:
  64402758 (100.00%) were unpaired; of these:
    1271927 (1.97%) aligned 0 times
    45320875 (70.37%) aligned exactly 1 time
    17809956 (27.65%) aligned >1 times
98.03% overall alignment rate
[bam_sort_core] merging from 16 files and 1 in-memory blocks...
64555367 reads; of these:
  64555367 (100.00%) were unpaired; of these:
    901311 (1.40%) aligned 0 times
    48336949 (74.88%) aligned exactly 1 time
    15317107 (23.73%) aligned >1 times
98.60% overall alignment rate
[bam_sort_core] merging from 16 files and 1 in-memory blocks...
68947224 reads; of these:
  68947224 (100.00%) were unpaired; of these:
    4955211 (7.19%) aligned 0 times
    45955461 (66.65%) aligned exactly 1 time
    18036552 (26.16%) aligned >1 times
92.81% overall alignment rate
[bam_sort_core] merging from 17 files and 1 in-memory blocks...

5.使用MACS獲得Chip-seq富集區(qū)
以下是官方網(wǎng)站對于常規(guī)峰和“寬”峰的一些解釋(https://www.encodeproject.org/chip-seq/histone/) 。有的文章說常規(guī)峰和寬峰分析過程會有些不一樣羊精，所以我特異查了一下H3K27Ac是什么峰囚玫，根據(jù)這個表來看屬于常規(guī)峰。

• For narrow-peak histone experiments, each replicate should have 20 million usable fragments.
• For broad-peak histone experiments, each replicate should have 45 million usable fragments.
• H3K9me3 is an exception as it is enriched in repetitive regions of the genome. Compared to other broad marks, there are few H3K9me3 peaks in non-repetitive regions of the genome in tissues and primary cells. This results in many ChIP-seq reads that map to a non-unique position in the genome. Tissues and primary cells should have 45 million total mapped reads per replicate.

image.png

• For narrow-peak histone experiments, each replicate should have 10 million usable fragments(https://www.encodeproject.org/data-standards/terms/#read-depth).
• For broad-peak histone experiments, each replicate should have 20 million usable fragments.

$ macs2 callpeak -c /media/yanfang/FYWD/CHIPSEQ/sambam/TE7_input.bam -t /media/yanfang/FYWD/CHIPSEQ/sambam/TE7_H3K27Ac.bam -m 10 30 -p 1e-5 -f BAM -g hs -n TE7_H3K27Ac 2>TE7_H3K27Ac.masc2.log
$ macs2 callpeak -c /media/yanfang/FYWD/CHIPSEQ/sambam/KYSE510_input.bam -t /media/yanfang/FYWD/CHIPSEQ/sambam/KYSE510_H3K27Ac.bam -m 10 30 -p 1e-5 -f BAM -g hs -n KYSE510_H3K27Ac 2>KYSE510_H3K27Ac.masc2.log

-c:對照組。
-t：處理組本昏。
--format:輸入的文件類型涌穆。
--輸出的文件前綴名稱。
--keep-dup對于重復(fù)序列的處理方式趁舀，1效果最好祝沸。指明MACS在染色體同一位置的reads（重復(fù)序列處理方式）。
--wig：輸出的文件類型奉狈。
--single-profile表示輸出單文件仁期。
--space是文獻(xiàn)中的要求。（Wiggle files were generated using read pileups for every 50 base pair bins）
--m:This parameter is used to select the regions within MFOLD range of high-confidence enrichment ratio against background to build model. The regions must be lower than upper limit, and higher than the lower limit of fold enrichment. DEFAULT:10,30 means using all regions not too low (>10) and not too high (<30) to build paired-peaks model. If MACS can not find more than 100 regions to build model, it will use the --shiftsize parameter to continue the peak detection.This parameter is used to select the regions within MFOLD range of high-confidence enrichment ratio against background to build model. The regions must be lower than upper limit, and higher than the lower limit of fold enrichment. 構(gòu)建雙峰模型時使用跛蛋，默認(rèn)[10,30]赊级。表示選擇那些倍數(shù)變化在10-30之間的富集區(qū)域，如果找不到100個區(qū)域構(gòu)建模型橡伞，并且你還設(shè)置了--fix-biomodal,它就會用--shiftsize參數(shù)繼續(xù)分析晋被。

call peak 后得到以下幾個文件類型：

image.png

打開peak.xls文件后墨微，表頭有幾行顯示：

image.png

6.構(gòu)建index文件
為了在IGV里查看peak翘县，需要把bam文件先構(gòu)建index谴分，即生成bam.bai文件：

$ samtools index TE7_H3K27Ac.bam TE7_H3K27Ac.bam.bai
$ samtools index TE7_input.bam TE7_input.bam.bai
$ samtools index KYSE510_H3K27Ac.bam KYSE510_H3K27Ac.bam.bai
$ samtools index KYSE510_input.bam KYSE510_input.bam.bai

順便插一句：[生信人必會之samtools](https://vip.biotrainee.com/d/117-samtools）
這篇文章大概是我看到過寫的最好最全面的samtool使用方法牺蹄。
7.用deeptools里的bamcoverage工具將bam文件標(biāo)準(zhǔn)化

$ bamCoverage -b TE7_H3K27Ac.bam -o TE7_H3K27Ac.bw -p 10 --normalizeUsing RPKM
$ bamCoverage -b TE7_input.bam -o TE7_input.bw -p 10 --normalizeUsing RPKM
$ bamCoverage -b KYSE510_H3K27Ac.bam -o KYSE510_H3K27Ac.bw -p 10 --normalizeUsing RPKM
$ bamCoverage -b KYSE510_input.bam -o KYSE510_input.bw -p 10 --normalizeUsing RPKM

將生成的bw文件（標(biāo)準(zhǔn)化后的）導(dǎo)入IGV：

這是文獻(xiàn)里IRF6基因附近的peak富集情況

這是我自己做的IRF6基因附近的富集情況

這是文獻(xiàn)里TP63基因附近的peak富集情況

這是我自己做的TP63基因附近peak的富集情況

其他文獻(xiàn)里涉及到的IGV的圖氓奈，都幾乎和我做出來的一模一樣鼎天，說明分析過程基本沒啥問題。

7.ROSE篩選super-enhancer
這部分是本次實(shí)踐的重點(diǎn)育勺，ROSE官網(wǎng)：http://younglab.wi.mit.edu/super_enhancer_code.html
有關(guān)這個ROSE軟件講解的比較好的中文版的文章：https://blog.csdn.net/oxygenjing/article/details/77862115
ROSE的最主要用法有ROSE_main.py和ROSE_geneMapper.py涧至。其中ROSE_main.py 用于尋找增強(qiáng)子而ROSE_geneMapper.py 用于為增強(qiáng)子相關(guān)的基因進(jìn)行注釋桑包。

首先看一下ROSE_mian.py 用法：

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

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

也可以像這篇文章一樣https://blog.csdn.net/oxygenjing/article/details/77862115蓖康，安裝ROSE

$ wget https://bitbucket.org/young_computation/rose/get/1a9bb86b5464.zip
$ unzip 1a9bb86b5464.zip

(2)安裝后蒜焊，首先準(zhǔn)備文件：peaks的gff文件。這個peak文件是MACS2 call peak之后生成的peaks文件鳖悠。
官方網(wǎng)站對gff文件的格式要求：
.gff must have the following columns:
1: chromosome (chr#)
2: unique ID for each constituent enhancer region
4: start of constituent
5: end of constituent
7: strand (+,-,.)
9: unique ID for each constituent enhancer region
所以要先從peak文件里提取這些要用的列优妙。這里要注意的是：1,2,4,5,7,9是gff文件的列數(shù)套硼，第3，6九妈，8列不是不要了雾鬼，而是要都寫成“ ”空格形式。
（剛開始我以為其他三列都不要了晶疼，所以在gff文件里只有6列翠霍，結(jié)果一直報錯蠢莺。后來按照某篇教程里寫的其他三列都換成"."，結(jié)果也一直報錯......整整折騰了我3天的時間蒋情，死活跑不通耸携。）

$ awk '{print $1"\t"$10"\t"" ""\t"$2"\t"$3"\t"" ""\t"".""\t"" ""\t"$10}' TE7_H3K27Ac_peaks.bed > TE7_peaks.gff

這里的\t是以制表符分割的意思夺衍。每一列用\t分割，并且用""括起來河劝。

(3)接下來創(chuàng)建一個虛擬的python2.7的環(huán)境，因?yàn)槲已b的是python3牌里，所以需要創(chuàng)建一個虛擬的環(huán)境來運(yùn)行ROSE牡辽。

source activate py27

(4)還需要將samtools的PATH設(shè)置到這個文件夾里敞临，因?yàn)镽OSE運(yùn)行時需要調(diào)用samtools挺尿。如果你的samtool沒有設(shè)置全局變量，是要多一步添加PATH的攀涵。但是因?yàn)槲业膕amtool已經(jīng)在全局變量里洽沟，可以直接在這個文件夾里打"samtools"就可以調(diào)用蜗细，所以我就沒有再設(shè)置環(huán)境變量了炉媒。

(5)之后我們需要bam文件（這個bam文件是需要先sort的吊骤。處理組，對照組）眷细，gff文件（處理組）溪椎，bam.bai文件（也就是bam的index文件沼侣，處理組和對照組）蛾洛。把這些文件全都放進(jìn)ROSE安裝的文件夾里，切記：運(yùn)行ROSE的時候也要在其安裝的文件夾里運(yùn)行扶供，我也不知道這軟件開發(fā)的人是咋想的。扳碍。笋敞。

然后一切準(zhǔn)備好之后，就是見證奇跡的時刻了趁餐。（雖然我這個奇跡到了第4天才看到后雷。贾漏。。之前3.9天都在解決報錯問題了）

python ROSE_main.py -g HG19 -i ./TE7/TE7_peaks.gff -r ./TE7/TE7_H3K27Ac_sorted.bam -c ./TE7/TE7_input_sorted.bam -o roseresult -s 12500 -t 2500

如果運(yùn)行時有報錯盒齿，例如有類似如下顯示：

Traceback (most recent call last):
  File "ROSE_bamToGFF.py", line 249, in <module>
    main()
  File "ROSE_bamToGFF.py", line 240, in main
    newGFF = mapBamToGFF(bamFile,gffFile,options.sense,int(options.extension),options.floor,options.rpm,options.matrix)
  File "ROSE_bamToGFF.py", line 76, in mapBamToGFF
    gffLocus = ROSE_utils.Locus(line[0],int(line[3]),int(line[4]),line[6],line[1])

請確保你的bam文件格式，gff文件格式翎承，samtool的環(huán)境變量都沒有問題的情況下叨咖，并且所有需要的文件在ROSE安裝的文件夾里啊胶，仍然出現(xiàn)這些報錯焰坪，你可以嘗試這樣修改源代碼：（參考官方給出的解決方法：https://bitbucket.org/young_computation/rose/issues/44/typeerror-issue）
FIX:
Make the following changes to ROSE_utils.py:
Line 602 (original): command = 'samtools view %s | head -n 1' % (bamFile)
Line 602 (new): command = './samtools view %s | head -n 1' % (bamFile)
Line 632 (original): command = 'samtools flagstat %s' % (self._bam)
Line 632 (new): command = './samtools flagstat %s' % (self._bam)
Line 657 (original): command = 'samtools view %s %s' % (self._bam,locusLine)
Line 657 (new): command = './samtools view %s %s' % (self._bam,locusLine)

最后運(yùn)行后會生成很多文件：

image.png

輸出文件如下（參考文章:https://blog.csdn.net/oxygenjing/article/details/77862115）：
1.gtf（這個文件在gff文件夾里）： 該文件為輸入gtf文件的副本某饰。

2.stitched.gtf（這個文件也在gff文件夾里）：該文件為通過STITCHING_DISTANCE將INPUT_CONSTITUENT_GFF拼接在一起創(chuàng)建的gff文件；文件列數(shù)如下：
chrom, name, [blank], start, end, [blank], [blank], strand, [blank], [blank], name
其中 name 字段的命名方式為：拼接起來的區(qū)域數(shù)+最左端區(qū)域ID诫尽。

3._MAPPED.gff（在mappedGFF文件夾里）: 每個bam文件通過bamToGFF的輸出文件炬守，包含以下列：
（成分ID减途，測試區(qū)域观蜗，平均讀取密度（單位為每百萬位元每百萬映射的單位讀數(shù)密度））

4.* _STITCHED * _MAPPED.gff（在mappedGFF文件夾里墓捻，處理組和對照組分別有一個文件）： 每個bam文件通過bamToGFF的輸出文件坊夫，該文件中對增強(qiáng)子區(qū)域進(jìn)行了拼接环凿，包含以下列：
（拼接增強(qiáng)子ID，測試區(qū)域羽杰，平均讀取密度（單位為百萬映射每單位拼接增強(qiáng)子數(shù)））

5.STITCHED_ENHANCER_REGION_MAP.txt： 利用bamToGFF計(jì)算后得到的拼接enhancer密度文件，包含以下列：
（拼接增強(qiáng)子ID惕澎，染色體唧喉，拼接增強(qiáng)子起始位置忍抽，拼接增強(qiáng)子末端位置，拼接數(shù)干跛，BAM信號等級驯鳖，BAM信號）

6. AllEnhancers.table.txt：enhancer列表久免，包含每個增強(qiáng)子的排名和是否為超級增強(qiáng)子阎姥，內(nèi)容包括：增強(qiáng)子ID，染色體泽腮，拼接增強(qiáng)子起始位點(diǎn)诊赊，拼接增強(qiáng)子末端鲸郊，拼接數(shù)秆撮，拼接成分大小，BAM的信號拨匆，BAM的等級惭每，是否為超增強(qiáng)子：是（1）否（0）

7._SuperEnhancers.table.txt： 超級enhancer的排名宏赘，為_AllEnhancers.table.txt 文件的子集峻汉。包含以下列：
（拼接增強(qiáng)劑ID扳埂，染色體，拼接增強(qiáng)子起始位點(diǎn)，拼接增強(qiáng)子末端，拼接數(shù)葱淳，縫合在一起的成分的大小途戒，RANKING_BAM的信號喷斋，RANKING_BAM的等級，超增強(qiáng)子的二進(jìn)制（1）與典型（0））

8._Enhancers_withSuper.bed 可以加載到UCSC瀏覽器中可視化的enhancer bed文件

其中那個png圖片就是我們需要的：所有enhancer的散點(diǎn)圖

所有enhancer的散點(diǎn)圖。其他有的教程這個橫坐標(biāo)和縱坐標(biāo)是我的兩倍，不知道是不是因?yàn)橛玫腗ACS版本不一樣造成的裙士。他們用的是MACS做的call peak。我用的是MACS2做的。

之后是KYSE510細(xì)胞系的：

python ROSE_main.py -g HG19 -i ./KYSE510_peaks.gff -r ./KYSE510_H3K27Ac_sorted.bam -c ./KYSE510_input_sorted.bam -o KYSE510_roseresult -s 12500 -t 2500

image.png

8.ROSE注釋enhancer

python ROSE_geneMapper.py -i ./TE7_roseresult/TE7_peaks_AllEnhancers.table.txt -g HG19 -o ./TE7_roseresult 2> ./TE7_roseresult/TE7_enhancer_anno.log
python ROSE_geneMapper.py -i ./KYSE510_roseresult/KYSE510_peaks_AllEnhancers.table.txt -g HG19 -o ./KYSE510_roseresult 2> ./KYSE510_roseresult/KYSE510_enhancer_anno.log

得到以下兩個表格：

image.png

1.ENHANCER_TO_GENE.txt： enhancer重疊基因边器、附近基因以及最近的基因列表
2.GENE_TO_ENHANCER.txt： 以每個基因?yàn)榱忻暮推湎嚓P(guān)的增強(qiáng)子位置信息列表（https://blog.csdn.net/oxygenjing/article/details/77862115）

9.利用R做GO富集分析
這部分我也是新手，不知道過程是否正確，反正是跑通了，對不對的還需再多看看虾宇。旭贬。。
首先在上一步里，我得到了_peaks_AllEnhancers_GENE_TO_ENHANCER.txt這個表，我用的這個表進(jìn)行的GO分析。

>library(clusterProfiler)
>library(DOSE)
>library(org.Hs.eg.db)
#TE7的enhancer GO富集分析
#首先讀取文件吁断，第一行作為列名
>enhancer.relate.gene<-read.table(file="TE7_peaks_AllEnhancers_GENE_TO_ENHANCER.txt",header=T)
#因?yàn)樵谶@個txt里的最后一列注明是否為super enhancer又兵，如果是則為1；如果不是則空白宅粥。所以取為1的所有行，即挑選出所有的super enhancer
>gene<-enhancer.relate.gene[which(enhancer.relate.gene$IS_SUPER=='1'),]
#查看super enhancer的數(shù)據(jù)集
> head(gene)
  GENE_NAME REFSEQ_ID                                PROXIMAL_ENHANCERS ENHANCER_RANKS IS_SUPER
1     LAMA5 NM_005560                           chr20:60924936-60956028              1        1
2     RPS21 NM_001024                           chr20:60924936-60956028              1        1
3 C20orf151 NM_080833                           chr20:60924936-60956028              1        1
4     ADRM1 NM_007002                           chr20:60924936-60956028              1        1
5   CABLES2 NM_031215                           chr20:60924936-60956028              1        1
6    MIR205 NR_029622 chr1:209528209-209591314,chr1:209607925-209609032         2,1441        1
#將這個super enhancer的table第一列提出來，因?yàn)榈谝涣惺腔蛎?gene_name<- gene[,1]
#查看基因名
head(gene_name)
[1] "LAMA5"     "RPS21"     "C20orf151" "ADRM1"     "CABLES2"   "MIR205" 
#調(diào)用R包
library("stringr") 
library("clusterProfiler")
#GO的BP分析
ego_bp <- enrichGO(gene_name,
                   OrgDb      = org.Hs.eg.db,
                   keyType    = 'SYMBOL',
                   ont        = "BP",
                   pAdjustMethod = "BH",
                   pvalueCutoff = 0.001,
                   )
#泡泡圖
dotplot(ego_bp, font.size=10)

泡泡圖

#GO圖
plotGOgraph(ego_bp)

GO圖

#bar圖
barplot(ego_bp,showCategory = 10,title="The GO_BP enrichment analysis of all super enhancers-related genes in TE7")+ 
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(ego_bp) str_wrap(ego_bp,width = 25))

image.png

KYSE510的分析過程和結(jié)果我就不放進(jìn)來了竹习，只需要把輸入的文件換掉就行，其他的一模一樣随夸。