模式植物背景基因集制作

一边学习,一边总结,一边分享!

写在前面

关于GO背景基因集文件的制作,我们在很早以前也发过。近两天,自己在分析时候,也是被搞了头疼。想重新制作一份GO背景基因集,进行富集分析。但是结果,也不如意。以及在制作的过程中,也是跟随着以前的教程制作,发现以前整理的教程比较乱,那么借此机会,也进行整理,重新进行记录。

本期教程

直接访问链接:https://mp.weixin.qq.com/s/08hAZs24mi_KBOa4QZRLdQ

前言

我们在做转录组数据分析时,大多数都会进行功能富集分析,预测目的基因所具有的的功能。富集工具常用到的R语言中clusterProfiler包,里面包含了上千个功能富集的背景数据文件,功能非常强大,目前已经更新到V4.0版本。

在agriGO数据库中下载

网址:http://systemsbiology.cau.edu.cn/agriGOv2/index.php

前期准备文件

  1. 所需注释的物种基因核酸序列或蛋白序列
  2. swissprot数据
  3. idmapping.tb.gz文件
  4. go-basic.obo文件

数据下载

你可以分别进去对应的网址下载最新的数据库即可。

  1. Swissprot数据库:https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/
  2. dimapping数据:https://ftp.proteininformationresource.org/databases/idmapping/
wget -o GO_database/swissprot.gz  https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/swissprot.gz
wget -o GO_database/go-basic.obo http://purl.obolibrary.org/obo/go/go-basic.obo
wget -o GO_database/idmapping.tb.gz https://ftp.proteininformationresource.org/databases/idmapping/idmapping.tb.gz

建库及文件提取

1. 使用diamond makedb建库

diamond makedb --in GO_database/swissprot.gz --threads 60 --db GO_database/swissprot

2. GO号与swissprot蛋白ID文件的提取

下载idmapping数据库

https://ftp.proteininformationresource.org/databases/idmapping/idmapping.tb.gz

解压idmapping.tb.gz文件

gunzip idmapping.tb.gz


该文件的第一列为数据库ID,第八列为GO_ID,是我们这一次要与未知的结果进行转换的关键部分。

提取idmapping.tb.gz文件ID~GO.list file

awk -v FS="\t" -v OFS="\t" '{print $1,$8}' idmapping.tb | grep "GO" > idmapping.GO.list

使用Python脚本提取

输出结果

3. GO term文件提取

下载go-basi.obo,GO_Term

http://purl.obolibrary.org/obo/go/go-basic.obo

原始go-basi.obo文件格式

脚本提取


输出结果


基因文件序列的准备

下载所需的背景基因序列,核酸序列或蛋白序列度都可以

我们这里以番茄基因组4.0版本为例子。

#!/bin/bash
## download the tomato reference geneome to 4.1 
wget https://solgenomics.net/ftp//tomato_genome/annotation/ITAG4.0_release/ITAG4.0_gene_models.gffwget https://solgenomics.net/ftp//tomato_genome/assembly/build_4.00/S_lycopersicum_chromosomes.4.00.fagffread ITAG4.0_gene_models.gff -T -o ITAG4.0_gene_models.gtf
##
gffread ITAG4.0_gene_models.gff -T -o ITAG4.0_gene_models.gtf
##
gffread -w tomato_4.0.fa -g S_lycopersicum_chromosomes.4.00.fa ITAG4.0_gene_models.gtf

注意:若你不想这用操作,下载蛋白序列即可

1. 比对

diamond blastx -d GO_database/swissprot.dmnd -q ../Tomato_4.0/tomato_4.00.fa -k 1 -e 0.00001 -o tomato.gene.m8

2. 筛选出最佳结果

这步,若你认为有必要进行,那就进行筛选。筛选的参数可以自己调整。

使用*.pl教程

die "perl $0 *.m8 *.m8.out\n" if(@ARGV!=2);
open IN, "$ARGV[0]" or die "can not open file: $ARGV[0]\n";
open OA, ">$ARGV[1]" or die "can not open file: $ARGV[1]\n";my ($line,@inf,%score_data,%m8_data,%order);
my $n=1;
while($line=<IN>){chomp $line;@inf=split /\t/,$line;if($inf[11]>$score_data{$inf[0]}){$score_data{$inf[0]}=$inf[11];$m8_data{$inf[0]}=$line;}else{next;}       $order{$line}=$n++;
}
foreach my $i (sort {$order{$a}<=>$order{$b}} keys %order){@inf=split /\t/,$i;if(exists $m8_data{$inf[0]}){print OA "$m8_data{$inf[0]}\n";}
}
close IN;
close OA;

