ChIP-seq基礎(chǔ)入門

理解ChIP-Seq

到了目前這個(gè)水平，我學(xué)習(xí)新的高通量數(shù)據(jù)分析流程時(shí)已經(jīng)不再考慮代碼應(yīng)該如何寫的問題了泡嘴。我更多要去考慮一個(gè)技術(shù)的目的和意義。

轉(zhuǎn)錄組主要研究的問題是基因在不同情況下的差異表達(dá)以及RNA結(jié)構(gòu)變化等赠群，而表觀組研究的問題是在基因序列不變的情況下离咐，基因的表達(dá)、調(diào)控和性狀發(fā)生了可遺傳變化的分子機(jī)制帮匾。也就是相同的DNA, RNA啄骇，蛋白質(zhì)經(jīng)過一定的修飾后會(huì)使得生物性狀發(fā)生了改變。

Jimmy在生信技能樹的論壇上發(fā)了一個(gè)帖子我對表觀的18個(gè)疑問, 目前他提出了22個(gè)問題瘟斜，還沒有改標(biāo)題缸夹。里面提到的就是ChIP-Seq僅僅是第一個(gè)表觀遺傳學(xué)領(lǐng)域比較成熟的技術(shù)而已，目前還有很多其他的技術(shù)螺句，比如說

  • DNA修飾: DNA甲基化免疫共沉淀技術(shù)（MeDIP), 目標(biāo)區(qū)域甲基化虽惭，全基因組甲基化（WGBS)，氧化-重亞硫酸鹽測序（oxBS-Seq), TET輔助重亞硫酸鹽測序（TAB-Seq）
  • RNA修飾: RNA甲基化免疫共沉淀技術(shù)（MeRIP)
  • 蛋白質(zhì)與核酸相互作用： RIP-Seq, ChIP-Seq, CLIP-Seq
  • 還有最近比較火的ATAC-Seq ATAC-seq能干啥蛇尚？

2013年P(guān)NAS上面有一篇文章叫做 Epigenetics: Core misconcept 講了幾點(diǎn)關(guān)于表觀遺傳學(xué),大家的理解錯(cuò)誤的地方芽唇。

本次主要是分析ChIP-Seq的高通量測序結(jié)果，因此取劫，先介紹什么是ChIP-Seq.
所謂的ChIP-Seq其實(shí)就是把ChIP實(shí)驗(yàn)做完得到的DNA不僅僅用來跑膠匆笤，還送去高通量測序了。重點(diǎn)在于ChIP勇凭，也就是染色體免疫共沉淀（Chromatin Immunoprecipitation）是用來解決什么科學(xué)問題的疚膊。畢竟數(shù)據(jù)分析不是目的，它只是解決問題的一種手段虾标。

ChIP被用來研究細(xì)胞內(nèi)DNA與蛋白質(zhì)相互作用寓盗，具體來說就是確定特定蛋白（如轉(zhuǎn)錄因子）是否結(jié)合特定基因組區(qū)域（如啟動(dòng)子或其它DNA結(jié)合位點(diǎn)）——可能定義順反組。ChIP還被用來確定基因組上與組蛋白修飾相關(guān)的特定位點(diǎn)（即組蛋白修飾酶類的靶標(biāo)） - 來自維基百科璧函。

ChIP的實(shí)驗(yàn)我從來沒有做過傀蚌，所以只能從網(wǎng)上找來大致的流程，簡單分為
第一步： 將蛋白交聯(lián)到DNA上蘸吓。 也就是保證蛋白和DNA能夠結(jié)合善炫，找到互作位點(diǎn)。
第二步： 通過超聲波剪切DNA鏈库继。
第三步： 加上附上抗體的磁珠用于免疫沉淀靶蛋白箩艺〈茏恚抗體很重要
第四步： 接觸蛋白交聯(lián)；純化DNA

實(shí)驗(yàn)流程

如果送去測序就是ChIP-Seq, 后續(xù)就是對應(yīng)分析流程艺谆，分析分為4步：

  • 質(zhì)量控制榨惰， 用到的是FastQC
  • 序列比對，Bowtie2或這BWA
  • peak calling, 建議用MACS
  • peak注釋静汤， 推薦Y叔的ChIPseeker

文章中ChIP-Seq解決的問題

