使用R语言的clusterProfiler对葡萄做GO富集分析的简单小例子

葡萄的参考基因组下载自NCBI，下载链接是https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/003/745/GCF_000003745.3_12X/

基本流程是
  • Hiast2 比对
  • samtools sam 转bam
  • stringtie组装转录本
  • gffcompare将stringtie输出的gtf文件与参考基因组的注释文件做比较得到一- 个merged.combine.gtf
  • 使用merged.combine.gtf 这个文件对每个样本计算表达量，输出文件存储到ballgown文件夹下，这一步用到的命令是 stringtie -e -B -p 8 -G merged.combined.gtf -o ballgown/L01/L01.gtf output_bam/L01.sorted.bam
image.png
  • 接下来是R语言的ballgown包读入数据获取基因和转录本的表达量代码是
library(ballgown)
library(genefilter)
library(dplyr)

pheno_data<-read.csv('pheno_data.txt',header=T)
grape_expr <- ballgown(dataDir = 'ballgown',
                    samplePattern = 'L0',
                    pData = pheno_data)

image.png
image.png

这一步可以拿到gene_id还有gene_name ,FPKM的表达量，cov对用的应该是reads count吧。

接下来是差异表达分析

代码是

grape_expr_filter<-subset(grape_expr,                          'rowVars(texpr(grape_expr))>1',                          genomesubset=T)results_genes <- stattest(grape_expr_filter,                          feature = 'gene',                          covariate = 'time_point',                          getFC = TRUE,                          meas = 'FPKM')#results_genes <- arrange(results_genes,pval)results_genes%>%  filter(abs(fc)>=2&pval<0.05) -> results_genes_diff

dim(results_genes_diff)head(results_genes_diff)

现在有了基因id

image.png

接下来是使用clusterProfiler做go注释

这里参考https://guangchuangyu.github.io/cn/2017/07/clusterprofiler-maize/#disqus_thread

首先把这个基因id对应着转换成 ENTREZID ，这里需要借助对应的gtf注释文件

这里只关注蛋白编码基因

grep 'gene_biotype 'protein_coding'' 12X_genomic.gtf > 12X_protein_coding.gtf
#python
from gtfparse import read_gtf
known_proteincoding = read_gtf('12X_protein_coding.gtf')
known_proteincoding.to_csv('all_protein_coding.csv')

GO富集分析的R语言代码

require(AnnotationHub)hub<-AnnotationHub() #这一步对网路有要求# aa<-query(hub,'zea')# aa$title# query(hub,'Malus domestica')bb<-query(hub,'Vitis vinifera')#bb$titlegrape<-hub[['AH85815']]# length(keys(grape))# columns(grape)protein_coding_all<-read.csv('all_protein_coding.csv',header=T)df<-merge(results_genes_diff,grape_expr_filter@expr$trans,by.x='id',by.y='gene_id')df1<-merge(df,protein_coding_all,by.x='gene_name',by.y='gene_id')dim(df1)gene_ids<-   df1$db_xref[!duplicated(df1$db_xref)]gene_ids<-stringr::str_replace(gene_ids,'GeneID:','')

library(clusterProfiler)bitr(keys(grape)[2],'ENTREZID',c('REFSEQ','GO',                                 'ONTOLOGY','GENENAME',                                 'SYMBOL'),grape)res = enrichGO(gene_ids,                OrgDb=grape, pvalueCutoff=0.05, qvalueCutoff=0.05)help(package='clusterProfiler')dotplot(res)

最后的结果是

image.png

欢迎大家关注我的公众号

小明的数据分析笔记本

小明的数据分析笔记本

数据分析和数据可视化有意思的简单小例子~石榴研究生的笔记本
351篇原创内容
公众号

小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！

YuLabSMU

南方医科大学生物信息学系余光创课题组
684篇原创内容
公众号
赞 (0)

相关推荐