systemPipeR 包中文翻译之 ChIP-seq 分析

作者: 热衷组培的二货潜 | 来源:发表于2019-06-08 21:00 被阅读54次

systemPipeR 包中文翻译之 ChIP-seq 分析
Bioconda软件安装神器：多版本并存、环境复制、环境导出
ATAC-seq 分析（上）
TCseq：基因表达趋势分析
学员分享-Chip-seq 实战分析流程
ATAC-seq分析练习
ChIP-Seq分析之ChIPQC结果含义
ChIP-seq数据比对注意事项
ChIP-seq数据质控与过滤
chip-seq实验方法及数据分析

原文英文链接：ChIP-Seq Workflow Template
作者：Daniela Cassol (danielac@ucr.edu) and Thomas Girke (thomas.girke@ucr.edu)

最近更新时间：2019.05.05

R包

systemPipeR 1.19.0

内容目录（我自己加在前面的）

1 简单介绍
2 工作环境
3 Read 预处理
4 比对
5 coverage 数据的实用程序( Utilities )
6 使用 MACS2 进行 Peak calling
7 用基因组背景注 Peaks
8 交集的 Peaks 的 reads 计数
9 差异结合分析
10 GO富集分析
11 Motif 分析
12 版本信息
13 基金
参考文献

1 简单介绍

用户希望在此提供有关其ChIP-Seq项目设计的背景信息。

1.1 背景和目标

本报告描述了几项ChIP-Seq实验的分析，这些生物实验研究了转录因子的DNA结合模式。

1.2 实验设计

通常，用户希望在此指定与其NGS研究分析相关的所有信息。这包括FASTQ文件，实验设计，参考基因组，基因注释等的详细描述。

2 工作环境

2.1 生成工作流程环境

将包含样本数据的工作流环境加载到当前工作目录中。样本数据描述见此

library(systemPipeRdata)
genWorkenvir(workflow = "chipseq")
setwd("chipseq")

可选择性的，我们也可以在linux环境中通过命令行运行：

Rscript -e "systemPipeRdata::genWorkenvir(workflow='chipseq')"

在 genWorkenvir 函数生成的工作流环境中，所有数据输入都存储在 data/ 目录中，所有分析结果将写入单独的 result/ 目录，而 systemPipeChIPseq.Rmd 脚本和 targets 文件应位于父目录中。R 从此父目录运行。其他参数文件存储在 param / 下。
为了处理真实数据，用户希望类似地组织自己的数据，并将所有测试数据替换为自己的数据。要在新数据上重新运行已建立的工作流，初始 targets 文件以及相应的 FASTQ 文件通常是用户需要提供的唯一输入。

2.2 运行工作流程

现在通过你的 R 中打开 R markdown 脚本 systemPipeChIPseq.Rmd ，例如：vim-r、Rstudio，然后按照如下进行运行。

2.2.1 在计算节点运行 `R`

在 Vim 中打开此工作流的 Rmd 文件并通过 F2（或其他）键附加已连接的 R 后，使用以下命令序列在计算机节点上运行 R 。

q("no") # 关闭节点上运行的R

srun --x11 --partition=short --mem=2gb --cpus-per-task 4 --ntasks 1 --time 2:00:00 --pty bash -l
module load R/3.4.2
R

现在检查您的 R 是否在群集的计算机节点上运行并评估您的环境。

system("hostname")  # 返回以 i 或 c 开头的计算节点的名称
getwd()  # 查看当前 R 的工作路径
dir()  # 返回当前工作路径下的文件

需要加载 systemPipeR 包以执行本报告中显示的分析步骤（H Backman和Girke 2016）。

library(systemPipeR)

如果适用，用户可以加载 systemPipeR 未提供的自定义功能。如果不是这样，请跳过此步骤。

source("systemPipeChIPseq_Fct.R")

3 Read 预处理

3.1 `targets` 文件提供的实验定义

targets 文件定义了分析工作流的所有 FASTQ 文件和样本比较。

targetspath <- system.file("extdata", "targets_chip.txt", package = "systemPipeR")
targets <- read.delim(targetspath, comment.char = "#")
targets[1:4, -c(5, 6)]

                  FileName SampleName Factor SampleLong SampleReference