這部分看看這篇文章想解決什么問題琅催，以及ChIP-Seq解決了什么問題。工具一定要為目的服務(wù)虫给。
PS: 可以先去把數(shù)據(jù)下載了藤抡，再看這個(gè)部分

文章的研究內(nèi)容：
PRC1(Polycomb repressive complex 1)與干細(xì)胞命運(yùn)決定有關(guān)。PRC1的組成非常的復(fù)雜抹估，原文這樣寫道：

The protein families that constitute the core of PRC1 contain several members: Cbx (Cbx2, Cbx4, Cbx6, Cbx7, or Cbx8); Ring1A or Ring1B; PHC (PHC1, PHC2, or PHC3); PCGF (PCGF1, PCGF2, PCGF3, PCGF4, PCGF5, or PCGF6); and RYBP or YAF2. Each combination establishes the subtype of PRC1 complexes

文章想解決的問題是：

  1. PRC1復(fù)合體兩類亞型的全基因組定位
  2. 兩類復(fù)合體調(diào)節(jié)的基因表達(dá)是否有差異
  3. Cbx7缠黍，RYBP是否有共同或獨(dú)特的生物學(xué)功能
  4. Cbx7，RYBP在染色質(zhì)的定位是否是相互依賴

得到的部分答案是：小鼠中有兩個(gè)互斥的PRC1復(fù)合體的變異體棋蚌，Cbx7和RYBP嫁佳。Cbx7和RYBP的在基因組上的結(jié)合在同一個(gè)位置，但是同樣互斥谷暮，不能共存蒿往。 Cbx7用于招募Ring1B結(jié)合到染色質(zhì)，RYBP增加PRC1酶活性湿弦。 RYBP所結(jié)合的基因瓤漏，有著較低水平Ring1B和H2AK119ub的水平，表達(dá)量高于Cbx7結(jié)合態(tài)颊埃。RYBP和Cbx7的在基因位點(diǎn)的結(jié)合情況還與代謝調(diào)節(jié)蔬充，細(xì)胞周期過程有關(guān)。

ChIP的用途：
為了解決PRC1亞型的定位問題班利， 作者對PRC1復(fù)合體的組成Cbx7, Ring1B, RYBP, 以及PRC2的Suz12f分別做了ChIP-Seq. 然后就是畫了很多韋恩圖看不同蛋白的靶基因的相互關(guān)系饥漫。還看了不同蛋白在基因的上的位置。

一般而言罗标，ChIP-Seq基本就是得到上面這些圖庸队，根據(jù)疑問再找到某一些基因看調(diào)控蛋白的結(jié)合和定位。本次實(shí)戰(zhàn)也就是嘗試做出這些圖闯割。

根據(jù)原文彻消，ChIP-Seq的流程如下：

  • 先用Bowtie（0.12.7）把RYBP, Cbx7, Ring1B, Suz12,和H2AK11Ub的ChIP-Seq產(chǎn)生短讀(reads) 比對到小鼠的參考基因組上（NCBIM37），參數(shù)允許seed聯(lián)配時(shí)有2個(gè)錯(cuò)配宙拉。
  • 使用MACS（1.4.1）尋找可能的結(jié)合位點(diǎn)宾尚，也就是基因組中大量短讀片段富集的區(qū)域。MACS比較這些區(qū)域在對照組和實(shí)驗(yàn)組的差異谢澈，選擇p值等于10e-5作為閾值煌贴。 之后PDR低于5%的Cbx7, Ring1B, Suz12和H2AK119ub的peak和PDR低于1%的RYBP的peak結(jié)合用于進(jìn)一步分析

這一步看起來就很復(fù)雜御板，之后會(huì)重點(diǎn)進(jìn)行學(xué)習(xí)

  • 在基因內(nèi)部或距離TSS 2.5kb的 peak 被認(rèn)為是靶基因。 每個(gè)TSS附近5Kb區(qū)域崔步，并且與一個(gè)或多個(gè)ChIP-Seq(RYBP,Ring1B, Cbx7, Suz12, H2AK119ub)結(jié)果關(guān)聯(lián)稳吮，用于計(jì)算ChIP-Seq表達(dá)譜和全ChIP-Seq覆蓋度。
  • 使用BEDtools 計(jì)算TSS附近ChIP-Seq表達(dá)譜

