由于之前沒有弄清楚這個(gè)步驟的具體目的奶甘，下來又重新看了一下，因?yàn)橐呀?jīng)拷貝了完整的數(shù)據(jù)钉赁，也有了前面call peaks的步驟完成，這一步就可以進(jìn)行了带膜，用R包ChIPQC對(duì)得到的peaks進(jìn)行質(zhì)控：
rm(list = ls())
options(stringsAsFactors = F)
library(ChIPQC)
# SampleID: 樣本ID
# Tissue, Factor, Condition: 不同的實(shí)驗(yàn)設(shè)計(jì)對(duì)照信息
# 三列信息必須包含在sampleSheet里咐扭，如果沒有某一列的信息設(shè)為NA袜爪。
# Replicate : 重復(fù)樣本的編號(hào)
# bamReads : 實(shí)驗(yàn)組BAM 文件的路徑（data/bams）
# ControlID : 對(duì)照組樣本ID
# bamControl :對(duì)照組樣本的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)
運(yùn)行到QCmetrics(exampleExp)報(bào)錯(cuò)
> 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]
解決方案應(yīng)該是去除染色體前面chr或者添加chr
另外今天在公眾號(hào)上看到一個(gè)推文ggVennDiagram 誕生記
就去下載了這個(gè)包醇疼，想重新畫一下韋恩圖
這個(gè)包需要一些向量或者列表直接畫圖倔毙，但是好像不可以對(duì)字符進(jìn)行作圖
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" ...
還是需要去好好讀一讀包的文檔再來做這個(gè)比較