1 ./data/SRR446027_1.fastq        M1A     M1  Mock.1h.A                
2 ./data/SRR446028_1.fastq        M1B     M1  Mock.1h.B                
3 ./data/SRR446029_1.fastq        A1A     A1   Avr.1h.A             M1A
4 ./data/SRR446030_1.fastq        A1B     A1   Avr.1h.B             M1B

3.2 `Read` 质量筛选与修剪

下面的示例演示如何使用 ShortRead 包提供的实用程序设计自定义读取预处理函数，然后将批处理模式下的 preprocessReads 应用于相应 SYSargs 实例（下面的 args 对象）中引用的所有 FASTQ 样本。systemPipeR 的介绍中提供了有关读取预处理的更多详细信息。

args <- systemArgs(sysma = "param/trim.param", mytargets = "targets_chip.txt")
filterFct <- function(fq, cutoff = 20, Nexceptions = 0) {
    qcount <- rowSums(as(quality(fq), "matrix") <= cutoff, na.rm = TRUE)
    fq[qcount <= Nexceptions]
    # 保留读取 Phred 分数 > = 阈值 cutoff 为 N 的 reads
    # 例外
}
preprocessReads(args = args, Fct = "filterFct(fq, cutoff=20, Nexceptions=0)", 
    batchsize = 1e+05)
writeTargetsout(x = args, file = "targets_chip_trim.txt", overwrite = TRUE)

3.3 `FASTQ` 质量报告

以下 seeFastq 和 seeFastqPlot 函数为一组 FASTQ 文件生成并绘制一系列有用的质量统计数据，包括每个循环质量 box plots，基本比例，基本级质量趋势，相对 k-mer 多样性，reads 的长度和分布，reads 高于质量 cutoff 和平均质量分布。结果将写入名为fastqReport.pdf的PDF文件中。

args <- systemArgs(sysma = "param/tophat.param", mytargets = "targets_chip.txt")
library(BiocParallel)
library(batchtools)
f <- function(x) {
    library(systemPipeR)
    args <- systemArgs(sysma = "param/tophat.param", mytargets = "targets_chip.txt")
    seeFastq(fastq = infile1(args)[x], batchsize = 1e+05, klength = 8)
}
moduleload(modules(args))  # Skip if a module system is not used
resources <- list(walltime = 120, ntasks = 1, ncpus = cores(args), 
    memory = 1024)
param <- BatchtoolsParam(workers = 4, cluster = "slurm", template = "batchtools.slurm.tmpl", 
    resources = resources)
fqlist <- bplapply(seq(along = args), f, BPPARAM = param)

pdf("./results/fastqReport.pdf", height = 18, width = 4 * length(fqlist))
seeFastqPlot(unlist(fqlist, recursive = FALSE))
dev.off()

Figure 1：18 个样本的 `FASTQ` 质量报告

4 比对

4.1 使用 `Bowtie2` 进行 `Reads` 的比对

通过 Bowtie2 将 reads 比对到参考基因组序列（Langmead, Salzberg 2012）。比对的参数设置在 bowtieSE.param 文件中定义。在 ChIP-Seq 实验中，通常过滤比对到多个位置的 reads 更为恰当。为实现此目的，用户希望在 bowtieSE.param 文件中删除参数设置 -k 50 non-deterministic 。

rgs <- systemArgs(sysma = "param/bowtieSE.param", mytargets = "targets_chip_trim.txt")

sysargs(args)[1]  # 第一个 FASTQ 文件的命令行参数

moduleload(modules(args))  # 如果不需要 module 则可以跳过

system("bowtie2-build ./data/tair10.fasta ./data/tair10.fasta")
# 对基因组文件建立索引

resources <- list(walltime = 120, ntasks = 1, ncpus = cores(args), 
    memory = 1024)

reg <- clusterRun(args, conffile = ".batchtools.conf.R", Njobs = 18, 
    template = "batchtools.slurm.tmpl", runid = "01", resourceList = resources)

getStatus(reg = reg)

