搞笑瞬间1：warning警告解决不了

搞笑瞬间2：就一个差异片段

我怀疑是样本的FRIP值太低了造成的

代码记录
library(DiffBind)
library(stats)

-------------构建包含bam文件和callpeak文件的数据框---------------

SampleID <- c(paste("Aging",1:3,sep = ""),paste("Young",1:3,sep = ""))
Tissue <- rep(c("Aging","Young"),c(3,3))
Factor <- rep("ATAC",times = 6)
Condition <- rep("same",times = 6)
Treatment <- rep("full-media",times = 6)
Replicate <- rep(1:3,length = 6)
PeakCaller <- rep("bed",times = 6)
bam_file_path <- "C:/文件夹数据重要一万分丢了就死人了/生信学习/otherdata/MLZ/callpeak"
bamReads <- c(paste(bam_file_path,"Aging1.last.bam",sep = "/"),
paste(bam_file_path,"Aging2.last.bam",sep = "/"),
paste(bam_file_path,"Aging3.last.bam",sep = "/"),
paste(bam_file_path,"Young1.last.bam",sep = "/"),
paste(bam_file_path,"Young2.last.bam",sep = "/"),
paste(bam_file_path,"Young3.last.bam",sep = "/"))
peak_file_path <- "C:/文件夹数据重要一万分丢了就死人了/生信学习/otherdata/MLZ/callpeak"
Peaks <- c(paste(peak_file_path,"Aging1_peaks.narrowPeak",sep = "/"),
paste(peak_file_path,"Aging2_peaks.narrowPeak",sep = "/"),
paste(peak_file_path,"Aging3_peaks.narrowPeak",sep = "/"),
paste(peak_file_path,"Young1_peaks.narrowPeak",sep = "/"),
paste(peak_file_path,"Young2_peaks.narrowPeak",sep = "/"),
paste(peak_file_path,"Young3_peaks.narrowPeak",sep = "/"))
samples <- data.frame(SampleID,Tissue,Factor,Condition,Treatment,Replicate,bamReads,Peaks,PeakCaller)

-----------------读取数据---------------------

ATAC <- dba(sampleSheet = samples,minOverlap = 2)

minoverlap表示该参数值为正整数n

表示将那些至少在n个样本中出现的peak纳入分析，其它peak舍弃

命令说明

1. 该命令将metadata中的所有样本信息一次性读入，并构建一个DBA对象

2. DBA对象可以理解为DiffBind用来存储信息的特定格式，本质上是一个S3类的列表，可以使用"$"符访问其中元素

3. 该DBA对象中存储大量信息，包括metadata中直接给出的，还有通过整合metadta得到的新的信息，其中重要的有"masks",在此例中，访问方式为h3k27ac$masks

4. 该命令在读入各样本peak之后，会根据minOverlap参数，去掉不符合条件的peak之后，再对所有剩余peak进行merge，有重叠的会被merge到一起形成新的peak，merge之后形成一个一致性peak（consensus peaksets）

5. 直接输入变量名可以查看该DBA对象的基本信息，要想查看该对象中所以信息，使用str(h3k27ac)

查看该DBA对象的基本信息

ATAC

6 Samples, 33984 sites in matrix (75392 total):

ID Tissue Factor Condition Treatment Replicate Intervals

1 Aging_1 Aging ATAC same full-media 1 18818

基本信息说明

1. 该信息中将peak文件和bam文件的位置隐去了，元数据中的其它信息都会展示在这里

2. 数据中添加了每个peak文件中的peak数目

3. 首行信息所表示的意思是：共读取了6个样本，将至少出现在三个样本中的所有peak进行merge之后共产生33984个peak，而在不筛选的情况下，所有peak在merge之后形成75392个peak

4. 虽然输入ATAC之后只显示了很少的信息，但实际上该数据中包含非常多的信息，使用str(ATAC)就可以查看所有信息

--------------------Occupancy analysis-------------------------------------

pdf("occupancy_plot.pdf")
plot(ATAC)
dev.off()

该热图所表示的样本之间的相关性是根据各样本的peak位置进行计算的

计算相关性时所用的peak，是根据dba函数中minOverlap参数读取并进行了merge之后的一致性peak

------------------Counting reads--------------------------------------------

ATAC <- dba.count(ATAC,minOverlap = 2)

这里的参数minOverlap的值，与dba中该参数的值可以不同，既可以比之前大，也可以比之前小。

查看基本信息

ATAC

6 Samples, 33869 sites in matrix:

ID Tissue Factor Replicate Caller Intervals FRiP

Aging_1 Aging ATAC same full-media 1 2615279 0.07

基本信息说明

1. 该信息与之前数据读取之后的结果十分相似，差别仅是多了一列内容

2. 最后一列表示在peak区域的reads在占所有bam文件中所有reads的比例，当然越高越高，表示背景越低

--------------根据reads count结果计算样本相关性并绘图--------------------------

pdf("affinity_plot.pdf")
plot(ATAC)
dev.off()

该热图所表示的样本之间的相关性是根据各样本在一致性peak上的信号强度（reads count数）计算的

------------------------------创建分组-----------------------------------

ATAC <- dba.contrast(ATAC,categories = DBA_TISSUE)
ATAC

Design: [~Tissue] | 1 Contrast:

Factor Group Samples Group2 Samples2

1 Tissue Young 3 Aging 3

多个分组信息

-------------------------------差异peak分析----------------------------------

ATAC <- dba.analyze(ATAC,method = DBA_ALL_METHODS)

默认情况是基于DEseq2, 可以设置参数method=DBA_EDGER选择edgeR，或者设置method=DBA_ALL_METHODS。每种方法都会评估差异结果的p-vaue和FDR。

用于差异分析的peak，来自于dba.count函数中进行了reads count计数的一致性peak

pdf(file="overlap_DESeq2_edgeR.pdf",width = 7,height = 7)
dba.plotVenn(ATAC,contrast=1,method=DBA_ALL_METHODS)
dev.off()
ATAC

6 Samples, 33096 sites in matrix:

Design: [~Tissue] | 1 Contrast:

Factor Group Samples Group2 Samples2 DB.edgeR DB.DESeq2

1 Tissue Young 3 Aging 3 200 1

找到的差异peak共有1个????但是DB。edgeR有200

默认情况下，peak具有显著差异的标准是 FDR<=0.05

comp1.edgeR <- dba.report(ATAC, method=DBA_EDGER, contrast = 1, th=1)
out <- as.data.frame(comp1.edgeR)
write.table(out, file="results_edgeR.txt", sep="\t", quote=F, col.names = NA)

edge.bed <- out[ which(out$FDR < 0.05),
c("seqnames", "start", "end", "strand", "Fold")]
write.table(edge.bed, file="results_edgeR_sig.bed", sep="\t", quote=F, row.names=F, col.names=F)

---------------根据差异peak绘图---------------------------

pdf("diff_plot.pdf")
plot(ATAC, contrast = 1)
dev.off()

------------结果报告和绘图------------------------------------

ATAC.DB <- dba.report(ATAC)
ATAC.DB.DF <- as.data.frame(ATAC.DB)
write.table(ATAC.DB.DF,file = "diff_peaks.txt",
quote = F,sep = "\t",
row.names = F,col.names = T)
pdf("volcano_plot.pdf")
dba.plotVolcano(ATAC)
dev.off()