win10下的RNA测序（二）

写在前面

这部分讲的是根据上游得到的counts数做差异分析以及通路富集

读入文件

options(stringsAsFactors = F)
rm(list = ls())
pwd = "E:\\process\\"
setwd(pwd)
library(readr)
library(dplyr)
data = read.csv("exp.txt",header = T,sep = "\t",
                fill=T,stringsAsFactors=F,skip = 1)
data = data[,c(1,7:15)]# 提取探针以及bam文件对应的counts数

先将探针转换为基因ID

这边转化ID，刚好之前写了个教程https://editor.csdn.net/md/?articleId=115487976，本着用代码处理的原则，则尝试了新方法

library("clusterProfiler")
library(dplyr)
gene = data$Geneid# %>% as.data.frame()
df = data.frame(gene)
#BiocManager::install("AnnotationDbi")
library("AnnotationDbi")
library("org.Mm.eg.db")
df$symbol <- mapIds(org.Mm.eg.db,
                    keys=gene,
                    column="SYMBOL",
                    keytype="ENSEMBL",
                    multiVals="first")
df = na.omit(df)

提取有对应基因的表达谱

如果一个基因对应多个探针的话，取最大值

## 2.1 subset the probe correspondant  to gene_symbols 
data = subset(data,Geneid%in%df$gene)
rownames(data) = data$Geneid 
df = df[match(df$gene,rownames(data)),]
data = data[,-1]
## 2.2 a gene correspondant different probe,you can choose mean or max 
data_finally = apply(data,2,function(x){
  tapply(x,df$symbol,max)
})

差异分析

# 加载包
library(edgeR)
library(limma)

# 设定组别
group_list = c(rep("PBS",3),rep("yangxing",3),rep("yinxing",3))
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))

## 3.1 filter some low express counts and transform to log2
dge <- DGEList(counts=data_finally)
keep <- filterByExpr(dge, design)
dge <- dge[keep,,keep.lib.sizes=FALSE]
dge <- calcNormFactors(dge)
logCPM <- cpm(dge, log=TRUE, prior.count=3)
rownames(design)=colnames(dge)

## 3.2 construct contrast matrix and search degs  
contrast_class=c("yangxing-PBS","yinxing-PBS","yangxing-yinxing")

## 3.3 define function to search degs 
search_limma_degs=function(data_subset,group_list,contrast_class){
    design <- model.matrix(~0+factor(group_list))
    colnames(design)=levels(factor(group_list))
    head(design)
    exprSet=as.matrix(data_subset)
    mode(exprSet)="numeric"
    rownames(design)=colnames(exprSet)
    head(design)
    
    # 2.2 contrast.matrix
    
    contrast.matrix<-makeContrasts(contrasts = contrast_class,
                                   levels = design)
    
    # 3.search the degs 
    # 3.1 degs 
    ##step1,lmfit起过滤的手续
    fit <- lmFit(exprSet,design)
    ##step2
    fit2 <- contrasts.fit(fit, contrast.matrix)
    fit2 <- eBayes(fit2)  ## default no trend !!!
    ##eBayes() with trend=TRUE
    
    data_degs=list()
    for (i in 1:length(contrast_class)){
      print(i)
      tempOutput1 = topTable(fit2, coef=i, number =Inf)
      data_degs[[i]] = na.omit(tempOutput1) 
    }
    return(data_degs)
  }
# search degs 
data_degs=search_limma_degs(data_subset = logCPM,group_list,contrast_class)
deg_yangxing_PBS = data_degs[[1]]
deg_yinxing_PBS = data_degs[[2]]
deg_yangxing_yinxing = data_degs[[3]]

通路富集与火山图以及相应图表的保存