运行

perl m8_best_pick.pl tomato.gene.m8 tomato.gene.m8.best.out

3. 提取最佳结果ID文件

使用Python脚本:

或,你可以使用wak命令提取就可以。

运行

python ../get_blastx_wiss_id.py 02.tomato.gene.best.m8 > 03.tomato.transcript.swissprot.list

结果文件:

4. 合并文件,获得目标基因-GO ID


运行:

python get_go_annotation.py GO_batabase/idmapping.GO.list 03.tomato.transcript.swissprot.list

结果文件:

5. 拆分文件


注意:
我们的基因ID中,没有以mRNA:Solyc00g500003.1.1命名,如Solyc00g500003.1.1.我们需要将split_with_go.py进行适当修改即可。

运行:

python ../split_with_go.py 04.tomato_go.annotation  05.tomato.4.0.Go.list

结果文件:


到这里基本结束了,你获得Gene ID与对应GO ID

富集分析

你可以使用相关的云平台做GO功能富集分析,例如使用基迪奥生信平台GO功能富集工具

在线网址:OmicShare Tools - 基迪奥生信云工具:

上传背景基因

云平台支持的背景文件的数据格式

<,>

自己进行GO term的提取

下载go-basi.obo,GO_Term

http://purl.obolibrary.org/obo/go/go-basic.obo

原始go-basi.obo文件格式

Python脚本:

运行:

python get_go_term.py go-basic.obo

使用R进行合并

library(clusterProfiler)
## 加载背景基因文件“gene-GO"
go_anno <- read.delim('tomato_go.annotation.new', header = FALSE, stringsAsFactors = FALSE)
names(go_anno) <- c("gene_id", "GO_ID")
head(go_anno)### 导入GO注释文件
go_class <- read.delim("go_term.list", header = F, stringsAsFactors = F)
names(go_class) <- c("GO_ID", "Description","Ontology")
head(go_class)## 合并背景基因
go_ann <- merge(go_anno, go_class, by = 'GO_ID', all.x = F)
head(go_ann)


开始富集分析:

# 导入差异基因
gene_list <- read.table("tomato.gene.5000.txt", stringsAsFactors = F)
head(gene_list)
names(gene_list) <- c("gene_id")
gene_select <- gene_list$gene_id## 富集分析
go_rich <- enricher(gene = gene_select,TERM2GENE = go_ann[c('GO_ID','gene_id')],TERM2NAME = go_ann[c('GO_ID','Description')],pvalueCutoff = 0.05,pAdjustMethod = 'BH',qvalueCutoff = 0.2,maxGSSize = 200) 
head(go_rich)

**柱状图**
barplot(go_rich,drop=T,showCategory = 10) 

**气泡图**
dotplot(go_rich)


网络图

enrichplot::cnetplot(go_rich,circular = F, colorEdge = T)


写在最后,为了方便,我将前面的步骤进行分别写在一个脚本中。只要前期的数据准备好,输入所需的物种序列的序列即可。

算是比较方便。

准备文件GO_database

  1. swissprot.gz
  2. go-basic.obo
  3. idmapping.tb

运行脚本:

sh 01.run.swissprot.sh

结果文件:

  1. go_term.list
  2. idmapping.GO.list

GO注释文件脚本:

sh 02_run.GO_enrichment_file.sh test.fa
  • test.fa为注释文件序列

**注意:**若你不更改blast的脚本,这里默认只支持核酸序列。

结果文件:

在结果文件中05_gene.GO.list即最终结果文件。

后面的分析与前面的一致。


若你不想制作,我们这里提供完整的GO_database文件夹中的文件。你只需要在此基础上,运行你所需的物种序列即可。

直接访问链接:https://mp.weixin.qq.com/s/08hAZs24mi_KBOa4QZRLdQ

往期文章:

1. 复现SCI文章系列专栏

2. 《生信知识库订阅须知》,同步更新,易于搜索与管理。


3. 最全WGCNA教程(替换数据即可出全部结果与图形)

  • WGCNA分析 | 全流程分析代码 | 代码一

  • WGCNA分析 | 全流程分析代码 | 代码二

  • WGCNA分析 | 全流程代码分享 | 代码三


