精华内容
下载资源
问答
  • ChIP-seq

    千次阅读 2016-07-05 20:57:14
    ChIP-seq

    ChIP-seq

    What is ChIP-Seq? Chromatin Immunoprecipitation coupled with tiling microarrays or Next-Generation Sequencing NGS. Main goal is to understand genome-wide protein-DNA interaction. Determine location of proteins bound to DNA. ChIP-seq is useful for detecting transcription factor binding sites and histone modification patterns. Common questions we have are which genes is this TF regulating and how do histone modifications affect expression.

    ChIP-seq workflow:
    这里写图片描述

    sequencing & alignment:

    Sequencing depth rules of thumb: > 10M reads for narrow peaks, > 20M for broad peaks
    Long & paired end useful but not essential – alignment in ambiguous regions
    Basic aligners generally adequate, e.g., no need to align splice junctions

    peak calling:

    Very large number of peak calling programs; some specialized for e.g., narrow vs. broad peaks. Commmonly used: MACS, PeakSeq, CisGenome, …
    Model based analysis of ChIP-Seq (MACS) is a ‘peak calling’ method for identifying genomic binding sites from read count data. It was developed in X. Shirley Liu’s lab at Dana Farber Cancer Institute. It uses condent peaks to model shift size. Dynamic Poisson distribution to model the tag distribution along the genome. There are several features available in MACS to conduct ChIP-Seq analysis:

    • call peak: Call peaks from alignment results.
    • bdgpeakcall: Call peaks from bedGraph output.
    • bdgbroadcall: Call broad peaks from bedGraph output.
    • bdgcmp: Deduct noise by comparing two signal tracks in bedGraph.
    • bdgdi : Dierential peak detection based on paired four bed graph les.
    • flterdup : Remove duplicate reads at the same position then convert acceptable format to BED format.
    • predicted: Predict d or fragment size from alignment results.
    • pileup: Pileup aligned reads with a given extension size fragment size or d in MACS language.
    • randsample: Randomly sample number/percentage of total reads.
    • refine peak : Experimental Take raw reads alignment rene peak summits and give scores measuring balance of forward- backward tags. Inspired by SPP.

    An example of using MACS for ChIP-Seq Data ChIP-Seq data from an experiment investigating the circadian recruitment of histone deacetylase 3 HDAC3 and its eect on liver metabolism.

    • Filter duplicates: remove duplicates to mitigate amplication bias:
    srun -n 1 --pty -p interact --mem2000 macs2 filterdup -i /n/stat115/labs/6/HDAC3.bam -g mm --keep-dup 1 -o HDAC3.bed 
    srun -n 1 --pty -p interact --mem2000 macs2 filterdup -i /n/stat115/labs/6/control.bam-g mm --keep-dup 1 -o control.bed  #control
    -i Intput file 
    -g Genome size 
    --keep-dup Duplicates to keep 
    -o Output file.
    • Sample down: Imbalance in read counts between treatment and control can bias results. To sample down the condition that has more reads:
    srun -n 1 --pty -p interact --mem2000 macs2 randsample -t control.bed -n 8757629 -o controlsamp.bed
    
    -t Intput file
    -n Number to sample
    -o Output file

    Installing Loading MACS
    1. Look at install page from Tao Liu’s GitHub page
    2. You can import it as a regular python package OR load it in a cluster environment.

    Example Cont’d Call peaks:

    srun -n 1 -t 60 --pty -p interact --mem2000 macs2 callpeak -t HDAC3.bed -c controlsamp.bed -f BED -g mm -n HDAC3peaks
    
    -t Treatment file
    -c Control file
    -f File format
    -n Output prefix

    Quality Assessment: ChIPQC
    Inputs: BAM files (raw data) and BED files (called peaks)

    experiment <- ChIPQC(samples) 
    ChIPQCreport(experiment) 

    Output: HTML report

    Downstream Analysis:

    • Annotation: what genes are my peaks near?
    • Differential representation: which peaks are over- or under-represented in treatment 1, compared to treatment 2?
    • Motif identification (peaks over known motifs?) and discovery
    • Integrative analysis, e.g., assoication of regulatory elements and expression

    Annotation: ChIPpeakAnno:
    Inputs I
    Peaks: RangedData (GRanges-like) peaks, e.g., from rtracklayer::import() BED files
    Annotation: RangedData representing gene boundaries, or query to biomaRt library(ChIPpeakAnno)

    annotated <- annotatePeakInBatch(peaks, AnnotationData=annotation) 

    Output: RangedData with annotations about near-by peaks.

    Differential Representation: DiffBind
    Inputs: called peaks and raw BED or BAM files

    library(DiffBind) 
    tamoxifen = dba(sampleSheet="tamoxifen.csv") 
    tamoxifen = dba.count(tamoxifen) 
    tamoxifen = dba.contrast(tamoxifen, categories=DBA_CONDITION) 
    tamoxifen = dba.analyze(tamoxifen) tamoxifen.DB = dba.report(tamoxifen)

    Outputs: diagnositics, visiualizations, and ‘top table’ of differentially expressed regions.

    ChIP-seq in Bioconductor: packages

    • Quality assessment – ChIPQC;
    • (Peak calling) – chipseq, PICS, triform, ChIPseqR, iSeq, …
    • Single sample summary / exploration – ChIPpeekAnno, chIPseeker
    • Differential representation – DiffBind, MMDiff , …
    • Motifs – MotifDb, TFBSTools (matching known motifs), motifRG, - MotIV, rGADEM BCRANK (motif discovery)
    • Integration with expression data – Rcade, epigenomix
    展开全文
  • 我开发了一个基于Snakemake的ChIP-seq管道: 。 和ATACseq管道: ChIP-seq的资源 : :来自ENCODE的元数据的汇编。 一个bioc包,用于访问ENCODE的元数据并下载原始文件。 论文: 。 序列为.sra格式,需要使用...
  • git clone https://github.com/niekwit/ChIP-Seq-pipeline.git 除了Picard(可以将其安装在/ home中的任何位置)之外,所有依赖项都应包含在环境变量中。 align设置 HISAT2和bwa索引文件必须从fasta文件中手动...
  • ChIP-Seq-开源

    2021-04-28 00:33:23
    ChIP-Seq软件提供了用于分析ChIP-seq数据和其他类型的质量基因组注释数据的方法。 最常见的分析任务包括位置相关性分析,峰检测以及将基因组划分为信号丰富和信号耗尽的区域。
  • 工作流程 自动化的ChIP-seq和ATAC-seq工作流程
  • MACS:ChIP-Seq的基于模型的分析 最新发布: Github: PyPI: Bioconda: Debian Med: 介绍 随着测序技术的改进,染色质免疫沉淀后再进行高通量测序(ChIP-Seq)成为研究全基因组蛋白质-DNA相互作用的工具。 ...
  • lcdb-wf 进行常见的高通量测序分析的工作流和工具的集合,以及相关的基础结构。 请参阅文档。 参考工作流程 DAG参考资料工作流程: RNA-seq工作流程 DAG用于RNA-seq工作流程: ...DAG for ChIP-seq工作流程:
  • 该存储库包括用于处理和分析 ChIP-Seq 数据的脚本。 bed2frip.sh 用于计算床格式中给定峰轮廓的读数分数 (FRiP) 的脚本。 有关此质量指标的其他信息,请参阅: ENCODE 和 modENCODE 联盟的 ChIP-seq 指南和实践。 ...
  • 结合分辨率放大器和合作互动定位器(BRACIL)是一种盲反卷积模型,将ChIP-seq覆盖范围与模体发现相结合。 BRACIL可在ChIP-seq数据中以高分辨率以及协作相互作用来预测结合位点的位置。
  • ChIP-seq analysis pipeline

    2019-04-11 09:32:24
    ChIP-seq analysis pipeline, bwa:mapping, samtools: samtobam, macs2:callpeak, bdgtobw, homer: findmotif; bedtools: annotation
  • CHANCE-HT 是一种 ChIP-seq 预处理软件,可过滤 IP 强度较弱的样本,识别严重偏向的对照实验和测序样本不足,检测批次效应,并对大型 ChIP-seq 数据集进行归一化。 CHANCE-HT 使用并行处理方法来规范化和过滤大量 ...
  • Chip-seq流程报告

    万次阅读 热门讨论 2017-02-01 21:17:01
    实验旨在了解Chip-seq的基本原理。通过模仿文献《Targeting superenhancerassociated oncogenes in oesophageal squamous cell carcinoma》的流程,学会利用NCBI和EBI数据库下载数据,熟悉Linux下的基本操作,并使用...

     

    一、摘要

    实验旨在了解Chip-seq的基本原理。通过模仿文献《Targeting super enhancer associated oncogenes in oesophageal squamous cell carcinoma》的流程,学会利用NCBI和EBI数据库下载数据,熟悉Linux下的基本操作,并使用R语言画图,用Python或者shell写脚本进行基本的数据处理,通过FastQC、Bowtie、Macs、samtools、ROSE等软件进行数据处理,并对预测结果进行分析讨论。

    二、材料

    1、硬件平台

    处理器:Intel(R) Core(TM)i7-4710MQ CPU @ 2.50GHz 2.50GHz

    安装内存(RAM):16.0GB

    2、系统平台

    Windows 8.1,Ubuntu

    3、软件平台

    ① Aspera connect ② FastQC ③ Bowtie

    ④ Macs 1.4.2 ⑤ IGV ⑥ ROSE

    4、数据库资源

    NCBI数据库:https://www.ncbi.nlm.nih.gov/;

    EBI数据库:http://www.ebi.ac.uk/;

    5、研究对象

    加入H3K27Ac 抗体处理过的TE7细胞系测序数据和其空白对照组

    加入H3K27Ac 抗体处理过的KYSE510细胞系和其空白对照组

    背景简介:食管鳞状细胞癌(OSCC)是一种侵袭性的恶性肿瘤,本文章通过高通量小分子抑制剂进行筛选,发现了一个高度有效的抗癌物,特异性的CDK7抑制剂THZ1。RNA-Seq显示,低剂量THZ1会对一些致癌基因的产生选择性抑制作用,而且,对这些THZ1敏感的基因组功能的进一步表征表明他们经常与超级增强子结合(SE)。ChIP-seq解读在OSCC细胞中,CDK7的抑制作用的机制。

    本文亮点:确定了在OSCC细胞中SE的位置,以及识别出许多SE有关的调节元件;并且发现小分子THZ1特异性抑制SE有关的转录,显示强大的抗癌性。

    文章PMID: 27196599

    三、方法

    实验数据获取流程如下:

    实验数据获取流程

    数据分析流程图如下:

    整体流程图

    1、Aspera软件下载及安装

    进入Aspera官网的Downloads界面,选中aspera connect server,点击Wwindows图标,选择v3.6.2版本,点击Download进行下载。

     

    图表 1 aspera的下载

    Linux下的安装配置参考博文:

    http://blog.csdn.net/likelet/article/details/8226368

    2、Chip-Seq数据下载

    1)选择NCBI的GEO DataSets数据库,输入GSE76861,打开GSM2039110、GSM2039111、2039112、GSM2039113获取它们对应的SRX序列号。

     

    图表 2 Chip-seq数据

     

    图表 3 获取SRA编号

    2)进入EBI,获取ascp下载地址

     

    图表 4 ascp下载地址

    3)使用aspera下载并解压

    aspera下载命令及gunzip解压命令(nohup+命令+&可以后台运行)

     

    3、FastQC质量检查

    3.1 FastQC的安装

    Ubuntu软件包内自带Fastqc

    故安装命令apt-get install fastqc

    3.2 使用FastQC进行质量检查

    fastqc命令:

    fastqc -o . -t 5 -f fastq SRR3101251.fastq &

    -o . 将结果输出到当前目录

    -t 5 表示开5个线程运行

    -f fastq SRR3101251.fastq 表示输入的文件

    (要分别对四个fastq文件执行四次)

    4、使用Bowtie对Reads进行Mapping

    4.1 Bowtie的安装

    Ubuntu软件包内自带bowtie

    故安装命令apt-get install bowtie

    4.2 下载人类参考基因组

    文献说序列比对到了人类参考基因组GRCh37/hg19上

    bowtie官网上面有人类参考基因组hg19已经建好索引的文件

     

    图表 5 bowtie hg19建好的索引

    再执行解压缩命令:unzip hg19.ebwt.zip

    4.3 使用bowtie进行比对

    bowtie命令:

    5、MACS寻找Peak富集区

    5.1 Macs14的安装

    至刘小乐实验室网站下载http://liulab.dfci.harvard.edu/MACS/Download.html

     

    解压后,切换到文件夹目录,执行

    python setup.py install

    5.2 使用Macs建模,寻找Peaks富集区

    MACS命令:

     

    6、IGV可视化

    6.1数据正规化normalised

    编写python程序对wig文件进行normalised

     

    对TE7_H3K27Ac和KYSE510_H3K27Ac的wig文件(即MACS后生成的treat文件夹里的wig文件)计算RPM

    RPM公式:(某位置的reads数目÷所有染色体上总reads数目)×1000000

    6.2 使用wigToBigWig转化格式

     

    6.3安装IGV(Integrative Genomics Viewer)对结果可视化

    从IGV官网下载windows版本http://software.broadinstitute.org/software/igv/download根据提示安装

    直接点击打开igv.jar或者对bat文件以管理员身份运行

    首先,载入hg19基因组;接着载入两个normalised后的bw文件即可

    7、ROSE鉴定Enhancer

    7.1 ROSE程序安装

    ROSE程序可以到http://younglab.wi.mit.edu/super_enhancer_code.html下载,并且有2.7G的示例数据

    7.2 数据预处理

     

    7.3运行ROSE程序

     

    7.4 进行基因注释

     

    7.5 编写R程序,绘制Enhancer及邻近基因

     

    图表 6 TE7.r程序

     

    图表 7 KYSE510.r程序

     

    四、结果

    1、Chip-Seq数据下载

    Chip-Seq数据下载并解压结果

     

    图表 8 Chip-Seq数据

     

    2、FastQC质量检查

    数据质量检查

     

    图表 9 质量检查文件

     

    图表 10 质量检查结果

     

    3、使用Bowtie对Reads进行Mapping

    3.1基因组文件

     

    图表 11人类参考基因组HG19索引

    3.2 Mapping结果

     

    图表 12 Mapping整体结果

    图表 13 生成的sam文件

     

    4、MACS寻找Peak富集区

    4.1MACS结果文件

     

    图表 14 TE7实验对照组结果

    图表 15 KYSE510实验对照组结果

    4.2 MACS结果解读

    Peaks.xls从左至右依次是:峰所在的染色体名称,峰的起始位置,峰的结束为止,峰的长度,峰的高度,贴上的reads标签个数,pvalue(表示置信度),峰的富集程度,FDR假阳性率(越小则峰越好)

     

    图表 16 Peaks.xls文件

    negative_peaks.xls当有对照组实验存在时,MACS会进行两次peak calling。第一次以实验组(Treatment)为实验组,对照组为对照组,第二次颠倒,以实验组为对照组,对照组为实验组。这个相当于颠倒过后计算出来的文件

     

    图表 17 negative_peaks.xls

    Peaks.bed文件相当于Peaks.xls的简化版,从左至右依次是:峰所在的染色体名称,峰的起始位置,峰的结束为止,峰的MACS名称,pvalue(表示置信度)

    图表 18  Peaks.bed文件

    summits.bed是峰顶文件,从左至右依次是:峰所在的染色体名称,峰顶的位置,峰的MACS名称,峰的高度

     

     

    图表 19 summits.bed文件

    MACS_wiggle文件夹下面分为control文件夹和treat文件夹,里面分别存了control组和treat组每隔50bp,贴上的reads数目。第一列为染色体上的位置;第二列为从第一列对应的位置开始,延伸50bp,总共贴上的标签(reads)个数。

     

    图表 20 wiggle文件夹下afterfiting_all.wig文件

    model.r文件可以使用R运行,绘制双峰模型的图片PDF

     

    图表 21 model.r文件

     

    图表 22 TE7双峰模型   图表 23 KYSE510双峰模型

     

    5、IGV对peaks可视化

    5.1Normalised后,wig文件与文献数据比较

     

    图表 24 peaks整体统计比较

    5.2 IGV peaks整体可视化

     

    图表 25 IGV可视化

    6、ROSE分析结果

    6.1 数据预处理结果

    Samtools将sam文件转化为bam文件,并且排序,再建立索引

     

    图表 26 bam文件和bai索引

    6.2 ROSE程序Enhancer分类结果

     

    图表 27 TE7 Enhancer分类结果

     

    图表 28 KYSE510 Enhancer分类结果

     

    peaks_AllEnhancers.table.txt文件从左到右分别是,Enhancer区域名称ID,染色体位置,Enhancer起始位置,结束位置,由多少个Enhancer缝合连接而成,Enhancer大小,Treat组峰高度,Control组峰高度,Enhancer大小排名,是否为Super Enhancer

     

    图表 29 peaks_AllEnhancers.table.txt文件

    peaks_Plot_points.png图片,纵坐标为peaks_AllEnhancers.table.txt中G,H列相减结果,及减掉对照组峰后的高度,横坐标为全部Enhancer的排名,越可能是SuperEnhancer则越靠图的右边。

     

    图表 30 TE7_peaks_Plot_points.png图表 31 KYSE510_peaks_Plot_points.png 

    6.3 基因注释结果

    AllEnhancers_ENHANCER_TO_GENE.txt第J列开始为离Enhancer最近的基因名称

    AllEnhancers_GENE_TO_ENHANCER.txt第1列为基因名,后面为邻近峰的名称

     

    图表 32 AllEnhancers_ENHANCER_TO_GENE.txt文件

     

    图表 33 AllEnhancers_GENE_TO_ENHANCER.txt

    五、讨论和结论

    1、结论

    1.1 FastQC质量检查

    FastQC 版本和机房小型机不同,为v0.10.1,因此检测结果略有区别。图表 8 质量检查结果显示,测序质量挺好,Per base sequence content、Per sequence GC content、Kmer Content出现警告更可能是由于测序方法本身存在的固有误差。

    1.2 bowtie整体覆盖度

    由图表 10 Mapping整体结果可以看出,四个fastq文件Mapping整体覆盖率都在90%以上,从另一方面说明数据质量很好

     

    1.3 ROSE辨别出的Super Enhancer

    由图表 29 TE7_peaks_Plot_points.png图表 28 KYSE510_peaks_Plot_points.png可以看出,在TE7细胞系中,找出了439个Super Enhancer,在KYSE510细胞系中,找出了823个Super Enhancer。 

     

    2、讨论

    由IGV可视化图可以看出,峰的高度和位置基本和文献相同。

     

    图表 34 IGV可视化图

    再用R程序根据ROSE程序结果,绘制和文献相同的图片,与文献的图片进行比较,可以看出来,基因的分布是相似的,就是具体位置和文献不是很一样。

     

    图表 35 本流程结果

     

     

    图表 36 文献结果

    在MACS结果中,有些很窄的峰高度明显比文献要低,这可能是因为bowtie时候,设置的参数使得多条reads比对上仅输出一次,使得峰高度减小。

    在ROSE结果中,MIR205HG没有标注出来,而文献中有此基因,经过检查,在相似位置ROSE程序有找到MIR205基因,这可能是基因注释文件和文献不同导致的。

     

    参考文献

    [1] Targeting super-enhancer-associated oncogenes in oesophageal squamous cell carcinoma PMID: 27196599

    展开全文
  • 1. 匹配不同特征的 Chip-Seq Data 2. 匹配 Chip-Seq Data 和输入数据 3. 查找 Seq 信息,方便 Chip-seq 数据的下载和分析
  • RNA-seqChIP-seq数据分析:课程资料 数据和会话设置 资料呈现 会话设置 序列,注释和索引 基因组序列(fasta) 注释(GTF文件): STAR指数 Bowtie2指数 笔录序列 原始数据(读取) RNA序列 原始读取-质量控制-...
  • DROMPA(DRaw和观察多种浓缩曲线和注释)是一种ChIP-seq管道工具,可满足各种需求,包括质量检查,分析和多个ChIP样品的可视化。 DROMPAplus用C ++编写,具有许多有价值的功能。 DROMPAplus: 接受多种地图文件...
  • ChiP-Seq-分析-复制 该项目是ChiP-Seq分析的复制,该实验是关于由独特表观遗传变化介导的终末红细胞生成过程中基因诱导和抑制的实验的 你可以在找到研究项目 数据 使用了来自项目数据的以下样本: 样品名称 GEO加入...
  • ChIP-seq 数据分析

    万次阅读 2018-09-11 18:04:46
    1 ChIP-Seq技术 1.1 概念 1.2 ChIP-seq技术原理 2 ChIP-Seq数据分析 2.1 数据下载 2.2 质量控制(data_assess) 2.2.1 质量评估 2.2.2 清理 reads 2.3 比对到参考基因组(mapping_analysis) fastq→sam→bam...

    ChIP-Seq仅仅是第一个表观遗传学领域比较成熟的技术而已,目前还有很多其他的技术,比如说

    DNA修饰: DNA甲基化免疫共沉淀技术(MeDIP), 目标区域甲基化,全基因组甲基化(WGBS),氧化-重亚硫酸盐测序(oxBS-Seq), TET辅助重亚硫酸盐测序(TAB-Seq)

    RNA修饰: RNA甲基化免疫共沉淀技术(MeRIP)

    蛋白质与核酸相互作用: RIP-Seq, ChIP-Seq, CLIP-Seq

    还有最近比较火的 ATAC-Seq ATAC-seq 能干啥?(http://www.biotrainee.com/thread-1218-1-1.html)

    1 ChIP-Seq技术

    1.1 概念

    染色质免疫共沉淀技术(Chromatin Immunoprecipitation,ChIP)也称结合位点分析法,研究体内蛋白质与DNA相互作用的一种方法,通常用于转录因子结合位点或组蛋白特异性修饰位点的研究。
    ChIP第二代测序技术相结合的ChIP-seq技术,能高效的在全基因组范围内检测与组蛋白、转录因子等互作的DNA片段。

    1.2 ChIP-seq技术原理

    在生理状态下,把细胞内的DNA与蛋白质交联(Crosslink)后裂解细胞,分离染色体,通过超声或酶处理将染色质随机切割;
    利用抗原抗体的特异性识别反应,将与目的蛋白相结合的DNA片段沉淀下来;
    再通过反交联(Reverse crosslink)释放结合蛋白的DNA片段;
    纯化;
    测序获得DNA片段的序列,最后将这些DNA片段比对到对应的参考基因组上。
    这里写图片描述

    2 ChIP-Seq数据分析

    2.1 数据下载

    GSE98149 (包含H3K9me3的全部阶段,H3K4me3和H3K27me3的zygote、E6.5 Epi、E6.5 Exe、E7.5 Epi、E7.5 Exe、E8.5 embryo、Esc)
    Title:Reprogramming of H3K9me3-dependent heterochromatin during mammalian early embryo development [ChIP-seq]
    Organism:Mus musculus

    for ((i=594;i<=670;i++));
    do
    wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP105/SRP105176/SRR5479$i/SRR5479$i.sra;
    done &
    # sratookit: .sra 文件 → fastq文件
    ls *sra |while read id;
    do
    /home/chen/sratoolkit.2.8.2-ubuntu64/bin/fastq-dump --gzip --split-3 $id;
    done &
    
    # 下载小鼠参考基因组的 index
    wget -c "ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip" &
    
    # 解压
    unzip mm10.zip &

    2.2 质量控制(data_assess)

    # Fastqc 进行质控
    ls *fq | while read id; do fastqc -t 4 $id; done &
    
    # multiqc:质控结果批量查看
    multiqc *fastqc.zip --export &
    ## trimmomatic 
    
    # 安装 trimmomatic
    wget -c http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.38.zip &
    unzip Trimmomatic-0.38.zip
    
    # 数据清理
    # -threads 设置多线程运行
    java -jar "/data/chen/biosoft/Trimmomatic-0.38/trimmomatic-0.38.jar" PE -threads 2 -phred33 \
    
    # 2个输入文件
    ${name}_1.fq.gz ${name}_2.trim.fq.gz \
    
    # 4个输出文件
    ${name}_R1.clean.fq.gz ${name}_R1.unpaired.fq.gz \
    ${name}_R2.clean.fq.gz ${name}_R2.unpaired.fq.gz \
    
    # ILLUMINACLIP:去接头
    # "$adapter"/Exome.fa :adapter 序列的 fasta 文件
    # 2:16 个碱基长度的种子序列中可以有 2 个错配
    # 30:采用回文模式时匹配得分至少为30 (约50个碱基)
    # 10:采用简单模式时匹配得分至少为10 (约17个碱基)
    ILLUMINACLIP:"$adapter"/Exome.fa:2:30:10 \
    
    # LEADING:3,从序列的开头开始去掉质量值小于 3 的碱基;
    # TRAILING:3,从序列的末尾开始去掉质量值小于 3 的碱基;
    # SLIDINGWINDOW:4:15,从 5' 端开始以 4 bp 的窗口计算碱基平均质量,
    # 如果此平均值低于 15,则从这个位置截断 read;
    # HEADCROP:<length> 在reads的首端切除指定的长度;
    # MINLEN:36, 如果 reads 长度小于 36 bp 则扔掉整条 read。
    LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 HEADCROP:10 MINLEN:36

    这里写图片描述
    这里写图片描述

    2.3 比对到参考基因组(mapping_analysis)

    Bowtie2 或 BWA

    # bowtie2 [options]* -x <bt2-idx> {-1 <m1> -2 <m2> | -U <r>} -S [<hit>]
    
    # -p/--threads NTHREADS 设置线程数. Default: 1
    # -q reads 是 fastq 格式的
    # -x <bt2-idx> index 路径
    # -1 <m1> 双末端测序的 _1.fastq 路径。可以为多个文件,并用逗号分开;多个文件必须和 -2 <m2> 中制定的文件一一对应。
    # -2 <m2> 双末端测序的 _2.fastq 路径.
    # -U <r> 非双末端测序的 fastq 路径。可以为多个文件,并用逗号分开。
    # -S <hit> 输出 Sam 格式文件。
    # -3/--trim3 <int> 剪掉3'端<int>长度的碱基,再用于比对。(default: 0).
    # 用fastqc看了看数据质量,发现3端质量有点问题,我就用了-3 5 --local参数,
    # --local 如果fq文件是没有经过 trim 的,可以用局部比对执行 soft-clipping,加上参数--local 。该模式下对read进行局部比对, 从而, read 两端的一些碱基不比对,从而使比对得分满足要求.
    
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479594_1.fastq -2 SRR5479594_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/MII_rep1.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479595_1.fastq -2 SRR5479595_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/MII_rep2.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479597_1.fastq -2 SRR5479597_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/sperm_rep1.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479598_1.fastq -2 SRR5479598_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/sperm_rep2.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479596_1.fastq -2 SRR5479596_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/sperm_input.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479599_1.fastq -2 SRR5479599_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/zygote_input.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479605_1.fastq -2 SRR5479605_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/zygote_rep1.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479606_1.fastq -2 SRR5479606_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/zygote_rep2.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479607_1.fastq -2 SRR5479607_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/2cell_rep1.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479608_1.fastq -2 SRR5479608_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/2cell_rep2.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479609_1.fastq -2 SRR5479609_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/4cell_rep1.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479610_1.fastq -2 SRR5479610_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/4cell_rep2.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479611_1.fastq -2 SRR5479611_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/8cell_rep1.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479612_1.fastq -2 SRR5479612_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/8cell_rep2.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479613_1.fastq -2 SRR5479613_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/morula_rep1.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479614_1.fastq -2 SRR5479614_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/morula_rep2.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479615_1.fastq -2 SRR5479615_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/ICM_rep1.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479616_1.fastq -2 SRR5479616_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/ICM_rep2.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479617_1.fastq -2 SRR5479617_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/TE_rep1.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479618_1.fastq -2 SRR5479618_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/TE_rep2.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479623_1.fastq -2 SRR5479623_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/ESC_rep1.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479624_1.fastq -2 SRR5479624_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/ESC_rep2.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479625_1.fastq -2 SRR5479625_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/TSC_rep1.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479626_1.fastq -2 SRR5479626_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/TSC_rep2.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479627_1.fastq -2 SRR5479627_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/TSC_rep3.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479628_1.fastq -2 SRR5479628_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E65Epi_input.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479634_1.fastq -2 SRR5479634_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E65Epi_rep1.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479635_1.fastq -2 SRR5479635_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E65Epi_rep2.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479636_1.fastq -2 SRR5479636_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E65Exe_input.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479642_1.fastq -2 SRR5479642_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E65Exe_rep1.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479643_1.fastq -2 SRR5479643_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E65Exe_rep2.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479644_1.fastq -2 SRR5479644_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E75Epi_input.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479650_1.fastq -2 SRR5479650_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E75Epi_rep1.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479651_1.fastq -2 SRR5479651_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E75Epi_rep2.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479652_1.fastq -2 SRR5479652_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E75Epi_rep3.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479653_1.fastq -2 SRR5479653_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E75Epi_rep4.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479654_1.fastq -2 SRR5479654_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E75Exe_input.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479661_1.fastq -2 SRR5479661_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E75Exe_rep1.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479662_1.fastq -2 SRR5479662_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E75Exe_rep2.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479663_1.fastq -2 SRR5479663_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E85Epi_input.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479669_1.fastq -2 SRR5479669_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E85Epi_rep1.sam &
    bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479670_1.fastq -2 SRR5479670_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E85Epi_rep2.sam &

    这里写图片描述
    这里写图片描述
    这里写图片描述
    这里写图片描述

    2.4 搜峰(Peak_calling)

    MACS2

    peaks calling:寻找可能的结合位点,即基因组中大量reads富集的区域。

    2.4.1 MACS2 核心: callpeak 用法

       # Example for regular peak calling
        macs2 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01
       # Example for broad peak calling
       macs2 callpeak -t ChIP.bam -c Control.bam --broad -g hs --broad-cutoff 0.1
    # 批量callpeak
    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 IgGold.bam -t cbx7.bam -q 0.05 -f BAM -g mm -n cbx7 2> cbx7.macs2.log &
    macs2 callpeak -c IgGold.bam -t RYBP.bam -q 0.01 -f BAM -g mm -n RYBP 2>RYBP.macs2.log &

    -t/–treatment FILENAME——处理组输入
    This is the only REQUIRED parameter for MACS. File can be in any supported format specified by –format option. Check –format for detail. If you have more than one alignment files, you can specify them as -t A B C. MACS will pool up all these files together.

    -c/–control——对照组输入
    The control or mock data file. Please follow the same direction as for -t/–treatment.

    -f/–format FORMAT——-t和-c提供文件的格式,目前MACS能够识别的格式有 “ELAND”, “BED”, “ELANDMULTI”, “ELANDEXPORT”, “ELANDMULTIPET” (双端测序), “SAM”, “BAM”, “BOWTIE”, “BAMPE”, “BEDPE”. 除”BAMPE”, “BEDPE”需要特别声明外,其他格式都可以用 AUTO自动检测。如果不提供这项,就是自动检测选择。
    Format of tag file, can be “ELAND”, “BED”, “ELANDMULTI”, “ELANDEXPORT”, “ELANDMULTIPET” (for pair-end tags), “SAM”, “BAM”, “BOWTIE”, “BAMPE” or “BEDPE”. Default is “AUTO” which will allow MACS to decide the format automatically. “AUTO” is also usefule when you combine different formats of files. Note that MACS can’t detect “BAMPE” or “BEDPE” format with “AUTO”, and you have to implicitly specify the format for “BAMPE” and “BEDPE”.

    -g/–gsize——基因组大小,默认提供了hs, mm, ce, dm选项,不在其中的话,比如说拟南芥,就需要自己提供了(拟南芥根据NCBI显示是119,667,750,也就是1.2e8)。
    PLEASE assign this parameter to fit your needs!
    It’s the mappable genome size or effective genome size which is defined as the genome size which can be sequenced. Because of the repetitive features on the chromsomes, the actual mappable genome size will be smaller than the original size, about 90% or 70% of the genome size. The default hs – 2.7e9 is recommended for UCSC human hg18 assembly. Here are all precompiled parameters for effective genome size:
    hs: 2.7e9 (人类是2.7e9,也就是2.7G)
    mm: 1.87e9
    ce: 9e7
    dm: 1.2e8

    -n/–name——输出文件的前缀名。表示实验的名字, 请取一个有意义的名字。
    The name string of the experiment. MACS will use this string NAME to create output files like ‘NAME_peaks.xls’, ‘NAME_negative_peaks.xls’, ‘NAME_peaks.bed’ , ‘NAME_summits.bed’, ‘NAME_model.r’ and so on. So please avoid any confliction between these filenames and your existing files.

    -B/–bdg
    会保存更多的信息在bedGraph文件中,如fragment pileup, control lambda, -log10pvalue and -log10qvalue scores。
    If this flag is on, MACS will store the fragment pileup, control lambda, -log10pvalue and -log10qvalue scores in bedGraph files. The bedGraph files will be stored in current directory named NAME+’_treat_pileup.bdg’ for treatment data, NAME+’_control_lambda.bdg’ for local lambda values from control, NAME+’_treat_pvalue.bdg’ for Poisson pvalue scores (in -log10(pvalue) form), and NAME+’_treat_qvalue.bdg’ for q-value scores from Benjamini–Hochberg–Yekutieli procedure http://en.wikipedia.org/wiki/False_discovery_rate#Dependent_tests

    -q/–qvalue——q值,也就是最小的PDR阈值, 默认是0.05。q值是根据p值利用BH计算,也就是多重试验矫正后的结果。
    The qvalue (minimum FDR) cutoff to call significant regions. Default is 0.01. For broad marks, you can try 0.05 as cutoff. Q-values are calculated from p-values using Benjamini-Hochberg procedure.

    -p/–pvalue——p值,指定 p 值后 MACS2 就不会用 q 值了。
    The pvalue cutoff. If -p is specified, MACS2 will use pvalue instead of qvalue.

    -m/–mfold——和MFOLD有关,而MFOLD和MACS预构建模型有关,默认是5:50,MACS会先寻找100多个peak区构建模型,一般不用改,因为你不懂。
    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:5,50 means using all regions not too low (>5) and not too high (<50) to build paired-peaks model. If MACS can not find more than 100 regions to build model, it will use the –extsize parameter to continue the peak detection ONLY if –fix-bimodal is set.

    2.4.2 callpeak 结果文件说明

    # (在当前目录下)统计 *bed 的行数(peak数)
    wc -l *bed
      2384 cbx7_summits.bed
      8342 ring1B_summits.bed
         0 RYBP_summits.bed
      7619 suz12_summits.bed
    # 在文件a中统计 hello 出现的行数:
    # grep hello a | wc -l
    # wc(word count)
    #-c 统计字节数。
    #-l 统计行数。line
    #-m 统计字符数。这个标志不能与 -c 标志一起使用。
    #-w 统计字数。一个字被定义为由空白、跳格或换行字符分隔的字符串。
    #-L 打印最长行的长度。

    callpeak会得到如下结果文件:

    NAME_summits.bed:Browser Extensible Data,记录每个peak的peak summits,话句话说就是记录极值点的位置。MACS建议用该文件寻找结合位点的motif。能够直接载入UCSC browser,用其他软件分析时需要去掉第一行。

    NAME_peaks.xls:以表格形式存放peak信息,虽然后缀是xls,但其实能用文本编辑器打开,和bed格式类似,但是以1为基,而bed文件是以0为基.也就是说xls的坐标都要减一才是bed文件的坐标。

    NAME_peaks.narrowPeakNAME_peaks.broadPeak 类似。后面4列表示为, integer score for display, fold-change,-log10pvalue,-log10qvalue,relative summit position to peak start。内容和NAME_peaks.xls基本一致,适合用于导入R进行分析。

    NAME_model.r:能通过$ Rscript NAME_model.r作图,得到是基于你提供数据的peak模型。

    .bdg:能够用 UCSC genome browser 转换成更小的 bigWig 文件。

    2.4.3 bdg file → wig file

    为了方便在IGV上查看ChIP-seq的结果和后期的可视化展示,需要把macs2的结果转化为bw提供给IGV。一共分为三步:

    第一步:使用 bdgcmp 得到 FE 或者 logLR 转化后的文件 (Run MACS2 bdgcmp to generate fold-enrichment and logLR track)

    macs2 bdgcmp -t H3K36me1_EE_rep1_treat_pileup.bdg -c H3K36me1_EE_rep1_control_lambda.bdg -o H3K36me1_EE_rep1_FE.bdg -m FE
    
    macs2 bdgcmp -t  H3K36me1_EE_rep1_treat_pileup.bdg -c H3K36me1_EE_rep1_control_lambda.bdg -o H3K36me1_EE_rep1_logLR.bdg  -m logLR -p 0.00001
    
    # 参数解释
    # -m FE 计算富集倍数降低噪音
    # -p 为了避免log的时候input值为0时发生error,给予一个很小的值

    第二步:预处理文件,下载对应参考基因组染色体长度文件

    使用conda安装以下三个处理软件:
    bedGraphToBigWig
    bedClip
    bedtools

    下载染色体长度文件:http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/chromInfo.txt.gz 并解压(针对human,其余物种的皆可以按照类似网址下载)

    写一个小小的sh脚本方便一步转化name.sh:

    #!/bin/bash
    # check commands: slopBed, bedGraphToBigWig and bedClip
    which bedtools &>/dev/null || { echo "bedtools not found! Download bedTools: <http://code.google.com/p/bedtools/>"; exit 1; }
    which bedGraphToBigWig &>/dev/null || { echo "bedGraphToBigWig not found! Download: <http://hgdownload.cse.ucsc.edu/admin/exe/>"; exit 1; }
    which bedClip &>/dev/null || { echo "bedClip not found! Download: <http://hgdownload.cse.ucsc.edu/admin/exe/>"; exit 1; }
    # end of checking
    if [ $# -lt 2 ];then
    echo "Need 2 parameters! <bedgraph> <chrom info>"
    exit
    fi
    F=$1
    G=$2
    bedtools slop -i ${F} -g ${G} -b 0 | bedClip stdin ${G} ${F}.clip
    LC_COLLATE=C sort -k1,1 -k2,2n ${F}.clip > ${F}.sort.clip
    bedGraphToBigWig ${F}.sort.clip ${G} ${F/bdg/bw}
    rm -f ${F}.clip ${F}.sort.clip

    chmod +x name.sh

    第三步:生成 bw 文件

    ./name.sh H3K36me1_EE_rep1_FE.bdg hg19.len
    ./name.sh H3K36me1_EE_rep1_logLR.bdg hg19.len

    最后得到产物,至于的使用哪一个作为输入文件大家就根据需要来吧

    H3K36me1EErep1_FE.bw
    H3K36me1EErep1_logLR.bw

    2.5 峰注释(Peak_anno)

    ChIPseeker

    ChIPseeker的功能分为三类:
    注释:提取peak附近最近的基因,注释peak所在区域。
    比较:估计ChIP peak数据集中重叠部分的显著性;整合GEO数据集,以便于将当前结果和已知结果比较。
    可视化:peak的覆盖情况;TSS区域结合的peak的平均表达谱和热图;基因组注释;TSS距离;peak和基因的重叠。

    # 加载ChIPseeker、基因组注释包和bed数据
    biocLite("ChIPseeker")
    biocLite("org.Mm.eg.db")
    biocLite("TxDb.Mmusculus.UCSC.mm10.knownGene")
    library("ChIPseeker")
    
    # 下载
    source ("https://bioconductor.org/biocLite.R")
    biocLite("ChIPseeker")
    biocLite("org.Mm.eg.db")
    biocLite("TxDb.Mmusculus.UCSC.mm10.knownGene")
    biocLite("clusterProfiler")
    biocLite("ReactomePA")
    biocLite("DOSE")
    
    #载入
    library("ChIPseeker")
    library("org.Mm.eg.db")
    library("TxDb.Mmusculus.UCSC.mm10.knownGene")
    txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
    library("clusterProfiler")
    
    # 读入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")

    Chip peaks coverage plot

    查看peak在全基因组的位置

    covplot(ring1B,chrs=c(“chr17”, “chr18”)) #specific chr

    ring1B

    suz12

    cbx7

    RYBP

    ring1B(chr17&18)

    ● Heatmap of ChIP binding to TSS regions

    promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
    tagMatrix <- getTagMatrix(ring1B, windows=promoter)
    tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color=”red”)

    Average Profile of ChIP peaks binding to TSS region

    ● Confidence interval estimated by bootstrap method

    plotAvgProf(tagMatrix, xlim=c(-3000, 3000), conf = 0.95, resample = 1000)

    peak的注释

    peak的注释用annotatePeak**,TSS (transcription start site) region 可以自己设定,默认是(-3000,3000),TxDb 是指某个物种的基因组,例如TxDb.Hsapiens.UCSC.hg38.knownGene, TxDb.Hsapiens.UCSC.hg19.knownGene for human genome hg38 and hg19, TxDb.Mmusculus.UCSC.mm10.knownGene and TxDb.Mmusculus.UCSC.mm9.knownGene for mouse mm10 and mm9.

    peakAnno <- annotatePeak(ring1B, tssRegion=c(-3000, 3000),
    TxDb=txdb, annoDb=”org.Mm.eg.db”)

    可视化 Pie and Bar plot

    plotAnnoBar(peakAnno)
    vennpie(peakAnno)
    upsetplot(peakAnno)

    饼图:

    条形图:

    upsetplot: upset技术适用于多于5个集合的表示情况。

    可视化TSS区域的TF binding loci

    plotDistToTSS(peakAnno,
    title=”Distribution of transcription factor-binding loci\nrelative to TSS”)

    多个peak的比较

    多个peak set注释时,先构建list,然后用lapply. list(name1=bed_file1,name2=bed_file2)RYBP的数据有问题,这里加上去,会一直报错。

    peaks <- list(cbx7=cbx7,ring1B=ring1B,suz12=suz12)
    promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
    tagMatrixList <- lapply(peaks, getTagMatrix, windows=promoter)
    plotAvgProf(tagMatrixList, xlim=c(-3000, 3000))
    plotAvgProf(tagMatrixList, xlim=c(-3000, 3000), conf=0.95,resample=500, facet=”row”)
    tagHeatmap(tagMatrixList, xlim=c(-3000, 3000), color=NULL)

    ChIP peak annotation comparision

    peakAnnoList <- lapply(peaks, annotatePeak, TxDb=txdb,
    tssRegion=c(-3000, 3000), verbose=FALSE)
    plotAnnoBar(peakAnnoList)
    plotDistToTSS(peakAnnoList)

    Overlap of peaks and annotated genes

    genes= lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
    vennplot(genes)

    展开全文
  • 基因组学技术的最新发展实现了通过RNA-seq进行转录谱分析和通过ChIP-seq进行蛋白质结合谱分析。 对这两种类型数据的综合分析使我们能够从基因组共定位和下游靶基因表达中研究TF和HM的相互作用。结果:我们提出了一...
  • BMC Bioinformatics 2014, 15:280 http://www.biomedcentral.com/1471-2105/15/280 SOFTWARE Open Access HiChIP: a high-throughput pipeline for integrative analysis of ChIP-Seq data 1,2 1 1 1 1
  • ChIP-seq 原理及实践

    2021-04-19 10:46:14
    包含RNA-seq/ChIP-seq/ATAC-seq/scRNA-seq等原理及实践 这个可能更全面,之后仔细看吧 简书:有数据可以跟着做
    展开全文
  • 随着染色质免疫共沉淀芯片(CHIP-chip)技术的发展,大量的表观基因组分析技术已经出现,如CHIP-seq、DNase I超灵敏位点测序(DNase-seq)、ATAC-seq等。单细胞测序的出现使下一代测序发生了革命性的变化,在单
  • 随着染色质免疫共沉淀芯片(CHIP-chip)技术的发展,大量的表观基因组分析技术已经出现,如CHIP-seq、DNase I超灵敏位点测序(DNase-seq)、ATAC-seq等。单细胞测序的出现使下一代测序发生了革命性的变化,在单
  • 琶音使我们能够有效地比较许多ChIP-seq数据集,这些数据集由从各种细胞,条件和生物体中收集的多种类型的DNA结合蛋白组成。 重要的是,这些相互作用模式广泛反映了结合事件的生物学特性。 为了生成称为琶音轮廓的...
  • Chip-seq流程详解
  • SUPERmerge是一种ChIP-seq读取堆积分析和注释算法,用于以单个碱基对的分辨率水平研究具有较宽染色质域的弥散组蛋白修饰ChIP-seq数据集的比对(BAM)文件。 SUPERmerge允许灵活调节各种读取堆积参数,从而揭示读取岛...
  • ABC --(来自 ChIP-Seq 的等位基因特异性结合) 在 ChIP-Seq 实验的比对读数中的已知杂合位置识别潜在的等位基因特异性结合事件。 ABC 至少需要两 (2) 个文件,一个 ChIP-Seq 实验的排序 BAM/SAM 文件和一个包含杂...
  • ChIP-Seq 常见问题解答

    2020-11-08 21:09:36
    关于ChIP-Seq 染色质免疫共沉淀(ChromatinImmunopreciption,ChIP)是研究蛋白质-DNA相互作用的经典实验方法,ChIP 与高通量测序的结合(ChIP-Seq)可以在全基因组范围内对蛋白质结合位点进行高效而准确地筛选与...
  • CHIP-Seq数据分析流程

    千次阅读 2020-11-10 20:11:44
    文章目录CHIP-Seq数据分析流程相关软件安装数据下载将SRA转化为fastq文件数据质控,过滤数据质控然后进行质控,过滤低质量reads,去接头比对下载小鼠参考基因组的索引和注释文件比对sam文件转bam为bam文件建立索引...
  • 1 ChIP-Seq技术 1.1 概念 1.2 ChIP-seq技术原理 2 ChIP-Seq数据分析 2.1 数据下载 2.2 质量控制(data_assess) 2.3 比对到参考基因组(mapping_analysis) 2.4 搜峰(Peak_calling) MACS2 2.4.1 MACS2 核心:...
  • ChIP-seq是进行体内检测TFBS的主要方法。然而,ChIP-seq依赖于抗体质量,这对低表达的蛋白质具有很大的挑战性。所以,ChIP-seq通常在规模上受到限制,难以进行高通量扩展。因此,只有少数转录因子可以获得结合位点...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 3,005
精华内容 1,202
关键字:

chip-seq