RNA-seq下游分析-1

database <- read.table(file = "D:/huage1/rnaseq1014/xiaohe/output.matrix", sep = "\t", header = T, row.names = 1)
library(magrittr)
library(dplyr)
database$ensembl_id<-rownames(database)
database$ensembl_id <- as.character(database$ensembl_id)
filtered_database <- database %>%  
  filter(!grepl("_PAR_Y$", database$ensembl_id))
names(filtered_database)[names(filtered_database) == 'ensembl_id'] <- 'v1'
library(tidyr)
database1 <- separate(filtered_database,v1,into = "ensembl_id",sep = "[.]")
names(database1)[5] <- "ensembl_id"
rownames(database1)=database1$ensembl_id
database2=database1[,c(1,2,3,4)]
colnames(database2) <- c("BHLHE40-rep1","BHLHE40-rep2","Control-rep1","Control-rep2")
library(DESeq2)
library(stringr)
library(dbplyr)
#BiocManager::install("DESeq2")
#BiocManager::install("dbplyr")
database <- round(as.matrix(database2))
condition <- factor(c(rep("treat",2),rep("control",2)))
coldata <- data.frame(row.names =colnames(database),condition)
library(tidyr)
dds <- DESeqDataSetFromMatrix(countData = database,colData = coldata,design = ~condition)
##countData用于说明数据来源，colData用于说明不同组数据的实验操作类型，design用于声明自变量，即谁和谁进行对比
nrow(dds)

## [1] 61498

dds2 <- dds[rowSums(counts(dds))>1,]
nrow(dds2)

## [1] 23051

#上面这步操作的目的是为了筛选数据。这步操作可以删除那些在所有样本中的表达计数总和不大于1的基因，因为这些基因的表达水平可能过于微弱，对于后续的分析可能不具有重要性或者可靠性。通过这种方式，可以减少噪声和潜在的误差，提高数据分析的准确性。
#聚类图
vsd <- vst(dds2, blind = FALSE)
sampleDists <- dist(t(assay(vsd)))
hc <- hclust(sampleDists, method = "ward.D2")
plot(hc, hang = -1)

#BiocManager::install("factoextra")
library(factoextra)
#install.packages("ggplot2")
res <- hcut(sampleDists, k = 2, stand = TRUE)
# Visualize 聚类
fviz_dend(res, 
          # 加边框
          rect = TRUE,
          # 边框颜色
          rect_border="cluster",
          # 边框线条类型
          rect_lty=2,
          # 边框线条粗细
          lwd=1.2,
          # 边框填充
          rect_fill = T,
          # 字体大小
          cex = 1,
          # 字体颜色
          color_labels_by_k=T,
          # 平行放置
          horiz=T)

vsd <- vst(dds2, blind = FALSE)
plotPCA(vsd, intgroup = "condition")

dds3 <- DESeq(dds2)#用于对dds数据进行运算及分析
#消除样本测序深度影响
normalized_counts <- as.data.frame(counts(dds3, normalized=TRUE))
#看一下基因是否发生变化
plotCounts(dds3, gene = "ENSG00000000003", intgroup=c("condition"))

#美化一下图
plotdata <- plotCounts(dds3, gene = "ENSG00000000003", intgroup=c("condition"),returnData = T)
library(ggplot2)
ggplot(plotdata,aes(x=condition,y=count,col=condition))+
  geom_point()+
  theme_bw()

library(DESeq2)
#MA图
dds4 <- results(dds3,contrast = c("condition","treat","control"),alpha = 0.05)
plotMA(dds4, ylim=c(-2,2))

#我们发现在左侧，有很多counts很小的基因，发生了很大的变化，但是没有明显意义。他们的counts很小，但波动性很大，对logFC产生了很大的影响。
#矫正后的MA图 在这句代码中，dd2 <- lfcShrink(dds, contrast=contrast, res=dd1)，lfcShrink是一个函数，它对数据集dds进行某种形式的"收缩"处理。这种处理可能涉及到统计假设检验中的标准化或者归一化等步骤。
#数据收缩：lfcShrink：两种数据特别需要：低表达量占比高的 & 数据特别分散的：
#normal 是DESeq2包原始的收缩估计量（shrikage estimator），自适应正态先验分布（adaptive normal prior）
#apeglm是apeglm包中的收缩估计量，自适应t先验分布（adaptive t prior）
#ashr是ashr包中的收缩估计量。
#BiocManager::install("apeglm")
library(apeglm)
#第一种校正方法
#resAsh <- lfcShrink(dds, coef="group_list_1_vs_0", type="ashr")
#plotMA(resAsh, ylim=c(-3,3))

##第二种
#resLFC <- lfcShrink(dds3, coef=2)
#plotMA(resLFC, ylim=c(-3,3))
#适合本文校正方法
resNorm <- lfcShrink(dds3, coef=2, type="normal")
plotMA(resNorm, ylim=c(-3,3))

resultsNames(dds3)

## [1] "Intercept"                  "condition_treat_vs_control"

summary(resNorm, alpha = 0.05)

## 
## out of 23051 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 49, 0.21%
## LFC < 0 (down)     : 29, 0.13%
## outliers [1]       : 0, 0%
## low counts [2]     : 8044, 35%
## (mean count < 11)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

#把差异分析的结果转化成data.frame的格式
library(dplyr)
library(tibble)
res <- resNorm %>% 
  data.frame() %>% 
  rownames_to_column("ensembl_id")
#基因id转换
library(AnnotationDbi)
library(org.Hs.eg.db)
res$symbol <- mapIds(org.Hs.eg.db,
                     keys=res$ensembl_id,
                     column="SYMBOL",
                     keytype="ENSEMBL",
                     multiVals="first")
res$entrez <- mapIds(org.Hs.eg.db,
                     keys=res$ensembl_id,
                     column="ENTREZID",
                     keytype="ENSEMBL",
                     multiVals="first")
#制作genelist
library(dplyr)
gene_df <- res %>% 
  dplyr::select(ensembl_id,log2FoldChange,symbol,entrez) %>% 
  ## 去掉NA
  filter(entrez!="NA") %>% 
  ## 去掉重复
  distinct(entrez,.keep_all = T)

geneList <- gene_df$log2FoldChange
names(geneList) = gene_df$entrez
geneList = sort(geneList, decreasing = TRUE)
head(geneList)

##     4210    10628    23581    56667      720     8322 
## 1.282975 1.257111 1.175282 1.062406 1.059902 1.022524

#GSEA分析
library(clusterProfiler)
gseaKEGG <- gseKEGG(geneList     = geneList,
                    organism     = 'hsa',
                    nPerm        = 1000,
                    minGSSize    = 20,
                    pvalueCutoff = 0.05,
                    verbose      = FALSE)
library(ggplot2)
dotplot(gseaKEGG,showCategory=4,split=".sign")+facet_grid(~.sign)

gseaKEGG_results <- gseaKEGG@result
library(enrichplot)
pathway.id = "hsa04060"
gseaplot2(gseaKEGG, 
          color = "red",
          geneSetID = pathway.id,
          pvalue_table = T)

library(pathview)
geneList_df <- data.frame(geneList) 
pv.out <- pathview(gene.data  = geneList,
                   pathway.id = "hsa05322",
                   species    = "hsa",
                   same.layer = T)
pv.out2 <- pathview(gene.data  = geneList,
                   pathway.id = "hsa05321",
                   species    = "hsa",
                   same.layer = T,
                   kegg.native = F)