waitForJobs(reg = reg)

writeTargetsout(x = args, file = "targets_bam.txt", overwrite = TRUE)

或者，可以在单个系统上顺序运行比对。

runCommandline(args)

检查是否已创建所有 BAM 文件

file.exists(outpaths(args))

4.2 `Reads` 和比对信息统计

以下概述了每个样品中的 reads，以及它们中有多少比对上了参考基因组。

read_statsDF <- alignStats(args = args)

write.table(read_statsDF, "results/alignStats.xls", row.names = FALSE, 
    quote = FALSE, sep = "\t")

read.delim("results/alignStats.xls")

4.3 创建用于在 `IGV` 中查看 `BAM` 文件的链接

symLink2bam 函数创建链接以查看基因组浏览器（如 IGV ）中的 BAM 比对文件，而无需将这些大文件移动到本地系统。相应的URL将写入具有在 urlfile 下指定的路径的文件，此处为 IGVurl.txt 。

symLink2bam(sysargs = args, htmldir = c("~/.html/", "somedir/"), 
    urlbase = "http://biocluster.ucr.edu/~tgirke/", urlfile = "./results/IGVurl.txt")

5 coverage 数据的实用程序( Utilities )

以下介绍几种对 ChIP-Seq 数据有用的程序。它们不是实际工作流程的一部分。

5.1 `Rle` 属性中储存的 `coverage` 信息

library(rtracklayer)
library(GenomicRanges)
library(Rsamtools)
library(GenomicAlignments)

aligns <- readGAlignments(outpaths(args)[1])
cov <- coverage(aligns)
cov

5.2 调整已比对上的 `reads`

trim(resize(as(aligns, "GRanges"), width = 200))

5.3 `Naive Peak calling`

islands <- slice(cov, lower = 15)
islands[[1]]

5.4 绘制定义区间的 `coverage` 信息

library(ggbio)
myloc <- c("Chr1", 1, 1e+05)
ga <- readGAlignments(outpaths(args)[1], use.names = TRUE, param = ScanBamParam(which = GRanges(myloc[1], 
    IRanges(as.numeric(myloc[2]), as.numeric(myloc[3])))))
autoplot(ga, aes(color = strand, fill = strand), facets = strand ~ 
    seqnames, stat = "coverage")

6 使用 MACS2 进行 Peak calling

6.1 在进行 `Peak calling` 之前将每个样本的重复文件 `BAM` 合并成一个 `BAM` 文件。

合并技术和/或 生物重复的 BAM 文件可以通过增加 reads覆盖度来提高 Peak calling 的灵敏度。mergeBamByFactor 函数根据 Factor 指定的分组信息合并 BAM 文件，此处为导入目标文件的 Factor 列。它还返回一个更新的 SYSargs 对象，该对象包含合并的 BAM 文件的路径以及没有重复的任何未合并的文件。如果不需要合并 BAM 文件，则可以跳过此步骤。

args <- systemArgs(sysma = NULL, mytargets = "targets_bam.txt")

args_merge <- mergeBamByFactor(args, overwrite = TRUE)

writeTargetsout(x = args_merge, file = "targets_mergeBamByFactor.txt", 
    overwrite = TRUE)

6.2 没有`input/reference`的`Peak calling`

MACS2 可以在有和没有 input 样本的情况下对ChIP-Seq 数据进行 Peak calling（Zhang et al.2008）。以下进行在args对象中指定的所有样本没有input的情况下 Peak calling 。注意，由于样本数据量小，需要在此处使用 nomodel 设置运行 MACS2 。对于实际数据集，用户希望在相应的 *.param 文件中删除此参数。

args <- systemArgs(sysma = "param/macs2_noinput.param", mytargets = "targets_mergeBamByFactor.txt")

sysargs(args)[1]  # 第一个 FASTQ 文件的命令行参数

runCommandline(args)
file.exists(outpaths(args))

writeTargetsout(x = args, file = "targets_macs.txt", overwrite = TRUE)

6.3 有`input/reference`的`Peak calling`