4. 精美图形绘制教程

  • 精美图形绘制教程

5. 转录组分析教程

转录组上游分析教程[零基础]


小杜的生信筆記,主要发表或收录生物信息学的教程,以及基于R的分析和可视化(包括数据分析,图形绘制等);分享感兴趣的文献和学习资料!!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/161229.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

排序【七大排序】

文章目录 1. 排序的概念及引用1.1 排序的概念1.2 常见的排序算法 2. 常见排序算法的实现2.1 插入排序2.1.1基本思想&#xff1a;2.1.2 直接插入排序2.1.3 希尔排序( 缩小增量排序 ) 2.2 选择排序2.2.1基本思想&#xff1a;2.2.2 直接选择排序:2.2.3 堆排序 2.3 交换排序2.3.1冒…

uniapp小程序中给web-view页面添加授权弹窗(使用cover-view组件覆盖实现该功能)

效果图&#xff1a; web-view是承载网页的容器。会自动铺满整个小程序页面&#xff0c;个人类型的小程序暂不支持使用。 再看下面一个提示&#xff1a; 每个页面只能有一个 web-view&#xff0c;web-view 会自动铺满整个页面&#xff0c;并覆盖其他组件。 也就是说&#xff0c;…

云安全—云计算基础

0x00 前言 学习云安全&#xff0c;那么必然要对云计算相关的内容进行学习和了解&#xff0c;所以云安全会分为两个部分来进行&#xff0c;首先是云计算先关的内容。 0x01 云计算 广泛传播 云计算最早大范围传播是2006年&#xff0c;8月&#xff0c;在圣何塞【1】举办的SES&a…

七大排序 (9000字详解直接插入排序,希尔排序,选择排序,堆排序,冒泡排序,快速排序,归并排序)

一&#xff1a;排序的概念及引入 1.1 排序的概念 1.1 排序的概念 排序&#xff1a;所谓排序&#xff0c;就是使一串记录&#xff0c;按照其中的某个或某些关键字的大小&#xff0c;递增或递减的排列起来的操作。 稳定性&#xff1a;假定在待排序的记录序列中&#xff0c;存在…

DITA-OT 4.0新特性 - PDF themes,定制PDF样式的新方法

随着DITA-OT 4.0的发布&#xff0c;它提供了一种新的定制PDF样式方法&#xff0c;这种方法就是PDF theme。这篇文章来聊一聊这种定制PDF输出的新方法和实验结果。 在进入PDF theme细节之前&#xff0c;为各位读者梳理一下DITA-OT将DITA和Markdown发布成PDF的几种方法。 - 1 …

element ui 下拉框 选择月份和天数

一、背景 目前做的管理系统项目&#xff0c;期望实现功能为&#xff1a;设置出账周期和出账日&#xff0c;考虑使用element ui下拉框实现功能 二、所用技术 vue2element ui 三、实现效果 四、具体代码 <template><popup-frame :title"批量设置出账日" …

OJ第四篇

文章目录 链表分割环形链表有效的括号 链表分割 链接: 链表分割 虽然这个题牛客网中只有C,但是无所谓&#xff0c;我们只要知道C是兼容C的就可以了 至于说这个题的思路&#xff0c;我们就弄两个链表&#xff0c;把小于x的结点放到一个链表中&#xff0c;剩下的放到另一个链表…

2023年10月工作经验及问题整理总结

目录 1.window自带的base64加密解密 2.ElementUI修改鼠标移动到表格的背景色 3.vscode保存时几万个eslint错误 4.Git 拉取Gitee仓库报错&#xff1a;“fatal: unable to access ": Failed to connect to 127.0.0.1 port 1080: Connection r... 4.1本地查看Git是否使用…

Java版本spring cloud + spring boot企业电子招投标系统源代码

项目说明 随着公司的快速发展&#xff0c;企业人员和经营规模不断壮大&#xff0c;公司对内部招采管理的提升提出了更高的要求。在企业里建立一个公平、公开、公正的采购环境&#xff0c;最大限度控制采购成本至关重要。符合国家电子招投标法律法规及相关规范&#xff0c;以及审…

spark stream入门案例:netcat准实时处理wordCount(scala 编程)

目录 案例需求 代码 结果 解析 案例需求&#xff1a; 使用netcat工具向9999端口不断的发送数据&#xff0c;通过SparkStreaming读取端口数据并统计不同单词出现的次数 -- 1. Spark从socket中获取数据&#xff1a;一行一行的获取 -- 2. Driver程序执行时&#xff0c…

