1. 获得peak集

image.png

这里的逻辑是：把四个样品合并，call peaks，然后获得peaks文件与前面idr 处理后的peaks进行overlap，都overlap的peaks，作为最终的peaks。

gsize=199000000

## callpeak and idr analysis of sample A
s1_r1=SRR6322534
s1_r2=SRR6322535
s1=SRR6322534_SRR6322535
input1=SRR6322538

# call peaks for replicte 1
macs2 callpeak -t ./${s1_r1}.ff.bam -c ./$input1.ff.bam -f BAM -n ${s1_r1} -g $gsize --keep-dup all -p 0.01

# call peaks for replicte 2
macs2 callpeak -t ./${s1_r2}.ff.bam -c ./$input1.ff.bam -f BAM -n ${s1_r2} -g $gsize --keep-dup all -p 0.01

# call peaks for combined dataset
macs2 callpeak -t ./${s1_r1}.ff.bam ./${s1_r2}.ff.bam -c ./$input1.ff.bam -f BAM -n ${s1} -g $gsize --keep-dup all -p 0.01

# run idr
idr --samples ${s1_r1}_peaks.narrowPeak ${s1_r2}_peaks.narrowPeak --peak-list ${s1}_peaks.narrowPeak --input-file-type narrowPeak --output-file ./${s1}.idr --rank p.value --soft-idr-threshold 0.05 --plot

# get idr produced narrowPeak file
cut -f 1-10 ${s1}.idr | sort -k1,1 -k2,2n -k3,3n >${s1}.idr.narrowPeak


## callpeak and idr analysis of sample B
s2_r1=SRR8423051
s2_r2=SRR8423052
s2=SRR8423051_SRR8423052
input2=SRR8423055

# call peaks for replicte 1
macs2 callpeak -t ./${s2_r1}.ff.bam -c ./$input2.ff.bam -f BAM -n ${s2_r1} -g $gsize --keep-dup all -p 0.01

# call peaks for replicte 2
macs2 callpeak -t ./${s2_r2}.ff.bam -c ./$input2.ff.bam -f BAM -n ${s2_r2} -g $gsize --keep-dup all -p 0.01

# call peaks for combined dataset
macs2 callpeak -t ./${s2_r1}.ff.bam ./${s2_r2}.ff.bam -c ./$input2.ff.bam -f BAM -n ${s2} -g $gsize --keep-dup all -p 0.01

# run idr
idr --samples ${s2_r1}_peaks.narrowPeak ${s2_r2}_peaks.narrowPeak --peak-list ${s2}_peaks.narrowPeak --input-file-type narrowPeak --output-file ./${s2}.idr --rank p.value --soft-idr-threshold 0.05 --plot

# get idr produced narrowPeak file
cut -f 1-10 ${s2}.idr | sort -k1,1 -k2,2n -k3,3n >${s2}.idr.narrowPeak

## Peaks in combined sample bams
macs2 callpeak -t ${s1_r1}.ff.bam ${s1_r2}.ff.bam ${s2_r1}.ff.bam ${$s2_r2}.ff.bam -c $input1.ff.bam $input2.ff.bam -f BAM -n ${s1}_${s2} -g $gsize --keep-dup all -p 0.01

获得最终的peaks文件

Cat ${s1_r1}_${s1_r2}.idr.narrowPeak ${s2_r1}_${s2_r2}.idr.narrowPeak | sort -k1,1 -k2,2n -k3,3n | bedtools intersect -a ${s1}_${s2}_peaks.narrowPeak -b - -F 0.5 -u >${s1}_${s2}_peaks.f.narrowPeak

## filter peaks against blacklist
bedtools intersect -a ${s1}_${s2}_peaks.f.narrowPeak -b ../ref/blacklist.bed -f 0.25 -v >temp && mv temp ${s1}_${s2}_peaks.f.narrowPeak

image.png

2. 统计peaks区域的counts，主要利用deeptools的 multiBamSummary

