2019-10-13-周末班-表观调控技术day1-2

本次线下班涉及内容NGS组学

RNA-seq
学徒第2周，RNA-seq数据分析实战训练 https://mubu.com/doc/38y7pmgzLg

ChIP-seq
学徒第4周，ChIP-seq数据分析实战训练 https://mubu.com/doc/11taEb9ZYg

ATAC-seq
学徒第5周-ATAC-seq数据处理实战 https://mubu.com/doc/2DG1mC2kdg
本次课程讲义：https://mubu.com/doc/3Bd4aieYug
参考在果蝇探索PRC复合物（逆向收费读文献2019-18）

准备工作

提前学习以下内容：

如果能下载就下载数据，网络原因我们拷贝数据练习

在Jimmy大神带领下学习文章图片复现

1，查看基因表达情况

rm(list = ls())
options(stringsAsFactors = F)
a=read.table('~/Downloads/10-12 tree/RNA/counts/all.counts.id.txt',header = T)
dim(a)

cg=a[a[,1]=='pho',7:16]
library(ggpubr)
library(stringr)
dat=data.frame(gene=as.numeric(cg),
               sample=names(cg),
               group=str_split(names(cg),'_',simplify = T)[,1]
               )
ggbarplot(dat,x='sample',y='gene',
          color = 'group',fill = 'group')

p <- ggbarplot(dat,x='sample',y='gene',
               color = 'group',fill = 'group')
pp <- p + theme(axis.title.x=element_text(face="italic"),
                axis.text.x = element_text(angle=50,vjust = 0.5))
pp

cg=a[a[,1]=='Spps',7:16]
library(ggpubr)
library(stringr)
dat=data.frame(gene=as.numeric(cg),
               sample=names(cg),
               group=str_split(names(cg),'_',simplify = T)[,1]
)
ggbarplot(dat,x='sample',y='gene',color = 'group',fill = 'group')

p <- ggbarplot(dat,x='sample',y='gene',
               color = 'group',fill = 'group')
pp <- p + theme(axis.title.x=element_text(face="italic"),
                axis.text.x = element_text(angle=50,vjust = 0.5))
pp

输出结果如下

image.png

由于对画图的包不熟练，下去还需要继续训练。
2019-10-13日修改,曾老师讲应该自行学习调整图例尽量让它好看一点，然后就下来学习了小洁老师的ggplot第八篇图片主题和格式修改

2查看相关性

rm(list = ls())
options(stringsAsFactors = F)
a=read.table('all.counts.id.txt',header = T)
dim(a)
dat=a[,7:16]

pheatmap::pheatmap(cor(dat))

library(stringr)
ac=data.frame(group=str_split(colnames(dat),'_',simplify = T)[,1])
rownames(ac)=colnames(dat)
M=cor(log(dat+1))
pheatmap::pheatmap(M,
                   annotation_col = ac) 
pheatmap::pheatmap(scale(M),
                   annotation_col = ac) 

cg=dat[,colnames(dat)[grepl('_1',colnames(dat))]]
library(ggpubr)
cg=log(cg+1)
pairs(cg)

直接将all.counts.id读入画热图

分组画热图

分组后进一步缩放美化

3，查看差异表达

做到这一步的时候我才知道我前面做的不对，之前都没有意识到，听课不是特别跟得上，之前读取数据我是把文件放在当前文件夹的，其实是可以读取其他目录的