ChIP-Seq數(shù)據(jù)分析的流程難點(diǎn)在于找到peak井濒，所以peak calling這一階段的軟件選擇比較重要，目前常用的是MACS2, 實(shí)際分析建議多用幾個(gè)軟件列林。目前常用的是MACS2, 實(shí)際分析建議多用幾個(gè)軟件瑞你。

數(shù)據(jù)下載和質(zhì)量控制

這部分下載ChIP-Seq數(shù)據(jù)，以及小鼠的參考基因組希痴，注釋者甲。

  • 下載測序數(shù)據(jù)，并用sratoolkit解壓縮
    文章提供的GEO是GSE42466砌创，可以在https://www.ncbi.nlm.nih.gov/geo/ 檢索到對應(yīng)的數(shù)據(jù)集
GEO
# 我習(xí)慣將數(shù)據(jù)保存到項(xiàng)目文件夾下的data中
mkdir -p GSE42466/data/chip-seq
cd GSE42466/data/chip-seq
for ((i=204;i<=209;i++));
do
  wget -4 -q ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP017/SRP017311/SRR620$i/SRR620$i.sra;
  ~/miniconda/envs/biostar/bin/fastq-dump –split-3 SRR620$i.sra;
done &
  • 下載小鼠參考基因組的索引和注釋文件, 這里用常用的mm10
# 索引大小為3.2GB虏缸， 不建議自己下載基因組構(gòu)建
mkdir referece && cd reference
wget -4 -q ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip
unzip mm10.zip

注釋文件到https://genome.ucsc.edu/cgi-bin/hgTables 選取保存。

image.png

或者用curl進(jìn)行下載嫩实。

curl 'https://genome.ucsc.edu/cgi-bin/hgTables?hgsid=603641571_MKLxB78UAcjaI3oj5YiFCA0Ms86o&boolshad.hgta_printCustomTrackHeaders=0&hgta_ctName=tb_refGene&hgta_ctDesc=table+browser+query+on+refGene&hgta_ctVis=pack&hgta_ctUrl=&fbQual=whole&fbUpBases=200&fbExonBases=0&fbIntronBases=0&fbDownBases=200&hgta_doGetBed=get+BED' > mm10_refSeq.bed

得到的fastq文件刽辙，一般而言都是上機(jī)得到的原始數(shù)據(jù)，也有可能是經(jīng)過質(zhì)量控制的結(jié)果甲献，但是無論如何都得自己去檢查一下宰缤。

## 需要安裝fastqc 和 multiqc
# 先獲取QC結(jié)果
ls *gz | while read id; do fastqc -t 4 $id; done
# multiqc
multiqc *fastqc.zip --pdf

根據(jù)質(zhì)控結(jié)果，發(fā)現(xiàn)序列的前后幾個(gè)堿基質(zhì)量不好晃洒，因此在比對的時(shí)候會(huì)主動(dòng)忽略慨灭。

QC

今后，可能會(huì)去下載其他實(shí)驗(yàn)的數(shù)據(jù)球及，還要進(jìn)行質(zhì)量控制氧骤，所以就寫一個(gè)Python腳本進(jìn)行流程化操作。腳本命名為sra2qc.py吃引，保存在我的GitHub上

序列比對

這一步非常簡單筹陵， 就是將得到fastq文件用bowtie2比對小鼠參考基因組上，

bowtie2 -p 6 -3 5 --local -x reference/mm10 -U ChIP-Seq/SRR620204.fastq | ~/miniconda3/bin/samtools sort -O bam -o ../analysis/alignment/ring1B.bam
bowtie2 -p 6 -3 5 --local -x reference/mm10 -U ChIP-Seq/SRR620205.fastq | ~/miniconda3/bin/samtools sort -O bam -o ../analysis/alignment/cbx7.bam
bowtie2 -p 6 -3 5 --local -x reference/mm10 -U ChIP-Seq/SRR620206.fastq | ~/miniconda3/bin/samtools sort -O bam -o ../analysis/alignment/suz12.bam
bowtie2 -p 6 -3 5 --local -x reference/mm10 -U ChIP-Seq/SRR620207.fastq | ~/miniconda3/bin/samtools sort -O bam -o ../analysis/alignment/RYBP.bam
bowtie2 -p 6 -3 5 --local -x reference/mm10 -U ChIP-Seq/SRR620208.fastq | ~/miniconda3/bin/samtools sort -O bam -o ../analysis/alignment/IgGold.bam
bowtie2 -p 6 -3 5 --local -x reference/mm10 -U ChIP-Seq/SRR620209.fastq | ~/miniconda3/bin/samtools sort -O bam -o ../analysis/alignment/IgG.bam