s1=SRR6322534
s2=SRR6322535
s3=SRR8423051
s4=SRR8423052
peak=./peak.narrowPeak

## for comparison between samples without replicate
cut -f 1-3 $peak >peak.1.bed  #把peaks转换成 bed 文件
multiBamSummary BED-file --bamfiles ../callpeak/downsample/$s2.ff.bam ../callpeak/downsample/$s3.ff.bam --BED peak.1.bed --smartLabels -p 4 --outRawCounts peak_counts.1.txt --extendReads 134


## for comparison between samples with replicates
cut -f 1-3 $peak >peak.2.bed
multiBamSummary BED-file --bamfiles ../callpeak/downsample/$s1.ff.bam ../callpeak/downsample/$s2.ff.bam ../callpeak/downsample/$s3.ff.bam ../callpeak/downsample/$s4.ff.bam --BED peak.2.bed --smartLabels -p 4 --outRawCounts peak_counts.2.txt --extendReads 134

The multiBamSummary command is part of the deepTools package and is used to generate summary statistics for multiple BAM files. Here's an explanation of the options in the example command you provided:

BED-file: The name of the BED file containing the genomic regions of interest.
--bamfiles: A list of BAM files to analyze.
--smartLabels: Automatically generate labels for the BAM files based on their file names.
-p: The number of threads to use for parallel processing.
--outRawCounts: The name of the output file containing the raw counts for each genomic region.
--extendReads: The number of base pairs to extend the reads in each direction.

3. R里面进行差异分析

1. read in sample information

library(DESeq2)
library(tidyverse)
library(pheatmap)

cat meta.2.txt
SRR6322534  rabbit, anti-IR_beta, sc-711    HepG2_IRb   19272489
SRR6322535  rabbit, anti-IR_beta, sc-711    HepG2_IRb   37475005
SRR8423051  rabbit, anti-IR_beta, sc-711    HepG2_IRb_starve    28559621
SRR8423052  rabbit, anti-IR_beta, sc-711    HepG2_IRb_starve    45616747

meta <- read_tsv("meta.2.txt", col_names = F)
meta <- meta[,c(1,3)]
colnames(meta) <- c("id", "group")
meta <- meta %>% column_to_rownames(var = "id")
meta$group <- factor(meta$group, levels = c("HepG2_IRb", "HepG2_IRb_starve"))
head(meta)
                    group
SRR6322534        HepG2_IRb
SRR6322535        HepG2_IRb
SRR8423051        HepG2_IRb_starve
SRR8423052        HepG2_IRb_starve

2. read in counts

counts <- read_tsv("peak_counts.2.txt", col_names = F, comment = "#")
head(counts)

# A tibble: 6 × 7
  X1          X2       X3    X4    X5    X6    X7
  <chr>    <dbl>    <dbl> <dbl> <dbl> <dbl> <dbl>
1 chr1   8878506  8879058    55    38    49    30
2 chr1   8880265  8880537    25    16    14    26
3 chr19  9768598  9768832    10    15    10    16
4 chr1   7942428  7942608    19     4    14    12
5 chr1   7961280  7961786    44    29    30    29
6 chr17 78187231 78187636    44    23    28    20



colnames(counts) <- c("chr", "start", "end", rownames(meta))
counts.dds <- counts[,-(1:3)] %>% as.data.frame()
rownames(counts.dds) <- paste(counts$chr, counts$start, counts$end, sep="-")
head(counts.dds)

                          SRR6322534   SRR6322535   SRR8423051   SRR8423052
chr1-8878506-8879058            55         38         49         30
chr1-8880265-8880537            25         16         14         26
chr19-9768598-9768832           10         15         10         16
chr1-7942428-7942608            19          4         14         12
chr1-7961280-7961786            44         29         30         29
chr17-78187231-78187636         44         23         28         20

3. calculate size factor

