riboseq的下游分析ribodiff，在R里进行GO分析和KEGG分析

1.导入数据

rm(list=ls())
options(stringsAsFactors = F)
library(tidyverse) # ggplot2 stringer dplyr tidyr readr purrr  tibble forcats
library(data.table) #多核读取文件
setwd("D:/87177/ngs/ADAR_PARCLIP_RIBOSEQ/ribodiff/quant")a1 <- fread('./result.txt',header = TRUE,data.table = F)

得到结果如上，所以需要调整row.names和最后一个col.names

fix(a1)
row.names(a1)<- a1$geneIDs####更换列名
counts <- a1[,2:ncol(a1)]###去除第一列

2.转换GENE ID和GENE NAME
设置环境，此处也包括GO的R包

library(clusterProfiler)
library(enrichplot)
library(tidyverse)
library(org.Hs.eg.db)###human
#library(org.Mm.eg.db)####mouse
library(DOSE)
library(pathview)

筛选需要的列

log2FC_cutoff = log2(10)
pvalue_cutoff = 0.001
padj_cutoff = 0.001
need_DEG <- counts[,c(2,3,6)]
head(need_DEG)
colnames(need_DEG) <- c('pvalue','padj','log2FoldChange')

得到结果如下：

转换ID

library(biomaRt)
library("curl")
my_mart <-useMart("ensembl")
datasets <- listDatasets(my_mart)
View(datasets)

选择数据集，并且将ID中的小数点去除

my_dataset <-useDataset("hsapiens_gene_ensembl",mart = my_mart)###选人类的数据集
inputdata<- as.data.frame(row.names(need_DEG))
newinput<- as.data.frame(gsub("\\.\\d*","",inputdata[,1]))#######去除小数点
colnames(newinput)<- "name"
row.names(need_DEG)<- newinput$name
####把rowname换掉

进行转换

listAttributes(my_dataset)
my_result <- getBM(attributes = c( "ensembl_transcript_id" , "external_gene_name" ,"description"),filters = "ensembl_transcript_id",values = newinput,mart = my_dataset)

得到结果如下：

然后把得到的result里面的genename，merge到数据里面

###DEG_SEQ <- cbind(inputdata,DEG_DEseq2)####给添加一列跟rowname一样的列
a1$geneIDs <- newinput$name###把a1的那一列的小数点去掉
colnames(my_result)<-c("geneIDs","gene_name","description")###把两个文件的列名改成一样的
final_res <- merge(my_result,a1, by = intersect(names(my_result), names(a1)))####进行合并

完成这一步骤，此处提供其他几个转换geneid的方法如下：
方法1：
ensembl_id转换与gene symbol基因名去重复的两种方法

方法2：
R|clusterProfiler-富集分析

方法3：
AnnoProbe：一个可以快速下载并处理GEO数据的好包

方法4：
uniport ID 转换为 gene symbol（ID转换）

方法5：
本文所选取的biomaRT,此处选用并非最好，根据自己的需求可以选择所需的方法。

方法6： TCGA和GTEx数据转换
tinyarray

去除重复

gene_name <-final_res$gene_name
table(duplicated(gene_name))
FALSE   TRUE 29461 162928

方法1：aggregate

final_res <- aggregate(final_res, by=list(gene_name), FUN=mean)####去除重复

方法2：!duplicated

ids <- final_res[!duplicated(final_res$gene_name),] ###提取不重复的

筛选出protein coding 的部分

library(AnnoProbe)
IDs <- need_DEG$gene_name
ID_type = "SYMBOL"
biotype <-annoGene(IDs,ID_type)
biotype <- subset(biotype,biotype$biotypes=="protein_coding")

进行下一步分析

row.names(ids)<-ids$gene_name
need_DEG <- ids[,c(5,6,9)]
head(need_DEG)
colnames(need_DEG) <- c('pvalue','padj','log2FoldChange')
log2FC_cutoff = log2(2)
pvalue_cutoff = 0.05
padj_cutoff = 0.05
gene_up=rownames(need_DEG[with(need_DEG,log2FoldChange>log2FC_cutoff & pvalue<pvalue_cutoff & padj<padj_cutoff),])
gene_down=rownames(need_DEG[with(need_DEG,log2FoldChange < -log2FC_cutoff & pvalue<pvalue_cutoff & padj<padj_cutoff),])
write.csv(gene_down,file = 'gene_down.csv')
write.csv(gene_up,file = 'gene_up.csv')####下载基因列表

得到的结果进行分析