我對比對結(jié)果進(jìn)行了簡單的統(tǒng)計(jì)际歼，發(fā)現(xiàn)比對率其實(shí)不太高. 雖然在RNA-Seq是要求80~90%的比對率惶翻，但是在ChIP-Seq上，我沒有具體了解過鹅心。

樣品 比對到單個(gè)區(qū)域 比對到多個(gè)區(qū)域
ring1B 43.76% 15.81%
cbx7 39.39% 17.49%
suz12 48.88% 17.17%
RYBP 58.76% 28.12%
IgGold 54.54% 27.21%
IgG 57.83% 25.45%

比對的結(jié)果吕粗，可以打開IGV查看，和RNA-Seq最大的區(qū)別在于你能看到明顯的峰旭愧，這就是peak颅筋。這類peak在全基因組上特別的多宙暇，我們無法用肉眼直接選擇，因此就需要借助專門的工具查找议泵。

peaks

一般ChIP-Seq都需要提供一個(gè)對照組占贫，也就是圖中IgG。很明顯的發(fā)現(xiàn)其他實(shí)驗(yàn)組有峰的區(qū)域在對照組中是沒有的先口。而實(shí)驗(yàn)組沒有峰的區(qū)域型奥，對照組也沒有，這樣子就能提出背景噪聲碉京。

peak calling

peak calling要用到MACS, 建議用ananconda安裝, 因?yàn)榭梢赃x擇合適的Python環(huán)境厢汹。

conda create -n macs python=2
source activate macs
pip install MACS2

macs2包含一系列的子命令，其中最主要的就是callpeak谐宙， 官方提供了使用實(shí)例

macs2 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01

一般而言烫葬，我們照葫蘆畫瓢，按照這個(gè)實(shí)例替換對應(yīng)部分就行了凡蜻，介紹一下各個(gè)參數(shù)的意義

  • -t: 實(shí)驗(yàn)組的輸出結(jié)果
  • -c: 對照組的輸出結(jié)果
  • -f: -t和-c提供文件的格式搭综，可以是"ELAND", "BED", "ELANDMULTI", "ELANDEXPORT", "ELANDMULTIPET" (for pair-end tags), "SAM", "BAM", "BOWTIE", "BAMPE" "BEDPE" 任意一個(gè)。如果不提供這項(xiàng)划栓，就是自動(dòng)檢測選擇兑巾。
  • -g: 基因組大小， 默認(rèn)提供了hs, mm, ce, dm選項(xiàng)茅姜， 不在其中的話闪朱，比如說擬南芥，就需要自己提供了钻洒。
  • -n: 輸出文件的前綴名
  • -B: 會(huì)保存更多的信息在bedGraph文件中奋姿，如fragment pileup, control lambda, -log10pvalue and -log10qvalue scores
  • -q: q值，也就是最小的PDR閾值素标， 默認(rèn)是0.05称诗。q值是根據(jù)p值利用BH計(jì)算，也就是多重試驗(yàn)矯正后的結(jié)果头遭。
  • -p： 這個(gè)是p值寓免，指定p值后MACS2就不會(huì)用q值了。
  • -m: 和MFOLD有關(guān)计维，而MFOLD和MACS預(yù)構(gòu)建模型有關(guān)袜香，默認(rèn)是5：50，MACS會(huì)先尋找100多個(gè)peak區(qū)構(gòu)建模型鲫惶，一般不用改蜈首，因?yàn)槟悴欢?/li>

原文寫到：MACS estimated the FDR by comparing the peaks obtained from the samples with those from the control samples, using the same p value cutoff (10e-5). 估計(jì)就是指文章用了q=0.05或0.01

macs2 callpeak -c IgGold.bam -t suz12.bam -q 0.05 -f BAM -g mm -n suz12 2> suz12.macs2.log &
macs2 callpeak -c IgGold.bam -t ring1B.bam -q 0.05 -f BAM -g mm -n ring1B 2> ring1B.macs2.log &
macs2 callpeak -c IgG.bam -t cbx7.bam -q 0.05 -f BAM -g mm -n cbx7 2> cbx7.macs2.log &
macs2 callpeak -c IgG.bam -t RYBP.bam -q 0.01 -f BAM -g mm -n RYBP 2>RYBP.macs2.log &