rm(list = ls())
options(stringsAsFactors = F)
a=read.table('../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)

# Firstly for DEseq2

if(T){
  exprSet=dat
  suppressMessages(library(DESeq2)) 
  (colData <- data.frame(row.names=colnames(exprSet), 
                         group_list=group_list) )
  dds <- DESeqDataSetFromMatrix(countData = exprSet,
                                colData = colData,
                                design = ~ group_list)
  dds <- DESeq(dds)
  png("qc_dispersions.png", 1000, 1000, pointsize=20)
  plotDispEsts(dds, main="Dispersion plot")
  dev.off()
  
  rld <- rlogTransformation(dds)
  exprMatrix_rlog=assay(rld) 
  #write.csv(exprMatrix_rlog,'exprMatrix.rlog.csv' )
  
  normalizedCounts1 <- t( t(counts(dds)) / sizeFactors(dds) )
  # normalizedCounts2 <- counts(dds, normalized=T) # it's the same for the tpm value
  # we also can try cpm or rpkm from edgeR pacage
  exprMatrix_rpm=as.data.frame(normalizedCounts1) 
  head(exprMatrix_rpm)
  #write.csv(exprMatrix_rpm,'exprMatrix.rpm.csv' )
  
  png("DEseq_RAWvsNORM.png",height = 800,width = 800)
  par(cex = 0.7)
  n.sample=ncol(exprSet)
  if(n.sample>40) par(cex = 0.5)
  cols <- rainbow(n.sample*1.2)
  par(mfrow=c(2,2))
  boxplot(exprSet, col = cols,main="expression value",las=2)
  boxplot(exprMatrix_rlog, col = cols,main="expression value",las=2)
  hist(as.matrix(exprSet))
  hist(exprMatrix_rlog)
  dev.off()
  
  library(RColorBrewer)
  (mycols <- brewer.pal(8, "Dark2")[1:length(unique(group_list))])
  cor(as.matrix(exprSet))
  # Sample distance heatmap
  sampleDists <- as.matrix(dist(t(exprMatrix_rlog)))
  #install.packages("gplots",repos = "http://cran.us.r-project.org")
  library(gplots)
  png("qc-heatmap-samples.png", w=1000, h=1000, pointsize=20)
  heatmap.2(as.matrix(sampleDists), key=F, trace="none",
            col=colorpanel(100, "black", "white"),
            ColSideColors=mycols[group_list], RowSideColors=mycols[group_list],
            margin=c(10, 10), main="Sample Distance Matrix")
  dev.off()
  
  cor(exprMatrix_rlog) 
  
  table(group_list)
  res <- results(dds, 
                 contrast=c("group_list","SppsKO","WT"))
  resOrdered <- res[order(res$padj),]
  head(resOrdered)
  DEG_SppsKO=as.data.frame(resOrdered)
  DEG_SppsKO=na.omit(DEG_SppsKO)
 
  table(group_list)
  res <- results(dds, 
                 contrast=c("group_list","PhoKO","WT"))
  resOrdered <- res[order(res$padj),]
  head(resOrdered)
  DEG_PhoKO=as.data.frame(resOrdered)
  DEG_PhoKO=na.omit(DEG_PhoKO)
  save(DEG_PhoKO,DEG_SppsKO,file = 'deg_output.Rdata')
}
ac=data.frame(group=group_list)
rownames(ac)=colnames(dat)

这一步输出结果较多：

qc_dispersions

qc-heatmap-samples

DEseq_RAWvsNORM

跟前面一步的作用差不多，多方验证，数据是可靠的。
差异分析：

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)

load(file = 'deg_output.Rdata')
library(ggpubr)
DEG_PhoKO=DEG_PhoKO[rownames(DEG_SppsKO),]
po=data.frame(SppsKO=DEG_SppsKO$log2FoldChange,
              PhoKO=DEG_PhoKO$log2FoldChange)
ggscatter(po,'SppsKO','PhoKO')
sp <- ggscatter(po,'SppsKO','PhoKO',
                add = "reg.line",  # Add regressin line
                add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
                conf.int = TRUE # Add confidence interval
)
# Add correlation coefficient
sp + stat_cor(method = "pearson", label.x = 3, label.y = 30)

输出结果如下：

4,较正后的相关性

5,根据基因排序后的样本相关性做比较

然后是今天的作业：完成查找韦恩图中重叠peaks基因座信息并转换成.bed文件
根据曾老师的代码找到了peaks所在染色体信息坐标，需要转换成peaks所在promoter区域所对应的基因然后这一步还不会，先把前面的代码放上来，后面再联系一下