valcano_pathway = function(nrDEG,i,contrast_class){
df=nrDEG
attach(nrDEG)
df$y= -log10(adj.P.Val)
df$state=ifelse(df$adj.P.Val>0.05,'NO',
                ifelse( df$logFC >2,'UP',
                        ifelse( df$logFC < -2,'DOWN','NO') )
)
a = table(df$state)
df$state1<-factor(df$state,
                      levels = names(a),
                      labels = paste(names(a),as.numeric(a)))
df$name=rownames(df)
head(df)
library(ggplot2)
library(ggpubr)

# valcano 
ggplot(df,aes(x=logFC,y=y))+
    geom_point(aes(color=state1))+
    labs(x="log2FC",y="-log10FDR")+
    #增加阈值线:分别对应FDR=0.05,|log2FC|=1
    geom_hline(yintercept=-log10(0.05),linetype=4)+
    geom_vline(xintercept=c(-2,2),linetype=4)+
    theme_classic() + ggtitle(contrast_class[i]) +
    theme(plot.title = element_text(hjust = 0.5)) +
  labs(color="Pvalue<0.05\n|log2FC|>2")
ggsave(filename = paste(contrast_class[i],"valcano.png"))
write.csv(df,file = paste(contrast_class[i],".csv"))
# pathway
if (T) {
  df$ENTREZID <- bitr(df$name, fromType = "SYMBOL",
             toType = c( "ENTREZID"),
             OrgDb = org.Mm.eg.db)$ENTREZID
  gene_up= df[df$state == 'UP','ENTREZID'] 
  gene_down=df[df$state == 'DOWN','ENTREZID'] 
  ##### up 
  
  go <- enrichGO(gene_up, OrgDb = "org.Mm.eg.db", ont="all") 
  library(ggplot2)
  library(stringr)
  barplot(go, split="ONTOLOGY",font.size =10)+ 
    facet_grid(ONTOLOGY~., scale="free") + 
    scale_x_discrete(labels=function(x) str_wrap(x, width=50))+
    ggsave(filename = paste(contrast_class[i],".gene_up_GO_all_barplot.png"))
  
  
  enrichKK <- enrichKEGG(gene         =  gene_up,
                         organism     = 'mmu',
                         pvalueCutoff = 0.1,
                         qvalueCutoff =0.1)
  enrichKK=DOSE::setReadable(enrichKK, OrgDb='org.Mm.eg.db',keyType='ENTREZID')
  barplot(enrichKK,showCategory=20)
  ggsave(filename = paste(contrast_class[i],".gene_up_kegg_all_barplot.png"))
  dotplot(enrichKK)
  ggsave(filename = paste(contrast_class[i],".gene_up_kegg_all_dotplot.png"))
  write.csv(go@result,file = paste(contrast_class[i],"gene_up_GO.csv"))
  write.csv(enrichKK@result,file = paste(contrast_class[i],"gene_up_KEGG.csv"))
  
  ###### down 
  go <- enrichGO(gene_down, OrgDb = "org.Mm.eg.db", ont="all") 
  library(ggplot2)
  library(stringr)
  barplot(go, split="ONTOLOGY",font.size =10)+ 
    facet_grid(ONTOLOGY~., scale="free") + 
    scale_x_discrete(labels=function(x) str_wrap(x, width=50))+
    ggsave(filename = paste(contrast_class[i],".gene_down_GO_all_barplot.png"))
  
  enrichKK <- enrichKEGG(gene         =  gene_down,
                         organism     = 'mmu',
                         pvalueCutoff = 0.1,
                         qvalueCutoff =0.1)
  enrichKK=DOSE::setReadable(enrichKK, OrgDb='org.Mm.eg.db',keyType='ENTREZID')
  barplot(enrichKK,showCategory=20)
  ggsave(filename = paste(contrast_class[i],".gene_down_kegg_all_barplot.png"))
  dotplot(enrichKK)
  ggsave(filename = paste(contrast_class[i],".gene_down_kegg_all_dotplot.png"))
  write.csv(go@result,file = paste(contrast_class[i],"gene_down_GO.csv"))
  write.csv(enrichKK@result,file = paste(contrast_class[i],"gene_down_KEGG.csv"))
  
  
  
}
}
valcano_pathway(nrDEG = deg_yangxing_PBS,1,contrast_class = contrast_class)
valcano_pathway(nrDEG = deg_yinxing_PBS,2,contrast_class = contrast_class)
valcano_pathway(nrDEG = deg_yangxing_yinxing,3,contrast_class = contrast_class)

结果展示

通路富集的另外食用方法

通路富集也可以通过kobas使用
http://kobas.cbi.pku.edu.cn/retrieve/visualization/?taskid=c74679b1b6e04162bf1dccbf43672018&app=gene_list使用