可以看下找到了多個(gè)peak，

wc -l *bed
  2384 cbx7_summits.bed
  8342 ring1B_summits.bed
     0 RYBP_summits.bed
  7619 suz12_summits.bed

結(jié)果顯示RYBP無peak，估計(jì)是數(shù)據(jù)上傳錯(cuò)誤欢策，我們可以下載作者上傳的peak數(shù)據(jù)吆寨。此外ring1B我找了8000多個(gè)，但是原文是7000個(gè)踩寇，或許要修改閾值啄清。

wget ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE42nnn/GSE42466/suppl/GSE42466_Suz12_peaks_10.txt.gz
gzip -d GSE42466_Suz12_peaks_10.txt.gz
mv GSE42466_RYBP_peaks_5.txt RYBP2_summits.bed

結(jié)果文件不只是上述提到的一類，還有如下格式：

  • NAME_peaks.xls: 以表格形式存放peak信息俺孙，雖然后綴是xls辣卒，但其實(shí)能用文本編輯器打開，和bed格式類似睛榄，但是以1為基添寺，而bed文件是以0為基.也就是說xls的坐標(biāo)都要減一才是bed文件的坐標(biāo)
  • NAME_peaks.narrowPeak NAME_peaks.broadPeak 類似。后面4列表示為懈费， integer score for display， fold-change博脑，-log10pvalue憎乙，-log10qvalue，relative summit position to peak start叉趣。內(nèi)容和NAME_peaks.xls基本一致泞边，適合用于導(dǎo)入R進(jìn)行分析。
  • NAME_summits.bed：記錄每個(gè)peak的peak summits疗杉，話句話說就是記錄極值點(diǎn)的位置阵谚。MACS建議用該文件尋找結(jié)合位點(diǎn)的motif。
  • NAME_model.r烟具，能通過$ Rscript NAME_model.r作圖梢什，得到是基于你提供數(shù)據(jù)的peak模型。

這一部分的工作主要是學(xué)習(xí)MACS2軟件使用方法和參數(shù)的含義朝聋，難度不大嗡午。目標(biāo)是找到候選靶基因。

結(jié)果注釋和可視化

當(dāng)我們得到上一步的輸出文件后冀痕，下一步就是得到文章的結(jié)果荔睹。用的是Y叔的ChIPseeker，一個(gè)充滿故事的R包言蛇，搜索一下就能看到Y(jié)叔的往事僻他。

ChIPseeker的功能分為三類:

  • 注釋：提取peak附近最近的基因， 注釋peak所在區(qū)域
  • 比較：估計(jì)ChIP peak數(shù)據(jù)集中重疊部分的顯著性腊尚；整合GEO數(shù)據(jù)集吨拗，以便于將當(dāng)前結(jié)果和已知結(jié)果比較
  • 可視化： peak的覆蓋情況；TSS區(qū)域結(jié)合的peak的平均表達(dá)譜和熱圖；基因組注釋丢胚；TSS距離翩瓜；peak和基因的重疊

完全符合我們的需要，都不需要用到bedtools和deeptools携龟。
在正式的工作開始之前兔跌，我們需要加載ChIPseeker, 基因組注釋包和bed數(shù)據(jù)

# biocLite("ChIPseeker")
# biocLite("org.Mm.eg.db")
# biocLite("TxDb.Mmusculus.UCSC.mm10.knownGene")
library("ChIPseeker")
# 根據(jù)自己存放文章具體地址修改，保證導(dǎo)入的數(shù)據(jù)是BED格式
ring1B <- readPeakFile("F:/Data/ChIP/ring1B_peaks.narrowPeak")
cbx7 <- readPeakFile("F:/Data/ChIP/cbx7_peaks.narrowPeak")
RYBP <- readPeakFile("F:/Data/ChIP/RYBP2_summits.bed")
suz12 <- readPeakFile("f:/Data/ChIP/suz12_peaks.narrowPeak")

得到最基本的peak數(shù)據(jù)后峡蟋，我們要做以下事情

  • 大量的韋恩圖坟桅，描述不同蛋白結(jié)合基因的關(guān)系
  • 不同蛋白在染色體結(jié)合位點(diǎn)與TSS的距離