gene_up_entrez <- as.character(na.omit(bitr(gene_up, #数据集fromType="ENSEMBLTRANS", #输入格式toType="ENTREZID", # 转为ENTERZID格式OrgDb="org.Hs.eg.db")[,2])) #"org.Hs.eg.db" "org.Mm.eg.db"
gene_down_entrez <- as.character(na.omit(bitr(gene_down, #数据集fromType="ENSEMBLTRANS", #输入格式toType="ENTREZID", # 转为ENTERZID格式OrgDb="org.Hs.eg.db")[,2])) #"org.Hs.eg.db" "org.Mm.eg.db"
gene_diff_entrez <- unique(c(gene_up_entrez ,gene_down_entrez ))#####kegg分析install.packages("R.utils")R.utils::setOption("clusterProfiler.download.method",'auto')
kegg_enrich_results <- enrichKEGG(gene  = gene_up_entrez,organism  = "hsa", #物种人类 hsa 小鼠mmupvalueCutoff = 0.05,qvalueCutoff = 0.2)kegg_enrich_results <- DOSE::setReadable(kegg_enrich_results, OrgDb="org.Hs.eg.db", keyType='ENTREZID')write.csv(kegg_enrich_results@result,'KEGG_gene_up_enrichresults.csv')
save(kegg_enrich_results, file ='KEGG_gene_up_enrichresults.Rdata')go_enrich_results <- enrichGO(gene = gene_up_entrez,OrgDb = "org.Hs.eg.db",ont   = "ALL"  ,     #One of "BP", "MF"  "CC"  "ALL" pvalueCutoff  = 0.05,qvalueCutoff  = 0.2,readable      = TRUE)write.csv(go_enrich_results@result, 'GO_gene_up_BP_enrichresults.csv')
save(go_enrich_results, file ='GO_gene_up_enrichresults.Rdata')

可视化画图

if (!require("BiocManager", quietly = TRUE))install.packages("BiocManager")BiocManager::install("AnnotationHub")
library(AnnotationHub)
genelist <- as.numeric(ids[,9])
names(genelist) <- row.names(ids)
genelist_entrez <- genelist
names(genelist_entrez) <- bitr(names(genelist),fromType="SYMBOL",toType="ENTREZID", OrgDb="org.Hs.eg.db")[,2]
genelist_entrez <- genelist_entrez[!is.na(names(genelist_entrez))]
kegg_enrich_results@result$Description[1:10]
i=3 select_pathway <-  kegg_enrich_results@result$ID[i]
pathview(gene.data     = genelist_entrez,pathway.id    = select_pathway,species       = 'hsa' ,      # 人类hsa 小鼠mmu kegg.native   = T,# TRUE输出完整pathway的png文件，F输出基因列表的pdf文件 new.signature = F, #pdf是否显示pathway标注limit         = list(gene=2.5, cpd=1) #图例color bar范围调整 )
pathview(gene.data     = genelist_entrez,pathway.id    = select_pathway,species       = 'hsa' ,      # 人类hsa 小鼠mmu kegg.native   = F,# TRUE输出完整pathway的png文件，F输出基因列表的pdf文件 new.signature = F, #pdf是否显示pathway标注limit         = list(gene=2.5, cpd=1)) #图例color bar范围调整

GO的有无环向图

go_enrich_results <- enrichGO(gene = gene_up_entrez,OrgDb = "org.Hs.eg.db",ont   = "BP"  ,     #One of "BP", "MF"  "CC"  "ALL" pvalueCutoff  = 0.05,qvalueCutoff  = 0.2,readable      = TRUE)####前面生成过一个result，但是当时ont的参数选择的是all，此处画图需要参数是bp
gop <- goplot(go_enrich_results, showCategory = 10)
ggsave(gop, filename = "goplot.pdf", width=10, height=10)

barplot与dotplot——最常用的可视化图形

### barplot
barp <- barplot(go_enrich_results, font.size=14, showCategory=10)+theme(plot.margin=unit(c(1,1,1,1),'lines'))
#如果enrichGO的ont为'ALL'
if (F) {barp <- barplot(go_enrich_results, split= "ONTOLOGY", font.size =14)+facet_grid(ONTOLOGY~., scale="free")+     theme(plot.margin=unit(c(1,1,1,1),'lines')) #调整图形四周留白大小}
ggsave(barp,filename = paste0(F,'.pdf'), width=10, height=10)### dotplot
dotp <- enrichplot::dotplot(go_enrich_results,font.size =14)+theme(legend.key.size = unit(10, "pt"),#调整图例大小plot.margin=unit(c(1,1,1,1),'lines'))#调整四周留白大小
#如果enrichGO的ont为'ALL'
if (F) {dotp <- enrichplot::dotplot(go_enrich_results,font.size =8,split = 'ONTOLOGY')+facet_grid(ONTOLOGY~., scale="free")+     theme(legend.key.size = unit(10, "pt"),#调整图例大小plot.margin=unit(c(1,1,1,1),'lines'))#调整四周留白大小}
ggsave(dotp,filename = paste0(T,'.pdf'),width =10,height =10)

cnetplot——Gene-Concept Network

### cnetplot: Gene-Concept Network
#构建含log2FC信息的genelist
genelist <- as.numeric(DEG_DEseq2[,2])
names(genelist) <- row.names(DEG_DEseq2)
###缺两个包
BiocManager::install("FlowSOM")
library(clusterProfiler)
library(ggnewscale)cnetp1 <- cnetplot(go_enrich_results,  foldChange = genelist,showCategory = 6,colorEdge = T,node_label = 'all',color_category ='steelblue')
cnetp2 <- cnetplot(go_enrich_results,  foldChange = genelist,showCategory = 6,node_label = 'gene',circular = T, colorEdge = T)
ggsave(cnetp1,filename ='cnetplot.pdf', width =12,height =10)
ggsave(cnetp2,filename = 'cnetplot_cir.pdf', width =15,height =10)

emapplot ——Enrichment Map

pt <- pairwise_termsim(go_enrich_results)
emapp <- emapplot(pt,showCategory = 30, node_label   = 'category')  # 'category', 'group', 'all' and 'none'.)
ggsave(emapp,filename = 'emapplot.pdf',width =10,height =10)

riboseq的下游分析ribodiff，在R里进行GO分析和KEGG分析相关推荐

使用R语言ggplot2包绘制pathway富集分析气泡图（Bubble图）:数据结构及代码
气泡图是在笛卡尔坐标系同加入大小的参数所形成的可以表示三个变量关系的图例.在对基因完成GO/KEGG分析后,使用气泡图可以直观的展示pathway.pvalue.count之间的关系.下面为使用R语言 ...
R 语言实现股票数据的预处理及分析
基于 R 语言的股票数据分析一.实验介绍 1.1 实验内容本实验是以股票数据作为分析背景,股票数据如何从雅虎财经板块上获取,观察股票每日价格和成交量数据开始,接着计算某一支股票数据中比较重要的日度 ...
R语言淮河流域水库水质数据相关性分析、地理可视化、广义相加模型GAM调查报告
最近我们被客户要求撰写关于水质数据的研究报告,包括一些图形和统计输出. 调查时间和地点采样地点:淮河流域一带,昭平台水库.白龟山水库.燕山水库.石漫滩水库.板桥水库.宿鸭湖水库.博山水库.南湾水库. ...
用plink做GWAS（PCA、关联分析）并用R绘图
用plink做GWAS(PCA.关联分析)并用R绘图 GWAS 一.观察初始数据二.质量控制样本缺失率和位点缺失率过滤(产生.imiss和lmiss文件) 计算MAF 数据清理三.检查亲缘关系 ...
R语言淮河流域水库水质数据相关性分析、地理可视化、广义相加模型GAM调查报告...
采样地点:淮河流域一带,昭平台水库.白龟山水库.燕山水库.石漫滩水库.板桥水库.宿鸭湖水库.博山水库.南湾水库.石山口水库.五岳水库.泼河水库.鲶鱼山水库(点击文末"阅读原文"获取 ...
R语言使用Rtsne包进行TSNE分析：通过数据类型筛选数值数据、scale函数进行数据标准化缩放、提取TSNE分析结果合并到原dataframe中（tSNE with Rtsne package）
R语言使用Rtsne包进行TSNE分析:通过数据类型筛选数值数据.scale函数进行数据标准化缩放.提取TSNE分析结果合并到原dataframe中(tSNE with Rtsne package) ...
R语言使用Rtsne包进行TSNE分析：提取TSNE分析结果合并到原dataframe中、可视化tsne降维的结果、并圈定降维后不匹配的数据簇（tSNE identifying mismatch)
R语言使用Rtsne包进行TSNE分析:提取TSNE分析结果合并到原dataframe中.可视化tsne降维的结果.并使用两个分类变量从颜色.形状两个角度来可视化tsne降维的效果.并圈定降维后不匹配 ...
R语言诊断试验数据处理与ROC分析实战案例2
R语言诊断试验数据处理与ROC分析实战案例2 目录 R语言诊断试验数据处理与ROC分析实战案例2 #ROC指标 #样例数据
R语言诊断试验数据处理与ROC分析实战案例1
R语言诊断试验数据处理与ROC分析实战案例1 目录 R语言诊断试验数据处理与ROC分析实战案例1 #ROC指标 #样例数据

riboseq的下游分析ribodiff，在R里进行GO分析和KEGG分析

riboseq的下游分析ribodiff，在R里进行GO分析和KEGG分析相关推荐

最新文章

热门文章