由于之前没有弄清楚这个步骤的具体目的，下来又重新看了一下，因为已经拷贝了完整的数据，也有了前面call peaks的步骤完成，这一步就可以进行了，用R包ChIPQC对得到的peaks进行质控：

rm(list = ls())
options(stringsAsFactors = F)
library(ChIPQC)
# SampleID: 样本ID
# Tissue, Factor, Condition: 不同的实验设计对照信息
# 三列信息必须包含在sampleSheet里，如果没有某一列的信息设为NA。
# Replicate : 重复样本的编号
# bamReads : 实验组BAM 文件的路径（data/bams）
# ControlID : 对照组样本ID
# bamControl :对照组样本的bam文件路径
# Peaks :样本peaks文件的路径
# PeakCaller ：peak类型的字符串，可以是raw,bed,narrow,macs等。

(bams=list.files(path = '~/Downloads/fly/CHIP-SEQ/align/',
                 pattern = '*.raw.bam$', full.names = T))
bams=bams[grepl('WT',bams)]
peaks=list.files(path = '~/Downloads/fly/CHIP-SEQ/peaks/',pattern = '*_raw_peaks.narrowPeak', full.names = T)
peaks=peaks[grepl('WT',peaks)]

library(stringr)
SampleID=sub('.raw.bam','',basename(bams))
Replicate=str_split(basename(bams),'_',simplify = T)[,3]
Factor=str_split(basename(bams),'_',simplify = T)[,1]


samples=data.frame(SampleID=SampleID,
                   Tissue='WT', 
                   Factor=Factor, 
                   Replicate=1,            
                   bamReads=bams,                           
                   Peaks=peaks) 

exampleExp = ChIPQC(samples,
                    chromosomes='2L',
                    annotaiton="dm6")
QCmetrics(exampleExp)  #shows a summary of the main QC metrics
ChIPQCreport(exampleExp)

运行到QCmetrics(exampleExp)报错

> QCmetrics(exampleExp)  #shows a summary of the main QC metrics
Error in names(res) <- c("Reads", "Map%", "Filt%", "Dup%", "ReadL", "FragL",  : 
  'names' attribute [9] must be the same length as the vector [7]

解决方案应该是去除染色体前面chr或者添加chr

另外今天在公众号上看到一个推文ggVennDiagram 诞生记
就去下载了这个包，想重新画一下韦恩图
这个包需要一些向量或者列表直接画图，但是好像不可以对字符进行作图

rm(list = ls())
options(stringsAsFactors = F)
a=read.table('~/Downloads/10-12 tree/RNA/counts/all.counts.id.txt',header = T)
dim(a)
dat=a[,7:16]
rownames(dat)=a[,1]
dat[1:4,1:4]
library(stringr)
group_list=str_split(colnames(dat),'_',simplify = T)[,1]
table(group_list)

load(file = 'deg_output.Rdata')
library(ggpubr)
colnames(DEG_PhoKO)
DEG_PhoKO$log=log(DEG_PhoKO$baseMean+1)
DEG_PhoKO$change=ifelse(DEG_PhoKO$padj>0.05,'stable',
                        ifelse(DEG_PhoKO$log2FoldChange > 0,'up','down'))
table(DEG_PhoKO$change)
ggscatter(DEG_PhoKO,x="log" ,y="log2FoldChange",color = 'change')
 
DEG_SppsKO$log=log(DEG_SppsKO$baseMean+1)
DEG_SppsKO$change=ifelse(DEG_SppsKO$padj>0.05,'stable',
                        ifelse(DEG_SppsKO$log2FoldChange > 0,'up','down'))
table(DEG_SppsKO$change)
ggscatter(DEG_SppsKO,x="log" ,y="log2FoldChange",color = 'change')

library(UpSetR)

SppsKO_up=rownames(DEG_SppsKO[DEG_SppsKO$change=='up',])
SppsKO_down=rownames(DEG_SppsKO[DEG_SppsKO$change=='down',])
PhoKO_up=rownames(DEG_PhoKO[DEG_PhoKO$change=='up',])
PhoKO_down=rownames(DEG_PhoKO[DEG_PhoKO$change=='down',])

allG=unique(c(SppsKO_up,SppsKO_down,PhoKO_up,PhoKO_down))

df=data.frame(allG=allG,
              SppsKO_up=as.numeric(allG %in% SppsKO_up),
              SppsKO_down=as.numeric(allG %in% SppsKO_down),
              PhoKO_up=as.numeric(allG %in% PhoKO_up),
              PhoKO_down=as.numeric(allG %in% PhoKO_down))

upset(df)
library(gplots)

input <- list(SppsKO_up,SppsKO_down,PhoKO_up,PhoKO_down)
venn(input)
class(input[[1]])
str(input[[1]])
library(ggVennDiagram)
x <- c(input[[1]],input[[2]],input[[3]],input[[4]])
ggVennDiagram(x)

upset

gplots

> x <- c(input[[1]],input[[2]],input[[3]],input[[4]])
> library(ggVennDiagram)
> ggVennDiagram(x)
Error in ggVennDiagram(x) : Only support 2-4 dimension venn diagram.
> class(input[[1]])
[1] "character"
> x <- c(as.numeric(input[[1]]),as.numeric(input[[2]]),as.numeric(input[[3]]),as.numeric(input[[4]]))
Warning messages:
1: NAs introduced by coercion 
2: NAs introduced by coercion 
3: NAs introduced by coercion 
4: NAs introduced by coercion 
> library(ggVennDiagram)
> ggVennDiagram(x)
Error in ggVennDiagram(x) : Only support 2-4 dimension venn diagram.
> class(input[[1]])
[1] "character"
> str(input[[1]])
 chr [1:656] "miple1" "Ady43A" "ry" "subdued" "CG31676" "HLH3B" "CG18278" "pst" "stops" ...

还是需要去好好读一读包的文档再来做这个比较