ChIP-seq實(shí)踐(H3K27Ac,enhancer的篩選和enhancer相關(guān)基因的GO分析)

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

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

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

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

image

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

本次實(shí)踐文章：
Targeting super-enhancer-associated oncogenes in oesophageal squamous cell carcinoma
用不同濃度的OSCC抗癌物THZ1處理兩種OSCC細(xì)胞系：KYSE510和TE7。
尋找藥物處理后的超級(jí)增強(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)行一個(gè)文件，都會(huì)在家目錄下的ncbi文件里生成一個(gè)臨時(shí)文件，非常占內(nèi)存。由于我的存儲(chǔ)空間沒有了，所以我就直接下載了fastq文件。下載fastq文件的網(wǎng)址是：https://www.ebi.ac.uk/ena/data/view/SRX1530372。再?gòu)?fù)制下fastq.gz文件的地址。
網(wǎng)上有的文章提到：盡量不要用wget或curl去下載sra文件，某些情況下會(huì)導(dǎo)致下載的文件不完整。所以這次我也嘗試了用ascp下載的方法。
這里說(shuō)一下用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í)時(shí)知道程序在干啥，方便查錯(cuò)。
-T 取消加密，否則有時(shí)候數(shù)據(jù)下載不了。
-i 提供私鑰文件的地址。
-l 設(shè)置最大傳輸速度，一般200m到500m，如果不設(shè)置，反而速度會(huì)比較低，可能有個(gè)較低的默認(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是單端測(cè)序。
我的代碼：

$ 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后面是:號(hào)，不是路徑
最后gz結(jié)尾有一個(gè)空格，然后是"."。如果沒有空格和點(diǎn)，則會(huì)報(bào)錯(cuò)。無(wú)法下載。

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

下面這些代碼是按照上一次實(shí)踐的方法寫的，直接生成bam文件。因?yàn)檫@篇文獻(xiàn)的測(cè)序結(jié)果質(zhì)量很好，所以沒有去掉3'端的幾個(gè)堿基
#!/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比對(duì)后生成sam文件的方法。有些文章里在進(jìn)行比對(duì)的時(shí)候是用的bowtie進(jìn)行比對(duì)的，并且設(shè)置了只保留比對(duì)一次的reads。我這里用的是bowtie2，所以沒有設(shè)置這個(gè)參數(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

#比對(duì)率如下：
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)站對(duì)于常規(guī)峰和“寬”峰的一些解釋(https://www.encodeproject.org/chip-seq/histone/) 。有的文章說(shuō)常規(guī)峰和寬峰分析過(guò)程會(huì)有些不一樣，所以我特異查了一下H3K27Ac是什么峰，根據(jù)這個(gè)表來(lái)看屬于常規(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
$ 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:對(duì)照組。
-t：處理組。
--format:輸入的文件類型。
--輸出的文件前綴名稱。
--keep-dup對(duì)于重復(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)建雙峰模型時(shí)使用，默認(rèn)[10,30]。表示選擇那些倍數(shù)變化在10-30之間的富集區(qū)域，如果找不到100個(gè)區(qū)域構(gòu)建模型，并且你還設(shè)置了--fix-biomodal,它就會(huì)用--shiftsize參數(shù)繼續(xù)分析。

call peak 后得到以下幾個(gè)文件類型：


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

順便插一句：[生信人必會(huì)之samtools](https://vip.biotrainee.com/d/117-samtools）
這篇文章大概是我看到過(guò)寫的最好最全面的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的圖，都幾乎和我做出來(lái)的一模一樣，說(shuō)明分析過(guò)程基本沒啥問(wèn)題。

7.ROSE篩選super-enhancer
這部分是本次實(shí)踐的重點(diǎn)，ROSE官網(wǎng)：http://younglab.wi.mit.edu/super_enhancer_code.html
有關(guān)這個(gè)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倉(cāng)庫(kù)中克隆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文件。這個(gè)peak文件是MACS2 call peak之后生成的peaks文件。
官方網(wǎng)站對(duì)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é)果一直報(bào)錯(cuò)。后來(lái)按照某篇教程里寫的其他三列都換成"."，結(jié)果也一直報(bào)錯(cuò)......整整折騰了我3天的時(shí)間，死活跑不通。）

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

這里的\t是以制表符分割的意思。每一列用\t分割，并且用""括起來(lái)。

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

source activate py27

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

(5)之后我們需要bam文件（這個(gè)bam文件是需要先sort的。處理組，對(duì)照組），gff文件（處理組），bam.bai文件（也就是bam的index文件，處理組和對(duì)照組）。把這些文件全都放進(jìn)ROSE安裝的文件夾里，切記：運(yùn)行ROSE的時(shí)候也要在其安裝的文件夾里運(yùn)行，我也不知道這軟件開發(fā)的人是咋想的。。。

然后一切準(zhǔn)備好之后，就是見證奇跡的時(shí)刻了。（雖然我這個(gè)奇跡到了第4天才看到。。。之前3.9天都在解決報(bào)錯(cuò)問(wèn)題了）

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)行時(shí)有報(bào)錯(cuò)，例如有類似如下顯示：

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])

請(qǐng)確保你的bam文件格式，gff文件格式，samtool的環(huán)境變量都沒有問(wèn)題的情況下，并且所有需要的文件在ROSE安裝的文件夾里，仍然出現(xiàn)這些報(bào)錯(cuò)，你可以嘗試這樣修改源代碼：（參考官方給出的解決方法：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)行后會(huì)生成很多文件：

image.png

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

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

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

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

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

  1. AllEnhancers.table.txt：enhancer列表，包含每個(gè)增強(qiáng)子的排名和是否為超級(jí)增強(qiáng)子，內(nèi)容包括：增強(qiáng)子ID，染色體，拼接增強(qiáng)子起始位點(diǎn)，拼接增強(qiáng)子末端，拼接數(shù)，拼接成分大小，BAM的信號(hào)，BAM的等級(jí)，是否為超增強(qiáng)子：是（1）否（0）

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

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

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

所有enhancer的散點(diǎn)圖。其他有的教程這個(gè)橫坐標(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

得到以下兩個(gè)表格：

image.png

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

9.利用R做GO富集分析
這部分我也是新手，不知道過(guò)程是否正確，反正是跑通了，對(duì)不對(duì)的還需再多看看。。。
首先在上一步里，我得到了_peaks_AllEnhancers_GENE_TO_ENHANCER.txt這個(gè)表，我用的這個(gè)表進(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)樵谶@個(gè)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
#將這個(gè)super enhancer的table第一列提出來(lái)，因?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的分析過(guò)程和結(jié)果我就不放進(jìn)來(lái)了，只需要把輸入的文件換掉就行，其他的一模一樣。

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。
禁止轉(zhuǎn)載，如需轉(zhuǎn)載請(qǐng)通過(guò)簡(jiǎn)信或評(píng)論聯(lián)系作者。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容