要使用有 input 的样本进行 Peak calling ，可以在初始 targets 文件的 SampleReference 列中最方便地指定它们。writeTargetsRef 函数使用此信息创建一个目标文件，用于运行带有相应 input 样本的 MACS2 。

writeTargetsRef(infile = "targets_mergeBamByFactor.txt", outfile = "targets_bam_ref.txt", 
    silent = FALSE, overwrite = TRUE)

args_input <- systemArgs(sysma = "param/macs2.param", mytargets = "targets_bam_ref.txt")

sysargs(args_input)[1]  # 第一个 FASTQ 文件的命令行参数

runCommandline(args_input)
file.exists(outpaths(args_input))

writeTargetsout(x = args_input, file = "targets_macs_input.txt", 
    overwrite = TRUE)

MACS2 的 Peak calling 结果是针对每个样本编写的，以分离结果目录中的文件。它们以具有 MACS2 使用的扩展名的相应文件命名。

6.4 鉴定共同的 `Peaks`

以下示例显示了如何使用 subsetByOverlaps 或olRanges 函数分别在共享最小绝对重叠和/或最小相对重叠的两个 Peak 的共同 Peak 。注意，后者是通过 source 它而导入的自定义函数。

source("http://faculty.ucr.edu/~tgirke/Documents/R_BioCond/My_R_Scripts/rangeoverlapper.R")

peak_M1A <- outpaths(args)["M1A"]

peak_M1A <- as(read.delim(peak_M1A, comment = "#")[, 1:3], "GRanges")

peak_A1A <- outpaths(args)["A1A"]

peak_A1A <- as(read.delim(peak_A1A, comment = "#")[, 1:3], "GRanges")

(myol1 <- subsetByOverlaps(peak_M1A, peak_A1A, minoverlap = 1))
# 返回任何有交集的区域

myol2 <- olRanges(query = peak_M1A, subject = peak_A1A, output = "gr")
# 返回 0L 长度的区域

myol2[values(myol2)["OLpercQ"][, 1] >= 50]
# 返回最小交集 50% 的区域

7 用基因组背景注 Peaks

7.1 使用 `ChIPpeakAnno` 包进行注释

以下分别使用 ChIPpeakAnno 和 ChIPseeker 包对带有基因组背景信息的鉴定到的 Peaks 进行注释（Zhu et al. 2010; Yu，Wang; He 2015）。

library(ChIPpeakAnno)
library(GenomicFeatures)

args <- systemArgs(sysma = "param/annotate_peaks.param", mytargets = "targets_macs.txt")

# txdb <- loadDb('./data/tair10.sqlite')
txdb <- makeTxDbFromGFF(file = "data/tair10.gff", format = "gff", 
    dataSource = "TAIR", organism = "Arabidopsis thaliana")

ge <- genes(txdb, columns = c("tx_name", "gene_id", "tx_type"))

for (i in seq(along = args)) {
    peaksGR <- as(read.delim(infile1(args)[i], comment = "#"), 
        "GRanges")
    annotatedPeak <- annotatePeakInBatch(peaksGR, AnnotationData = genes(txdb))
    df <- data.frame(as.data.frame(annotatedPeak), as.data.frame(values(ge[values(annotatedPeak)$feature, 
        ])))
    write.table(df, outpaths(args[i]), quote = FALSE, row.names = FALSE, 
        sep = "\t")
}

writeTargetsout(x = args, file = "targets_peakanno.txt", overwrite = TRUE)

为每个 Peak 集写入 Peak 注释结果，以将结果目录中的文件分开。它们以 annotate_peaks.param 文件中指定扩展名的相应 Peak 文件命名，此处为 * .peaks.annotated.xls 。

7.2 使用 `ChIPseeker` 包来注释

除了换成 ChIPseeker 来进行注释以外其他步骤都一样

library(ChIPseeker)

for (i in seq(along = args)) {
    peakAnno <- annotatePeak(infile1(args)[i], TxDb = txdb, verbose = FALSE)
    df <- as.data.frame(peakAnno)
    write.table(df, outpaths(args[i]), quote = FALSE, row.names = FALSE, 
        sep = "\t")
}