為了完成上面說的兩件事情，我們需要注釋peak蕊蝗，準(zhǔn)備TSS所在位置然后把peak比對到這些區(qū)域.
注： 原文提供不同蛋白的交集的基因TSS區(qū)域仅乓，所以要先做注釋。

關(guān)于peak注釋一定要多說幾句：

  • 啟動(dòng)子區(qū)域目前沒有明確定義蓬戚，一般認(rèn)為基因起始轉(zhuǎn)錄位點(diǎn)的上下游區(qū)域是啟動(dòng)子區(qū)域夸楣，文章認(rèn)為是TSS2.5 kb 以內(nèi)都是靶基因
  • 注釋有兩類，genomic annotation和nearest gene annotation. 前者是研究直接調(diào)控對象子漩，后者是做基因表達(dá)調(diào)控
  • 還有一類注釋就是簡單看peak上下游的范圍的基因

根據(jù)文章得知需要距離TSS一定距離的注釋豫喧，以及peak上下5 Kb 的注釋

Genes with a peak within the gene body, or within 2.5 kb from the TSS, were considered to be target genes. An area of 5 kb surrounding each TSS that was associated to one or more ChIP-seq (RYBP, Ring1B, Cbx7, Suz12, or H2AK119ub) was used to calculate the ChIP-seq profile and the whole ChIP-seq coverage.

ChIPseeker負(fù)責(zé)注釋的就是annotatePeak，提供了3類注釋方式幢泼，符合文章的要求

# 對于非向量的循環(huán)紧显，推薦構(gòu)建list，然后用lapply
peaks <- list(ring1B, cbx7, RYBP, suz12)
# 注釋
peakAnnoList <- lapply(peaks, annotatePeak, tssRegion=c(-2500,2500), TxDb=txdb, 
                       addFlankGeneInfo=TRUE, flankDistance=5000)

結(jié)果可以用as.GRanges或者as.data.frame查看缕棵。

as.GRanges(peakAnnoList[[1]])
as.data.frame(peakAnnoList[[1]])

得到的結(jié)果會(huì)在原先的bed文件中增加額外的注釋信息孵班，你可以和之前導(dǎo)入的數(shù)據(jù)進(jìn)行比較，比如說geneId, transcriptId,distanceToTSS,flank_txIds等招驴。這些都保存后續(xù)作圖所需要的信息篙程。我們可以先看不同蛋白的靶位點(diǎn)的注釋情況，

plotDistToTSS(peakAnnoList)
image.png

原文大量的是不同蛋白結(jié)合基因的韋恩圖忽匈，我們可以用更有趣的UpSet的點(diǎn)圖

genes= lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
library(UpSetR)
upset(fromList(genes))
不同蛋白靶基因的集合關(guān)系

UpSet技術(shù)適用于多于5個(gè)集合表示情況房午，對于兩個(gè)蛋白的靶基因的比較，大家更加喜歡韋恩圖的丹允。

vennplot(genes[1:2])
簡單的韋恩圖

感覺沒有原文那么高大上郭厌，主要原因就是沒有上色，可以自行搜索找合適的軟件作圖雕蔽。關(guān)于不同蛋白靶基因的集合關(guān)系圖基本上就是這樣子折柠。

下面是根據(jù)靶基因間的交集，繪制TSS距離圖批狐，比如說四個(gè)蛋白的重疊基因扇售。

# 準(zhǔn)備TSS區(qū)域前塔，然后將peak映射到這些區(qū)域，得到tagMatrix
fourgenes <- intersect(genes[[1]],intersect(genes[[2]],intersect(genes[[3]],genes[[4]])))
target <- select(txdb, keys=fourgenes, columns=c("TXID","GENEID"), keytype = "GENEID")
promoter <- promoters(txdb, upstream = 2500, downstream = 2500)
target_promoter <- promoter[promoter$tx_id %in% target$TXID]
tagMatrixList <- lapply(peaks, getTagMatrix, window=target_promoter)
plotAvgProf(tagMatrixList , xlim = c(-3000,2999))
TSS距離圖

感覺和原圖有那么一點(diǎn)像承冰，還可以畫個(gè)熱圖

tagHeatmap(tagMatrixList, xlim=c(-2999, 3000), color=NULL)
熱圖

基本上原文的圖华弓，我都能畫出來，至于想不想就不管了困乒。