colnames(df)
# "SppsKO_up"   "SppsKO_down" "PhoKO_up"    "PhoKO_down" 
SppsKO_up_uniq=df[df[,2]==1 & df[,3]==0 & df[,4]==0 & df[,5]==0 ,1] 
SppsKO_down_uniq=df[df[,2]==0 & df[,3]==1 & df[,4]==0 & df[,5]==0 ,1]
PhoKO_up_uniq=df[df[,2]==0 & df[,3]==0 & df[,4]==1 & df[,5]==0 ,1]
PhoKO_down_uniq=df[df[,2]==0 & df[,3]==0 & df[,4]==0 & df[,5]==1 ,1]
SppsKO_up_PhoKO_up=df[df[,2]==1 & df[,3]==0 & df[,4]==1 & df[,5]==0 ,1]
SppsKO_down_PhoKO_down=df[df[,2]==0 & df[,3]==1 & df[,4]==0 & df[,5]==1 ,1]
  
library(org.Dm.eg.db)
t1=toTable(org.Dm.egCHRLOC)
t2=toTable(org.Dm.egCHRLOCEND)
t3=toTable(org.Dm.egSYMBOL)
 
#　作业 

require(TxDb.Dmelanogaster.UCSC.dm6.ensGene  )
txdb <- TxDb.Dmelanogaster.UCSC.dm6.ensGene
gs=as.data.frame(genes(txdb))
library(org.Dm.eg.db)
t100=toTable(org.Dm.egFLYBASE)
t101=toTable(org.Dm.egSYMBOL)
t200=merge(t100,t101,by='gene_id')
colnames(gs)[6]="flybase_id"
t300=merge(t200,gs,by="flybase_id")

SppsKO_up_uniq
SppsKO_up_bed=t300[match(SppsKO_up_uniq[SppsKO_up_uniq %in% t300$symbol ],
                      t300$symbol),4:6]
write.table(SppsKO_up_bed,sep = '\t',
            col.names = F,file = 'SppsKO_up.bed'
            ,quote = F,row.names = F)
SppsKO_down_uniq
SppsKO_down_bed=t300[match(SppsKO_down_uniq[SppsKO_down_uniq %in% t300$symbol ],
                         t300$symbol),4:6]
write.table(SppsKO_down_bed,sep = '\t',
            col.names = F,file = 'SppsKO_down.bed'
            ,quote = F,row.names = F)

4，call peaks

第一步先对得到的每个bed文件进行注释

bedPeaksFile= 'oldBedFiles/Cg_WT.narrowPeak.bed'; 
bedPeaksFile
## loading packages
require(ChIPseeker)
# BiocManager::install("TxDb.Dmelanogaster.UCSC.dm3.ensGene")
# BiocManager::install("org.Dm.eg.db")
require(TxDb.Dmelanogaster.UCSC.dm3.ensGene  )
txdb <- TxDb.Dmelanogaster.UCSC.dm3.ensGene
require(clusterProfiler) 
peak <- readPeakFile( bedPeaksFile )  
keepChr= !grepl('Het',seqlevels(peak)) 
seqlevels(peak, pruning.mode="coarse") <- seqlevels(peak)[keepChr]
cat(paste0('there are ',length(peak),' peaks for this data' ))
peakAnno <- annotatePeak(peak, tssRegion=c(-3000, 3000), 
                         TxDb=txdb, annoDb="org.Dm.eg.db") 
peakAnno_df <- as.data.frame(peakAnno)
sampleName=basename(strsplit(bedPeaksFile,'\\.')[[1]][1])
print(sampleName)
plotAnnoPie(peakAnno) 
plotAnnoBar(peakAnno)

输出结果如下

第二步对所有的peaks进行注释

require(ChIPseeker)
# BiocManager::install("TxDb.Dmelanogaster.UCSC.dm3.ensGene")
# BiocManager::install("org.Dm.eg.db")
require(TxDb.Dmelanogaster.UCSC.dm3.ensGene  )
txdb <- TxDb.Dmelanogaster.UCSC.dm3.ensGene
require(clusterProfiler) 