writeTargetsout(x = args, file = "targets_peakanno.txt", overwrite = TRUE)

ChIPseeker 软件包提供总结性的图。此处采用一个样本来进行展示。

peak <- readPeakFile(infile1(args)[1])

covplot(peak, weightCol = "X.log10.pvalue.")
peakHeatmap(outpaths(args)[1], TxDb = txdb, upstream = 1000, 
    downstream = 1000, color = "red")

plotAvgProf2(outpaths(args)[1], TxDb = txdb, upstream = 1000, 
    downstream = 1000, xlab = "Genomic Region (5'->3')", ylab = "Read Count Frequency")

8 交集的 Peaks 的 reads 计数

countRangeset 函数是一个方便的包装器，可以在多个范围集（这里是 Peaks 范围集）上迭代执行读取计数。在内部，使用 GenomicAlignments 包中的 summarizeOverlaps 函数执行读取计数。生成的计数表直接保存到文件中，每个峰值集一个。

library(GenomicRanges)

args <- systemArgs(sysma = "param/count_rangesets.param", mytargets = "targets_macs.txt")

args_bam <- systemArgs(sysma = NULL, mytargets = "targets_bam.txt")

bfl <- BamFileList(outpaths(args_bam), yieldSize = 50000, index = character())

countDFnames <- countRangeset(bfl, args, mode = "Union", ignore.strand = TRUE)

writeTargetsout(x = args, file = "targets_countDF.txt", overwrite = TRUE)

9 差异结合分析

runDiff 函数使用 edgeR 或 DESeq2 以批处理模式对多个计数表执行差异结合分析（ Robinson，McCarthy, Smyth 2010; Love，Huber, Anders 2014）。在内部，它调用函数 run_edgeR 和 run_DESeq2 。它还使用 dbrfilter 参数下提供的 Fold change 和 FDR 阈值 cutoff 返回下游 filterDEGs 函数的过滤结果和图。

args_diff <- systemArgs(sysma = "param/rundiff.param", mytargets = "targets_countDF.txt")

cmp <- readComp(file = args_bam, format = "matrix")

dbrlist <- runDiff(args = args_diff, diffFct = run_edgeR, targets = targetsin(args_bam), 
    cmp = cmp[[1]], independent = TRUE, dbrfilter = c(Fold = 2, 
        FDR = 1))

writeTargetsout(x = args_diff, file = "targets_rundiff.txt", 
    overwrite = TRUE)

10 GO富集分析

以下对注释的 Peak 集进行 GO 富集分析。

args <- systemArgs(sysma = "param/macs2.param", mytargets = "targets_bam_ref.txt")

args_anno <- systemArgs(sysma = "param/annotate_peaks.param", 
    mytargets = "targets_macs.txt")

annofiles <- outpaths(args_anno)

gene_ids <- sapply(names(annofiles), function(x) unique(as.character(read.delim(annofiles[x])[, 
    "gene_id"])), simplify = FALSE)

load("data/GO/catdb.RData")

BatchResult <- GOCluster_Report(catdb = catdb, setlist = gene_ids, 
    method = "all", id_type = "gene", CLSZ = 2, cutoff = 0.9, 
    gocats = c("MF", "BP", "CC"), recordSpecGO = NULL)

11 Motif 分析

11.1 从基因组中获得 `Peaks` 区间的DNA序列

已知DNA结合 motif 的富集分析或新基因的重新发现需要所鉴定的 Peak 区域的 DNA 序列。为了解析来自参考基因组的相应序列，可以使用来自 Biostrings 包的 getSeq 函数。以下示例解析每个 Peak 集的序列，并将结果保存到单独的 FASTA 文件中，每个 Peak 集对应一个。此外，FASTA 文件中的序列如某些 Motif 发现工具（例如 BCRANK ）所预期的那样通过按 p 值递增来排序。

library(Biostrings)
library(seqLogo)
library(BCRANK)

args <- systemArgs(sysma = "param/annotate_peaks.param", mytargets = "targets_macs.txt")

rangefiles <- infile1(args)