推薦閱讀

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 人面猴
    序言：七十年代末娜搂，一起剝皮案震驚了整個(gè)濱河市迁霎，隨后出現(xiàn)的幾起案子，更是在濱河造成了極大的恐慌百宇，老刑警劉巖考廉，帶你破解...
    沈念sama閱讀 206,013評論 6贊 481
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件，死亡現(xiàn)場離奇詭異携御，居然都是意外死亡昌粤，警方通過查閱死者的電腦和手機(jī)，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,205評論 2贊 382
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進(jìn)店門啄刹，熙熙樓的掌柜王于貴愁眉苦臉地迎上來婚苹，“玉大人，你說我怎么就攤上這事鸵膏。” “怎么了怎炊？”我有些...
    開封第一講書人閱讀 152,370評論 0贊 342
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵谭企，是天一觀的道長。 經(jīng)常有香客問我评肆，道長债查，這世上最難降的妖魔是什么？ 我笑而不...
    開封第一講書人閱讀 55,168評論 1贊 278
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任瓜挽，我火速辦了婚禮盹廷，結(jié)果婚禮上，老公的妹妹穿的比我還像新娘久橙。我一直安慰自己俄占，他們只是感情好，可當(dāng)我...
    茶點(diǎn)故事閱讀 64,153評論 5贊 371
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布淆衷。 她就那樣靜靜地躺著缸榄，像睡著了一般。 火紅的嫁衣襯著肌膚如雪祝拯。 梳的紋絲不亂的頭發(fā)上甚带，一...
    開封第一講書人閱讀 48,954評論 1贊 283
  • 城市分裂傳說
    那天，我揣著相機(jī)與錄音，去河邊找鬼鹰贵。 笑死晴氨，一個(gè)胖子當(dāng)著我的面吹牛，可吹牛的內(nèi)容都是我干的碉输。 我是一名探鬼主播籽前，決...
    沈念sama閱讀 38,271評論 3贊 399
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼，長吁一口氣：“原來是場噩夢啊……” “哼腊瑟！你這毒婦竟也來了聚假？” 一聲冷哼從身側(cè)響起，我...
    開封第一講書人閱讀 36,916評論 0贊 259
  • 萬榮殺人案實(shí)錄
    序言：老撾萬榮一對情侶失蹤闰非，失蹤者是張志新（化名）和其女友劉穎膘格，沒想到半個(gè)月后，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體财松，經(jīng)...
    沈念sama閱讀 43,382評論 1贊 300
  • ?護(hù)林員之死
    正文 獨(dú)居荒郊野嶺守林人離奇死亡瘪贱，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 35,877評論 2贊 323
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了辆毡。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片菜秦。...
    茶點(diǎn)故事閱讀 37,989評論 1贊 333
  • 活死人
    序言：一個(gè)原本活蹦亂跳的男人離奇死亡，死狀恐怖舶掖，靈堂內(nèi)的尸體忽然破棺而出球昨，到底是詐尸還是另有隱情，我是刑警寧澤眨攘，帶...
    沈念sama閱讀 33,624評論 4贊 322
  • ?日本核電站爆炸內(nèi)幕
    正文 年R本政府宣布主慰，位于F島的核電站，受9級特大地震影響鲫售，放射性物質(zhì)發(fā)生泄漏共螺。R本人自食惡果不足惜，卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 39,209評論 3贊 307
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一情竹、第九天 我趴在偏房一處隱蔽的房頂上張望藐不。 院中可真熱鬧，春花似錦秦效、人聲如沸雏蛮。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,199評論 0贊 19
  • 一樁弒父案阱州，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽底扳。三九已至，卻和暖如春贡耽，著一層夾襖步出監(jiān)牢的瞬間衷模，已是汗流浹背鹊汛。 一陣腳步聲響...
    開封第一講書人閱讀 31,418評論 1贊 260
  • 情欲美人皮
    我被黑心中介騙來泰國打工， 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留阱冶，地道東北人刁憋。 一個(gè)月前我還...
    沈念sama閱讀 45,401評論 2贊 352
  • 代替公主和親
    正文 我出身青樓，卻偏偏與公主長得像木蹬，于是被迫代替她去往敵國和親至耻。 傳聞我的和親對象是個(gè)殘疾皇子，可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 42,700評論 2贊 345

推薦閱讀更多精彩內(nèi)容