bedPeaksFile        = 'oldBedFiles/Cg_WT.narrowPeak.bed'; 
bedPeaksFile
anno_bed <- function(bedPeaksFile){
  peak <- readPeakFile( bedPeaksFile )  
  keepChr= !grepl('Het',seqlevels(peak)) 
  seqlevels(peak, pruning.mode="coarse") <- seqlevels(peak)[keepChr]
  cat(paste0('there are ',length(peak),' peaks for this data' ))
  peakAnno <- annotatePeak(peak, tssRegion=c(-3000, 3000), 
                           TxDb=txdb, annoDb="org.Dm.eg.db") 
  peakAnno_df <- as.data.frame(peakAnno)
  sampleName=basename(strsplit(bedPeaksFile,'\\.')[[1]][1])
  return(peakAnno_df)
}

f=list.files(path = 'oldBedFiles/',pattern = 'WT',full.names = T)
tmp = lapply(f, anno_bed)
head(tmp[[1]])
df=do.call(rbind,lapply(tmp, function(x){
  #table(x$annotation)
  num1=length(grep('Promoter',x$annotation))
  num2=length(grep("5' UTR",x$annotation))
  num3=length(grep('Exon',x$annotation))
  num4=length(grep('Intron',x$annotation))
  num5=length(grep("3' UTR",x$annotation))
  num6=length(grep('Intergenic',x$annotation))
  return(c(num1,num2,num3,num4,num5,num6 ))
}))
colnames(df)=c('Promoter',"5' UTR",'Exon','Intron',"3' UTR",'Intergenic')
rownames(df)=unlist(lapply(f, function(x){
  basename(strsplit(x,'\\.')[[1]][1])
}))
library(reshape2)
df2=melt(apply(df,1,function(x) x/sum(x)))
colnames(df2)=c('dis','sample','fraction')
library(ggpubr)
ggbarplot(df2, "sample", "fraction",
          fill = "dis", color = "dis", palette = "jco" )

image.png

此步骤的作业是重新阅读文献，找到筛选条件优化结果

5，重新找相关性

先用deeptools找到相关性表格的数据

# https://deeptools.readthedocs.io/en/develop/content/tools/multiBigwigSummary.html
# https://deeptools.readthedocs.io/en/develop/content/tools/plotCorrelation.html
multiBigwigSummary  bins -b  *WT*bw  -o wt_results.npz -p 8

plotCorrelation -in wt_results.npz  \
--corMethod spearman --skipZeros \
--plotTitle "Spearman Correlation of Read Counts" \
--whatToPlot heatmap --colorMap RdYlBu --plotNumbers \
--plotFileFormat pdf \
-o heatmap_SpearmanCorr_readCounts.pdf   \
--outFileCorMatrix SpearmanCorr_readCounts.tab

multiBigwigSummary  bins -b    *WT*bw  -o wt_merge_results.npz -p 8

plotCorrelation -in wt_merge_results.npz  \
--corMethod spearman --skipZeros \
--plotTitle "Spearman Correlation of Read Counts" \
--whatToPlot heatmap --colorMap RdYlBu --plotNumbers \
--plotFileFormat pdf \
-o merge_heatmap_SpearmanCorr_readCounts.pdf   \
--outFileCorMatrix merge_SpearmanCorr_readCounts.tab

这一步我还没能运行，因为对deeptools的multiBigwigSummary函数也不是太懂，就先用曾老师跑完的数据用R重现相关性分析

rm(list = ls())
options(stringsAsFactors = F)
a=read.table('SpearmanCorr_readCounts.tab')
pheatmap::pheatmap(a)
library(stringr)
ac=data.frame(group=str_split(rownames(a),'_',simplify = T)[,1])
rownames(ac)=colnames(a)
M=a
pheatmap::pheatmap(M,
                   annotation_col = ac) 
 
a=read.table('merge_SpearmanCorr_readCounts.tab')
pheatmap::pheatmap(a)

得到结果如下

只有Ez_WT相关性是最好的

merge后的相关性

6, ChIPQC,由于前面的步骤没有走完，这一步也没有完成

library(ChIPQC)
samples = read.csv(file.path(system.file("extdata", package="ChIPQC"),
                             "example_QCexperiment.csv"))
