【表观调控实战】五、DEG分析与peaks注释分析

这里是佳奥，我们继续RNA-Seq分析的结果探究。

第三步，DEG-差异分析。

1 RNA-Seq结果分析：DEG

##导入矩阵
rm(list = ls())
options(stringsAsFactors = F)
a=read.table('../figure-01-check-gene-expression//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)

> table(group_list)
group_list
 PhoKO SppsKO     WT 
     3      4      3 

##第一步Firstly for DEseq2
##第一张图，归一化前后比较
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()
##第一张图绘制结束，文件名为DEseq_RAWvsNORM.png

DEseq_RAWvsNORM.png

归一化后表达量成一条线了。

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()

qc-heatmap-samples.png

##绘制火山图
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')

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')##第二个基因的上下调基因火山图

QQ截图20220824105516.png

##韦恩图
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)

QQ截图20220824110305.png

##中间不断调试，两个基因的相关度
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)

QQ截图20220824111304.png

2 ChIP-Seq结果分析：单个peaks注释

第四步，peaks-distribution。

先对peaks文件进行注释。

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)

QQ截图20220824131758.png

QQ截图20220824131806.png

自动化流程

 ---
title: "Cg_Wt"
author: "jmzeng1314@163.com"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: html_document
---
 
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(warning = F)
knitr::opts_chunk$set(message  = F) 
```

 * [我的博客](http://www.bio-info-trainee.com/)
 * [我们的论坛](http://www.biotrainee.com/forum.php)
 * [捐赠我](http://www.bio-info-trainee.com/donate)

##背景介绍
这里面描述一下背景

##读入peaks
 ```{r}
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' ))
```
 
##peaks性质

##ChIP Peaks over Chromosomes.
首先查看这些peaks在各个染色体的分布，全局浏览
```{r, fig.height=10} 
covplot(peak, weightCol="V5") 
```

##Heatmap of ChIP peaks binding to TSS regions
然后查看这些peaks在所有基因的启动子附近的分布情况，热图模式
```{r} 
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
tagMatrix <- getTagMatrix(peak, windows=promoter) 
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")
```

##Then Average Profile of ChIP peaks binding to TSS region
然后查看这些peaks在所有基因的启动子附近的分布情况，信号强度曲线图
```{r} 
plotAvgProf(tagMatrix, xlim=c(-3000, 3000), 
            xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency")
```

##peaks的注释
注释结果如下表，鼠标滑动可以查看全部详细信息：
```{r} 
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])
write.csv(peakAnno_df,paste0(sampleName,'_peakAnno_df.csv'))
DT::datatable(peakAnno_df,
                  extensions = 'FixedColumns',
                  options = list(
                    #dom = 't',
                    scrollX = TRUE,
                    fixedColumns = TRUE
                  ))
```

##可以对peaks的性质做一些可视化，如下：

```{r} 
#png('Pie-summarize the distribution of peaks over different type of features.png')
plotAnnoPie(peakAnno)
#png('Bar-summarize the distribution of peaks over different type of features.png')
plotAnnoBar(peakAnno)
#png('vennpie-summarize the distribution of peaks over different type of features.png')
#vennpie(peakAnno)
```

##还可以查看peaks的长度分布，只统计长度在1000bp以下的peaks

```{r}
peaksLength=abs(peakAnno_df$end-peakAnno_df$start)
peaksLength=peaksLength[peaksLength<1000]  
hist(peaksLength, breaks = 50, col = "lightblue", xlim=c(0,1000),xlab="peak length", main="Histogram of peak length") 
```

##peaks相关基因的注释
##这里可以把peaks先分类再注释，也可以直接拿所有peaks相关基因去富集分析，如果要分类，可以根据:

- Promoter
- 5’ UTR
- 3’ UTR
- Exon
- Intron
- Downstream
- Intergenic

##但是如果peaks本来就不多，那么分类后基因太少，注释可能并没有意义，这里只给出所有peaks相关基因的注释结果。

在R Studio中打开这个Cg_WT.Rmd文件，点击Knit便可以批量运行。

QQ截图20220824134424.png

分析结束后会出现一个.html的报表。（在当前目录）

QQ截图20220824134543.png

想要批量处理的话，修改如下：

##修改title
 ---
title: "Cg_Wt"
author: "jmzeng1314@163.com"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: html_document
---

##修改读取的文件
```{r}
bedPeaksFile        = 'oldBedFiles/Cg_WT.narrowPeak.bed'; 

##修改.Rmd文件名
Cg_Wt.Rmd

然后运行，得到新的.html报表。

3 批量peaks注释

##annotation-for-each-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)


##anno-for-all-peaks批量注释
##barplot
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" )

QQ截图20220824141446.png

下一篇我们继续专注peaks的分析。

我们下一篇再见！

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

【表观调控 实战】五、DEG分析与peaks注释分析

1 RNA-Seq结果分析：DEG

2 ChIP-Seq结果分析：单个peaks注释

自动化流程

3 批量peaks注释

推荐阅读更多精彩内容

【表观调控实战】五、DEG分析与peaks注释分析