for (i in seq(along = rangefiles)) {
    df <- read.delim(rangefiles[i], comment = "#")
    peaks <- as(df, "GRanges")

    names(peaks) <- paste0(as.character(seqnames(peaks)), "_", 
        start(peaks), "-", end(peaks))

    peaks <- peaks[order(values(peaks)$X.log10.pvalue., decreasing = TRUE)]

    pseq <- getSeq(FaFile("./data/tair10.fasta"), peaks)
    names(pseq) <- names(peaks)
    writeXStringSet(pseq, paste0(rangefiles[i], ".fasta"))
}

11.2 使用 `BCRANK` 发现 `Motif`

Bioconductor 软件包 BCRANK 是 ChIP-Seq 实验 Peak 区域中从头发现 DNA bingding motif 的众多工具之一。给出的示例将此方法应用于第一个 Peak 样本集，并绘制最显著 motif 的图。

set.seed(0)

BCRANKout <- bcrank(paste0(rangefiles[1], ".fasta"), restarts = 25, 
    use.P1 = TRUE, use.P2 = TRUE)

toptable(BCRANKout)

topMotif <- toptable(BCRANKout, 1)

weightMatrix <- pwm(topMotif, normalize = FALSE)

weightMatrixNormalized <- pwm(topMotif, normalize = TRUE)

pdf("results/seqlogo.pdf")
seqLogo(weightMatrixNormalized)
dev.off()

Figure 2: 被 BCRANK 鉴定的其中一个 motif

12 版本信息

sessionInfo()

## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.2 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.10-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.10-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices
## [6] utils     datasets  methods   base     
## 
## other attached packages:
##  [1] systemPipeR_1.19.0          ShortRead_1.43.0           
##  [3] GenomicAlignments_1.21.0    SummarizedExperiment_1.15.0
##  [5] DelayedArray_0.11.0         matrixStats_0.54.0         
##  [7] Biobase_2.45.0              BiocParallel_1.19.0        
##  [9] Rsamtools_2.1.0             Biostrings_2.53.0          
## [11] XVector_0.25.0              GenomicRanges_1.37.0       
## [13] GenomeInfoDb_1.21.0         IRanges_2.19.0             
## [15] S4Vectors_0.23.0            BiocGenerics_0.31.0        
## [17] BiocStyle_2.13.0           
## 
## loaded via a namespace (and not attached):
##  [1] Category_2.51.0          bitops_1.0-6            
##  [3] bit64_0.9-7              RColorBrewer_1.1-2      
##  [5] progress_1.2.0           httr_1.4.0              
##  [7] Rgraphviz_2.29.0         backports_1.1.4         
##  [9] tools_3.6.0              R6_2.4.0                
## [11] DBI_1.0.0                lazyeval_0.2.2          
## [13] colorspace_1.4-1         withr_2.1.2             
## [15] tidyselect_0.2.5         prettyunits_1.0.2       
## [17] bit_1.1-14               compiler_3.6.0          
## [19] graph_1.63.0             formatR_1.6             
## [21] rtracklayer_1.45.0       bookdown_0.9            
## [23] scales_1.0.0             checkmate_1.9.3         
## [25] genefilter_1.67.0        RBGL_1.61.0             
## [27] rappdirs_0.3.1           stringr_1.4.0           
## [29] digest_0.6.18            rmarkdown_1.12          
## [31] AnnotationForge_1.27.0   pkgconfig_2.0.2         
## [33] htmltools_0.3.6          BSgenome_1.53.0         
## [35] limma_3.41.0             rlang_0.3.4             
## [37] RSQLite_2.1.1            GOstats_2.51.0          
## [39] hwriter_1.3.2            dplyr_0.8.0.1           
## [41] VariantAnnotation_1.31.0 RCurl_1.95-4.12         
## [43] magrittr_1.5             GO.db_3.8.2             
## [45] GenomeInfoDbData_1.2.1   Matrix_1.2-17           
## [47] Rcpp_1.0.1               munsell_0.5.0           
## [49] stringi_1.4.3            yaml_2.2.0              
## [51] edgeR_3.27.0             zlibbioc_1.31.0         
## [53] plyr_1.8.4               grid_3.6.0              
## [55] blob_1.1.1               crayon_1.3.4            
## [57] lattice_0.20-38          splines_3.6.0           
## [59] GenomicFeatures_1.37.0   annotate_1.63.0         
## [61] hms_0.4.2                batchtools_0.9.11       
## [63] locfit_1.5-9.1           knitr_1.22              
## [65] pillar_1.3.1             rjson_0.2.20            
## [67] base64url_1.4            codetools_0.2-16        
## [69] biomaRt_2.41.0           XML_3.98-1.19           
## [71] glue_1.3.1               evaluate_0.13           
## [73] latticeExtra_0.6-28      data.table_1.12.2       
## [75] BiocManager_1.30.4       gtable_0.3.0            
## [77] purrr_0.3.2              assertthat_0.2.1        
## [79] ggplot2_3.1.1            xfun_0.6                
## [81] xtable_1.8-4             survival_2.44-1.1       
## [83] tibble_2.1.1             pheatmap_1.0.12         
## [85] AnnotationDbi_1.47.0     memoise_1.1.0           
## [87] brew_1.0-6               GSEABase_1.47.0