加固数据安全:Java助力保护Excel文件,让数据无懈可击

摘要&#xff1a;本文由葡萄城技术团队于CSDN原创并首发。转载请注明出处&#xff1a;葡萄城官网&#xff0c;葡萄城为开发者提供专业的开发工具、解决方案和服务&#xff0c;赋能开发者。 前言 Excel文件保护是常用的一种功能&#xff0c;文件保护主要有三种&#xff1a; 添…

单目3D自动标注

这里介绍两种 1. 基于SAM的点云标注 Seal&#xff1a;是一个多功能的自监督学习框架&#xff0c;能够通过利用视觉基础模型的现成知识和2D-3D的时空约束分割自动驾驶数据集点云 Scalability&#xff1a;可拓展性强&#xff0c;视觉基础模型蒸馏到点云中&#xff0c;避免2D和…

SaaS人力资源管理系统的Bug

SaaS人力资源管理系统的Bug Bug1【18】 这里我是直接把代码复制过来的&#xff0c;然后就有一个空白 这是因为它的代码有问题&#xff0c;原本的代码如下所示 <el-table-column fixed type"index" label"序号" width"50"></el-table…

Linux性能优化--性能追踪:受CPU限制的应用程序(GIMP)

10.0 概述 本章包含了一个例子&#xff1a;如何用Linux性能工具在受CPU限制的应用程序中寻找并修复性能问题。 阅读本章后&#xff0c;你将能够&#xff1a; 在受CPU限制的应用程序中明确所有的CPU被哪些源代码行使用。用1trace和oprofile弄清楚应用程序调用各种内部与外部函…

KNN-近邻算法 及 模型的选择与调优(facebook签到地点预测)

什么是K-近邻算法&#xff08;K Nearest Neighbors&#xff09; 1、K-近邻算法(KNN) 1.1 定义 如果一个样本在特征空间中的k个最相似(即特征空间中最邻近)的样本中的大多数属于某一个类别&#xff0c;则该样本也属于这个类别。 来源&#xff1a;KNN算法最早是由Cover和Hart提…

【JVM面试】从JDK7 到 JDK8, JVM为啥用元空间替换永久代?

系列文章目录 【JVM系列】第一章 运行时数据区 【面试】第二章 从JDK7 到 JDK8, JVM为啥用元空间替换永久代&#xff1f; 大家好&#xff0c;我是青花。拥有多项发明专利&#xff08;都是关于商品、广告等推荐产品&#xff09;。对广告、Web全栈以及Java生态微服务拥有自己独到…

VSS、VDD、VBAT、VSSA

引言 在学习设计TM32时&#xff0c;发现芯片除了GPIO引脚外还会引出许多引脚&#xff0c;以STM32F407ZGT6为例除了GPIO引脚还会有以下引脚 如VSS、VDD、VBAT、VSSA、NRST、VREF、VDDA、VCAP_1、VCAP_2、PDR_ON这些引脚。他们有何作用&#xff0c;电路设计中应如何连接&#x…

迅为RK3588开发板使用RKNN-Toolkit-lite2运行测试程序

1 首先也需要部署运行环境&#xff0c;将库文件放入 RK3588 开发板上&#xff0c;我们将网盘资料“iTOP-3588 开发 板 \02_ 【 iTOP-RK3588 开 发 板 】 开 发 资 料 \12_NPU 使 用 配 套 资 料 \05_Linux_librknn_api\librknn_api\aarch64”路径下的文件通过U盘拷贝到开发板的…

LLM-RAG-WEB 大模型+文件+可视化界面

注意&#xff1a;这里只是简单实现了功能和界面&#xff0c;文件对话也暂时只支持一个文件&#xff0c;如果跳到模型对话再切换回文件对话会文件会删除重置会话&#xff0c;但模型对话切换回来时保留之前会话的。 1、代码&#xff08;使用步骤说明在链接里&#xff09; 参考下…

SQLite4Unity3d安卓 在手机上创建sqlite失败解决

总结 要在Unity上运行一次出现库&#xff0c;再打包进APK内 问题 使用示例代码的创建库 var dbPath string.Format("Assets/StreamingAssets/{0}", DatabaseName); #else// check if file exists in Application.persistentDataPathvar filepath string.Format…