這里是佳奧卵迂，我們繼續(xù)RNA-Seq分析的結(jié)果探究裕便。
第三步，DEG-差異分析见咒。
1 RNA-Seq結(jié)果分析：DEG
##導(dǎo)入矩陣
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()
##第一張圖繪制結(jié)束挂疆，文件名為DEseq_RAWvsNORM.png
DEseq_RAWvsNORM.png
歸一化后表達(dá)量成一條線了改览。
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')##第一個(gè)基因的上下調(diào)基因火山圖
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')##第二個(gè)基因的上下調(diào)基因火山圖
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
##中間不斷調(diào)試缤言，兩個(gè)基因的相關(guān)度
load(file = 'deg_output.Rdata')
library(ggpubr)
DEG_PhoKO=DEG_PhoKO[rownames(DEG_SppsKO),]##兩個(gè)數(shù)據(jù)的基因縱列順序要一一對(duì)應(yīng)
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結(jié)果分析：?jiǎn)蝹€(gè)peaks注釋
第四步宝当，peaks-distribution。
先對(duì)peaks文件進(jìn)行注釋胆萧。
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
自動(dòng)化流程
---
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)
* [捐贈(zèng)我](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性質(zhì)
##ChIP Peaks over Chromosomes.
首先查看這些peaks在各個(gè)染色體的分布庆揩，全局瀏覽
```{r, fig.height=10}
covplot(peak, weightCol="V5")
```
##Heatmap of ChIP peaks binding to TSS regions
然后查看這些peaks在所有基因的啟動(dòng)子附近的分布情況，熱圖模式
```{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在所有基因的啟動(dòng)子附近的分布情況跌穗，信號(hào)強(qiáng)度曲線圖
```{r}
plotAvgProf(tagMatrix, xlim=c(-3000, 3000),
xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency")
```
##peaks的注釋
注釋結(jié)果如下表订晌，鼠標(biāo)滑動(dòng)可以查看全部詳細(xì)信息：
```{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
))
```
##可以對(duì)peaks的性質(zhì)做一些可視化，如下：
```{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的長(zhǎng)度分布蚌吸，只統(tǒng)計(jì)長(zhǎng)度在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相關(guān)基因的注釋
##這里可以把peaks先分類再注釋锈拨，也可以直接拿所有peaks相關(guān)基因去富集分析，如果要分類羹唠，可以根據(jù):
- Promoter
- 5’ UTR
- 3’ UTR
- Exon
- Intron
- Downstream
- Intergenic
##但是如果peaks本來就不多奕枢，那么分類后基因太少，注釋可能并沒有意義佩微，這里只給出所有peaks相關(guān)基因的注釋結(jié)果缝彬。
在R Studio中打開這個(gè)Cg_WT.Rmd文件，點(diǎn)擊Knit便可以批量運(yùn)行哺眯。
QQ截圖20220824134424.png
分析結(jié)束后會(huì)出現(xiàn)一個(gè).html的報(bào)表谷浅。（在當(dāng)前目錄）
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
然后運(yùn)行，得到新的.html報(bào)表壳贪。
3 批量peaks注釋
##annotation-for-each-bed單個(gè)注釋
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
下一篇我們繼續(xù)專注peaks的分析陵珍。
我們下一篇再見！