13 基金

该项目得到了美国国立卫生研究院（NIH）和国家科学基金会（NSF）的资助。

参考文献

H Backman, Tyler W, and Thomas Girke. 2016. “systemPipeR: NGS workflow and report generation environment.” BMC Bioinformatics 17 (1):388. https://doi.org/10.1186/s12859-016-1241-0.
Langmead, Ben, and Steven L Salzberg. 2012. “Fast Gapped-Read Alignment with Bowtie 2.” Nat. Methods 9 (4). Nature Publishing Group:357–59. https://doi.org/10.1038/nmeth.1923.
Love, Michael, Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for RNA-seq Data with DESeq2.” Genome Biol. 15 (12):550. https://doi.org/10.1186/s13059-014-0550-8.
Robinson, M D, D J McCarthy, and G K Smyth. 2010. “EdgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.” Bioinformatics 26 (1):139–40. https://doi.org/10.1093/bioinformatics/btp616.
Yu, Guangchuang, Li-Gen Wang, and Qing-Yu He. 2015. “ChIPseeker: An R/Bioconductor Package for ChIP Peak Annotation, Comparison and Visualization.” Bioinformatics 31 (14):2382–3. https://doi.org/10.1093/bioinformatics/btv145.
Zhang, Y, T Liu, C A Meyer, J Eeckhoute, D S Johnson, B E Bernstein, C Nussbaum, et al. 2008. “Model-Based Analysis of ChIP-Seq (MACS).” Genome Biol. 9 (9). https://doi.org/10.1186/gb-2008-9-9-r137.
Zhu, Lihua J, Claude Gazin, Nathan D Lawson, Hervé Pagès, Simon M Lin, David S Lapointe, and Michael R Green. 2010. “ChIPpeakAnno: A Bioconductor Package to Annotate ChIP-seq and ChIP-chip Data.” BMC Bioinformatics 11 (11~may):237. https://doi.org/10.1186/1471-2105-11-237.

