Diffbind差异片段分析
搞笑瞬间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()