samples
exampleExp = ChIPQC(samples,annotaiton="hg19")
QCmetrics(exampleExp)  #shows a summary of the main QC metrics
ChIPQCreport(exampleExp)

不太清楚是不是这个包不完整
报错如下

> library(ChIPQC)
> samples = read.csv(file.path(system.file("extdata", package="ChIPQC"),
+                              "example_QCexperiment.csv"))
> samples
  SampleID  Tissue Factor Replicate            bamReads                           Peaks
1   CTCF_1    A549   CTCF         1 reads/SRR568129.bam peaks/SRR568129_chr22_peaks.bed
2   CTCF_2    A549   CTCF         2 reads/SRR568130.bam peaks/SRR568130_chr22_peaks.bed
3   cMYC_1    A549   cMYC         1 reads/SRR568131.bam peaks/SRR568131_chr22_peaks.bed
4   cMYC_2    A549   cMYC         2 reads/SRR568132.bam peaks/SRR568132_chr22_peaks.bed
5   E2F1_1 HeLa-S3   E2F1         1 reads/SRR502355.bam peaks/SRR502355_chr22_peaks.bed
6   E2F1_2 HeLa-S3   E2F1         2 reads/SRR502356.bam peaks/SRR502356_chr22_peaks.bed
> exampleExp = ChIPQC(samples,annotaiton="hg19")
CTCF_1 A549 CTCF   1 bed
Error in if (file.info(peaks)$size > 0) { : 
  missing value where TRUE/FALSE needed

需要进一步学习这个包
尝试注释掉这个例子继续向下运行没有运行成功

library(ChIPQC)
# samples = read.csv(file.path(system.file("extdata", package="ChIPQC"),
#                              "example_QCexperiment.csv"))
# samples
# exampleExp = ChIPQC(samples,annotaiton="hg19")
# QCmetrics(exampleExp)  #shows a summary of the main QC metrics
# ChIPQCreport(exampleExp)

# 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/align/peaks-single/',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)

由于没有运行call peaks的上游分析path = '~/Downloads/fly/CHIP-SEQ/align/peaks-single/'这个路径中缺少上一步生成的peaks-single/文件
继续进行下一步，因为曾老师把生成的数据传输给我们了

7.重新注释

bedPeaksFile= 'Ez_SppsKO_1_raw_peaks.narrowPeak'; 
bedPeaksFile
## loading packages
require(ChIPseeker)
# BiocManager::install("TxDb.Dmelanogaster.UCSC.dm6.ensGene")
# BiocManager::install("org.Dm.eg.db")
require(TxDb.Dmelanogaster.UCSC.dm6.ensGene  )
txdb <- TxDb.Dmelanogaster.UCSC.dm6.ensGene
require(clusterProfiler) 
peak <- readPeakFile( bedPeaksFile ) 
seqlevels(peak)
keepChr= seqlevels(peak) %in% c('2L','2R','3L','3R','4','X','Y','M')
seqlevels(peak, pruning.mode="coarse") <- seqlevels(peak)[keepChr]
seqlevels(peak, pruning.mode="coarse") <- paste0('chr',seqlevels(peak))
seqlevels(peak)
cat(paste0('there are ',length(peak),' peaks for this data' ))
peakAnno <- annotatePeak(peak, tssRegion=c(-3000, 3000), 
                         TxDb=txdb, annoDb="org.Dm.eg.db") 
peakAnno_df <- as.data.frame(peakAnno)
sampleName=basename(strsplit(bedPeaksFile,'\\.')[[1]][1])
print(sampleName)
plotAnnoPie(peakAnno) 
plotAnnoBar(peakAnno)

结果如下

require(ChIPseeker) 
require(TxDb.Dmelanogaster.UCSC.dm6.ensGene  )
txdb <- TxDb.Dmelanogaster.UCSC.dm6.ensGene
require(clusterProfiler)  
anno_bed <- function(bedPeaksFile){
  peak <- readPeakFile( bedPeaksFile )  
  seqlevels(peak)
  keepChr= seqlevels(peak) %in% c('2L','2R','3L','3R','4','X','Y','M')
  seqlevels(peak, pruning.mode="coarse") <- seqlevels(peak)[keepChr]
  seqlevels(peak, pruning.mode="coarse") <- paste0('chr',seqlevels(peak))
  seqlevels(peak)
  cat(paste0('there are ',length(peak),' peaks for this data' ))
  peakAnno <- annotatePeak(peak, tssRegion=c(-3000, 3000), 
                           TxDb=txdb, annoDb="org.Dm.eg.db") 
  peakAnno_df <- as.data.frame(peakAnno)
  sampleName=basename(strsplit(bedPeaksFile,'\\.')[[1]][1])
  return(peakAnno_df)
}

f=list.files(path = './',pattern = 'WT',full.names = T)
f
tmp = lapply(f, anno_bed)
head(tmp[[1]])
df=do.call(rbind,lapply(tmp, function(x){
  #table(x$annotation)
  num1=length(grep('Promoter',x$annotation))
  num2=length(grep("5' UTR",x$annotation))
  num3=length(grep('Exon',x$annotation))
  num4=length(grep('Intron',x$annotation))
  num5=length(grep("3' UTR",x$annotation))
  num6=length(grep('Intergenic',x$annotation))
  return(c(num1,num2,num3,num4,num5,num6 ))
}))
colnames(df)=c('Promoter',"5' UTR",'Exon','Intron',"3' UTR",'Intergenic')
rownames(df)=unlist(lapply(f, function(x){ 
  (strsplit(basename(x),'\\.')[[1]][1])
}))
library(reshape2)
df2=melt(apply(df,1,function(x) x/sum(x)))
colnames(df2)=c('dis','sample','fraction')
library(ggpubr)
ggbarplot(df2, "sample", "fraction",
          fill = "dis", color = "dis", palette = "jco"  )

重注释后的分布

8, 重新画韦恩图

rm(list=ls())
require(ChIPseeker) 
require(TxDb.Dmelanogaster.UCSC.dm3.ensGene  )
txdb <- TxDb.Dmelanogaster.UCSC.dm3.ensGene
require(clusterProfiler) 
(bedFiles=list.files(pattern = '*bed'))

library(ChIPpeakAnno)
fs=lapply(bedFiles,function(x){
  peak <- readPeakFile( x )  
  keepChr= !grepl('Het',seqlevels(peak)) 
  seqlevels(peak, pruning.mode="coarse") <- seqlevels(peak)[keepChr]
  peak
})
library(stringr)
ol <- findOverlapsOfPeaks(fs[[1]],fs[[2]],fs[[3]] )
png('3_factors_overlapVenn.png')
makeVennDiagram(ol,
                NameOfPeaks=str_split(bedFiles,'_',simplify = T)[,1],
                TxDb=txdb)
dev.off()

image.png

作业：由于没有好好读文献，之前用的参考基因组不是文章中所用的需要进一步优化
进行基因组的坐标转换，这一步真的有点难，一个是转换的方式比较多，有点混淆，另一个是对果蝇的基因组不熟悉。
明天再继续吧……

最后编辑于：2019.10.13 21:17:28

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 204,189评论 6赞 478
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 85,577评论 2赞 381
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 150,857评论 0赞 337
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 54,703评论 1赞 276
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 63,705评论 5赞 366
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,620评论 1赞 281
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 37,995评论 3赞 396
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,656评论 0赞 258
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 40,898评论 1赞 298
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,639评论 2赞 321
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,720评论 1赞 330
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,395评论 4赞 319
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 38,982评论 3赞 307
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,953评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,195评论 1赞 260
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 44,907评论 2赞 349
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,472评论 2赞 342

2019-10-13-周末班-表观调控技术day1-2

准备工作

在Jimmy大神带领下学习文章图片复现

1，查看基因表达情况

2查看相关性

3，查看差异表达

4，call peaks

5，重新找相关性

6, ChIPQC,由于前面的步骤没有走完，这一步也没有完成

7.重新注释

8, 重新画韦恩图

推荐阅读更多精彩内容