systemPipeR 包中文翻译之 ChIP-seq 分析
原文英文链接：ChIP-Seq Workflow Template作者：Daniela Cassol (danie...
Bioconda软件安装神器：多版本并存、环境复制、环境导出
seq你想知道的全在这）、ChIP-seq分析（ChIP-seq基本分析流程）、单细胞测序分析 (重磅综述：三万...
ATAC-seq 分析（上）
之前写过ChIP-seq的基础分析：ChIP-seq数据分析（一）：从raw reads到peaks，和此文有很多...
TCseq：基因表达趋势分析
TCseq包可以应用于转录组，单细胞转录组，ATAC-seq，Chip-seq等的表观基因组时序型数据分析。该包主...
学员分享-Chip-seq 实战分析流程
Chip-seq 实战分析流程本次实践是根据Jimmy 和洲更两位大神的笔记，还有Y叔的Chipseeker包，...
ATAC-seq分析练习
这篇文章来练习一下ATAC-seq分析。ATAC-seq和CHIP-seq的分析非常相似，CHIP-seq检测的是...
ChIP-Seq分析之ChIPQC结果含义
刘小泽写于2020.6.20首先是ChIP-Seq分析的前言介绍部分：1：了解ChIP-seq的实验流程 2：继...
ChIP-seq数据比对注意事项
刘小泽写于2020.5.6首先是ChIP-seq分析的前言介绍部分：1：了解ChIP-seq的实验流程2：继续了解...
ChIP-seq数据质控与过滤
刘小泽写于2020.4.15首先是ChIP-seq分析的前言介绍部分：1：了解ChIP-seq的实验流程2：继续了...
chip-seq实验方法及数据分析
chip-seq是用来分析蛋白和DNA互作的方法。chip-seq结合染色质免疫共沉淀和高通量测序技术来找到DNA...

systemPipeR 包中文翻译之 ChIP-seq 分析

R包

内容目录（我自己加在前面的）

1 简单介绍

2 工作环境

3 Read 预处理

4 比对

5 coverage 数据的实用程序( Utilities )

6 使用 MACS2 进行 Peak calling

7 用基因组背景注 Peaks

8 交集的 Peaks 的 reads 计数

9 差异结合分析

10 GO富集分析

11 Motif 分析

12 版本信息

13 基金

参考文献

1 简单介绍

1.1 背景和目标

1.2 实验设计

2 工作环境

2.1 生成工作流程环境

2.2 运行工作流程

2.2.1 在计算节点运行 R

3 Read 预处理

3.1 targets 文件提供的实验定义

3.2 Read 质量筛选与修剪

3.3 FASTQ 质量报告

4 比对

4.1 使用 Bowtie2 进行 Reads 的比对

4.2 Reads 和比对信息统计

4.3 创建用于在 IGV 中查看 BAM 文件的链接

5 coverage 数据的实用程序( Utilities )

5.1 Rle 属性中储存的 coverage 信息

5.2 调整已比对上的 reads

5.3 Naive Peak calling

5.4 绘制定义区间的 coverage 信息

6 使用 MACS2 进行 Peak calling

6.1 在进行 Peak calling 之前将每个样本的重复文件 BAM 合并成一个 BAM 文件。

6.2 没有input/reference的Peak calling

6.3 有input/reference的Peak calling

6.4 鉴定共同的 Peaks

7 用基因组背景注 Peaks

7.1 使用 ChIPpeakAnno 包进行注释

7.2 使用 ChIPseeker 包来注释

8 交集的 Peaks 的 reads 计数

9 差异结合分析

10 GO富集分析

11 Motif 分析

11.1 从基因组中获得 Peaks 区间的DNA序列

11.2 使用 BCRANK 发现 Motif

12 版本信息

13 基金

参考文献

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

2.2.1 在计算节点运行 `R`

3 `Read` 预处理

3.1 `targets` 文件提供的实验定义

3.2 `Read` 质量筛选与修剪

3.3 `FASTQ` 质量报告

4.1 使用 `Bowtie2` 进行 `Reads` 的比对

4.2 `Reads` 和比对信息统计

4.3 创建用于在 `IGV` 中查看 `BAM` 文件的链接

5 `coverage` 数据的实用程序( Utilities )

5.1 `Rle` 属性中储存的 `coverage` 信息

5.2 调整已比对上的 `reads`

5.3 `Naive Peak calling`

5.4 绘制定义区间的 `coverage` 信息

6 使用 `MACS2` 进行 `Peak calling`

6.1 在进行 `Peak calling` 之前将每个样本的重复文件 `BAM` 合并成一个 `BAM` 文件。

6.2 没有`input/reference`的`Peak calling`

6.3 有`input/reference`的`Peak calling`

6.4 鉴定共同的 `Peaks`

7 用基因组背景注 `Peaks`

7.1 使用 `ChIPpeakAnno` 包进行注释

7.2 使用 `ChIPseeker` 包来注释

8 交集的 `Peaks` 的 `reads` 计数

11 `Motif` 分析

11.1 从基因组中获得 `Peaks` 区间的DNA序列

11.2 使用 `BCRANK` 发现 `Motif`