library(csaw)
dpath <- paste("../../chip/callpeak/downsample/", rownames(meta), ".ff.bam", sep="")
bins <- windowCounts(dpath, bin=T, width=10000, BPPARAM=MulticoreParam(nrow(meta)))
nf <-normFactors(bins, se.out = F)

4. create dds

dds <- DESeqDataSetFromMatrix(countData=counts.dds, colData=meta, design=~group)
#dds <- estimateSizeFactors(dds) #这个是DESE2自带的size factor计算工具，由于我们提供给DeSEq2的是peaks里面的reads counts，因此不适合利用其自带的size factore 函数计算，当然我们如果提供了全基因组以bin'单位的raw reads counts是可以的，但是后续分析还得转换成peaks对应的counts进行分析，这样比较麻烦一点。

sizeFactors(dds) <- nf #前面的size factor赋值

5. 计算差异 peaks

dds <- DESeq(dds)
res <- results(dds)

存取数据

res.out <- res %>% as.data.frame() %>% rownames_to_column(var = "name") %>%
  mutate(
    change = case_when(
      padj >= 0.05 | is.na(padj) ~ "stable",
      padj < 0.05 & log2FoldChange < 0 ~ "down",
      padj < 0.05 & log2FoldChange > 0 ~ "up")
  )
res.out <- data.frame(counts[,1:3], res.out)
write.table(res.out, file = "diff.2.txt", sep = "\t", quote = F, row.names = F, col.names = T)

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

---

如果数据计算重复性欠佳，可以考虑用下面的方法对数据进一步标准化，进行相关性计算，

## signal transformation
rld <- rlog(dds, blind=F)
rldMat <- assay(rld)

The rlog function in the DESeq2 package performs a variance-stabilizing transformation on the count data in a DESeqDataSet object. The resulting transformed data can be used for downstream analysis such as clustering, visualization, and differential expression analysis. The blind argument in the rlog function controls whether the transformation is performed in a "blind" manner, meaning that the sample information is not used to estimate the dispersion parameters. By default, blind is set to TRUE, which means that the dispersion parameters are estimated from the data without taking into account the sample information. Setting blind to FALSE allows the dispersion parameters to be estimated using the sample information, which can improve the accuracy of the transformation.

The assay function in the SummarizedExperiment package extracts the assay data from a SummarizedExperiment object. In the case of a DESeqDataSet object that has been transformed using the rlog function, the assay data is the transformed count data. The assay function in the SummarizedExperiment package extracts the assay data from a SummarizedExperiment object. In the case of a DESeqDataSet object that has been transformed using the rlog function, the assay data is the transformed count data.

# distance between samples
png(file="dist.png", width = 400, height = 400)
pheatmap(as.matrix(dist(t(rldMat))), cluster_rows = T, cluster_cols = T)
dev.off()

image.png

pca analysis

plotPCA(rld, intgroup = "group")

image.png

correlation between samples using transformed data

pheatmap(cor(rldMat), cluster_rows = T, cluster_cols = T, display_numbers = T)
# scatter plot

image.png

ggplot(as.data.frame(rldMat), aes(x = SRR8423051, y = SRR8423052)) +
  geom_point(size = 0.5) + xlim(0,8) + ylim(0,8) + theme_cowplot()

image.png

correlation between samples using normalized counts (DeSeq2 按照size factor normalized，很多文章用的都是这个数据作图的，感觉效果不如transformed data好)

counts.norm <- counts(dds, normalized = T)
counts.norm.cor <- cor(log2(counts.norm), method = "pearson")
png(file="cor.png", width = 400, height = 400)
pheatmap(counts.norm.cor, cluster_rows = T, cluster_cols = T, display_numbers = T)
dev.off()

ggplot(as.data.frame(log2(counts.norm)), aes(x = SRR8423051, y = SRR8423052)) +
  geom_point(size = 0.5) + xlim(0,8) + ylim(0,8) + theme_cowplot()
save.image("diff.2.RData")

image.png

参考：基因课------表观基因组学之 ChIP-seq 数据分析