尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

生物数据库与数据格式详解:从NCBI到FASTQ/GFF

生物数据库与数据格式详解:从NCBI到FASTQ/GFF 简介一套面向生物信息学初学者和高校教师的公开课获奖课件系统讲解常用生物数据库与数据格式解决入门时‘数据多、格式多、数据库多’难以下手的问题。课件先介绍生物数据库的整体背景再逐类拆解FastA、FastQ、GFF、GenBank等常见序列与注释格式的书写规范和适用场景并详细展示GenBank文件的LOCUS、DEFINITION、ACCESSION、FEATURES、ORIGIN等主要段落。随后的数据库模块覆盖NCBI、EBI、DDBJ三大常用序列数据库以及Gene Ontology、KEGG、InterPro、UCSC基因组浏览器、Ensembl等基因功能与基因组数据库便于读者理解不同数据库的定位和选用逻辑。课件对FastA描述行与序列行、FastQ四行质量记录、GFF注释列含义、GenBank核心字段均配有示例便于对照理解同时通过最新生物数据库列表和典型检索场景帮助学习者建立从格式识别到数据库选择的分析思路。整套内容共1个pptx文件压缩包约9.49MB下载后可直接用于课堂演示、备课或自学目前已有88人学习下载适合作为生信课程讲义、组会入门分享或考前复习的参考资料。1. 一份讲清生物数据库与数据格式的课件先看它能不能帮你少走弯路做生信的人最常遇到的不是算法不会写而是手里一堆数据文件不知道怎么归位从 NCBI 拉下来一个.gbff打开全是 LOCUS、FEATURES、ORIGIN不知道哪个字段对应基因注释从测序公司拿到.fastq四行一条记录里面那个 ASCII 质量值到底怎么换算想对比一下 UCSC 和 Ensembl 的基因注释又发现两个数据库的坐标习惯不一样。这份「常用生物数据库和数据格式」公开课获奖课件本质上是把数据库全景和数据格式这两条线合到一套 PPT 里讲——先讲序列数据存在哪里NCBI、EBI、DDBJ再讲注释信息去哪里查GO、KEGG、InterPro最后把 FASTA、FASTQ、GenBank、GFF 这几类最常见的格式逐字段拆开。它适合刚入门的生信学生快速建立「数据库—格式—工具」的对应关系也适合科研助理拿来做备课底稿或者新人培训材料属于那种看完就能对着 NCBI 页面操作的实用型课件。2. 从三大核酸库到功能注释库生物数据库的全景与选库逻辑拿到一个物种的序列先去哪个库查直接决定你后面注释工作的效率。课件把数据库分成序列库、功能库和基因组库三层这个分层本身就是一套选库逻辑。2.1 三大核酸序列数据库NCBI、EBI、DDBJ 的分工与同步NCBI美国国家生物技术信息中心、EBI欧洲生物信息学中心和 DDBJ日本 DNA 数据库是国际核酸序列数据库联盟的三方成员日常数据是每日同步的。也就是说绝大多数情况下你在 NCBI 查得到的序列在 EBI 和 DDBJ 也有对应条目三者不是竞争关系而是互为镜像。实际使用中我一般遵循「谁的下拉页面顺手就用谁」的原则NCBI 的 Nucleotide 页面适合按 accession 号精确检索EBI 的 ENA 适合批量下载DDBJ 的优势在于亚洲地区的提交和部分特有数据集。课件里反复强调一个概念这些数据库的核心单位是记录record不是文件。每条记录对应一个唯一的检索号accession number比如 U49845这是序列记录在数据库中的唯一指针不管你怎么下载、怎么转换格式这个号不能丢。实际操作中经常有人把文件名改了结果回头自己都分不清这条序列来自哪条记录这点后面避坑章节还会展开。2.2 功能注释数据库GO、KEGG、InterPro 各管一段序列库只回答「是什么序列」功能库回答「这个序列的产物干什么」。课件把功能注释进一步划成三块数据库全称回答的问题典型应用场景Gene OntologyGO基因本体数据库基因产物在分子功能、细胞组分、生物过程中的角色差异基因富集分析KEGG京都基因与基因组百科全书代谢通路、信号通路、酶与反应关系通路富集、代谢网络分析InterPro蛋白家族与结构域数据库蛋白序列属于哪个家族、含哪些结构域蛋白功能预测、结构域注释选库逻辑很简单你想知道基因参与什么通路KEGG 优先你想知道蛋白有没有跨膜结构域、锌指结构InterPro 优先你想做整体的功能分类和富集统计GO 是地基。课件里把 InterPro 单独列出来讲是有道理的——它在功能注释流程里经常被当成「黑匣子」其实它背后聚合了 Pfam、PANTHER、PRINTS 等十几个子库同一个蛋白在不同子库里比对结果不一致时InterPro 会给出综合判定这个综合判定的置信度比单个子库的结果更值得写进论文。2.3 基因组数据库UCSC 与 Ensembl 的互补关系基因组层面课件聚焦 UCSC Genome Browser 和 Ensembl 两个入口。两者的区别在于UCSC 以基因组浏览器为核心交互方式擅长可视化浏览和 Table Browser 批量下载Ensembl 以注释为核心对基因结构的注释体系更完整且提供 BioMart 数据挖掘接口。做比较基因组或者要看某个基因在多个物种里的同源关系Ensembl 的 Compara 管线是首选做变异位点可视化、看 ChIP-seq 峰图在基因组上的分布UCSC 更直接。需要特别注意的是这两个数据库对同一基因的坐标注释可能不一致。原因在于它们使用的参考基因组版本和转录本注释版本可能不同比如 UCSC 习惯用chr1作为染色体命名Ensembl 可能用1直接拿两个数据库的坐标混用是翻车高发区。课件里提到的最新生物数据库列表Nucleic Acids Research 每年发布实际上是一个很好的导航工具NAR 每年一月刊会汇总当年所有主流生物数据库的更新情况和访问入口按分子类型、物种、功能类型分门别类比搜索引擎更可靠。3. 把 FASTA、FASTQ、GBFF、GFF 读明白四种格式的字段与解析要点格式问题看起来简单实际上是生信入行最常见的「第一道坎」。课件把这四种格式摊开讲我按自己的使用经验补上解析层面的细节。3.1 FASTA描述行与序列体的基本规则FASTA 是应用范围最广的序列格式格式本身极其简单以开头的描述行后面跟着若干行序列字符。课件里强调的「描述行没有统一标准」这点很关键——你可以在后面写 accession 号也可以写物种名加基因名甚至写一段叙述文字但解析工具通常会取第一个空白字符前的部分作为序列 ID。# 统计一个 FASTA 文件里有多少条序列 grep -c ^ example.fasta # 提取每条序列的 ID 和长度用于快速检查 awk /^/{if(seq) print id, length(seq); id$1; sub(^, , id); seq} !/^/{seqseq$0} END{print id, length(seq)} example.fasta提示grep -c ^统计的是以开头的行数等于序列条数。awk脚本在遇到新的行时输出上一条序列的 ID 和长度这里$1取的是描述行第一个字段所以描述行里如果写的是sp|P12345|PROTEIN_NAME最终提取的 ID 是sp|P12345|PROTEIN_NAME去掉后的完整内容而不是截断到|之前——截断逻辑需要额外指定分隔符。FASTA 格式不限定扩展名.fa、.fasta、.fna核酸、.faa蛋白都常见解析时不依赖扩展名判断内容类型。3.2 FASTQ四行结构与质量值编码理解FASTQ 是测序下机的原始格式和 FASTA 最大的区别是多了质量值。每四条线构成一条记录第一行以开头是序列标识第二行是碱基序列第三行以开头一般重复标识或留空第四行是质量值字符串字符数量必须与第二行碱基数量完全一致。质量值编码是新手最容易懵的地方。第四行的每个 ASCII 字符对应一个 Phred 质量分数公式是Q ASCII码 - 偏移量。Illumina 1.8 之后的平台统一使用 Phred33 编码也就是 ASCII 33 对应 Q0老式的 Solexa 和早期 Illumina 用 Phred64。判断编码方式的一个办法是看质量值字符的范围如果出现!ASCII 33到5ASCII 53这种低位字符基本就是 Phred33如果字符普遍在BASCII 66以上大概率是 Phred64。# 用 Python 快速读一条 FASTQ 并检查质量值范围 from pathlib import Path def read_fastq_record(fh): # 读取四行构成一条完整记录 header fh.readline().strip() if not header: return None seq fh.readline().strip() plus fh.readline().strip() qual fh.readline().strip() return header, seq, qual with open(sample.fastq, r) as fh: header, seq, qual read_fastq_record(fh) ascii_codes [ord(ch) for ch in qual] print(f最低字符: {chr(min(ascii_codes))} (ASCII {min(ascii_codes)})) print(f最高字符: {chr(max(ascii_codes))} (ASCII {max(ascii_codes)}))提示这个检查脚本不重复造轮子它读入第一条记录后取质量值字符串中最小和最大的 ASCII 码。如果最小值在 33 附近说明是 Phred33如果最小值在 64 以上就要考虑 Phred64。实际生产中可以直接用fastqc软件包的报告替代这个脚本但在没有外部工具的环境里这个检查逻辑能应急。3.3 GenBank/GBFF三段式结构与 LOCUS 行逐字段拆解GenBank 格式文件扩展名通常为.gb或.gbff是一种带完整注释的序列格式比 FASTA 复杂得多。课件把它拆成三部分头部描述符、FEATURES注释特性和 ORIGIN序列本身每条记录以//结尾。这个三段式的理解非常关键——头部回答「这条记录是什么」FEATURES 回答「序列上有什么」ORIGIN 才是真正的序列内容。LOCUS 行是头部第一行包含六个关键字段课件里逐一强调了字段位置内容说明第一项LOCUS 名称仅作检索码使用无实际生物学意义第二项序列长度单条记录长度GenBank 规定不能超过 350kb第三项分子类型DNA、RNA、ss-DNA 等必须为单一分子类型第四项GenBank 分类码由 3 个字母组成如 PLN植物、BCT细菌、VRL病毒第五项拓扑结构linear线性或 circular环状第六项修订日期格式 DD-MMM-YYYY有时仅为数据首次公开日期FEATURES 部分用的是Location/Qualifiers语法。比如CDS 1..206表示编码区从第 1 个碱基之后某个位置开始表示不完整到第 206 个碱基结束。每个 CDS 特性下面会跟/codon_start、/product、/translation等限定符——/translation里就是蛋白序列这是后期提取 CDS 序列的重要来源。# 从 GenBank 文件里提取所有 CDS 的 product 名称和坐标 grep -E ^ {5}CDS|/product example.gb提示这里的grep正则利用了 GenBank 缩进规则——特性关键词CDS、gene、mRNA固定在五级缩进位置限定符/product前有 21 个空格。不同来源的文件缩进可能略有差异稳妥做法是用 BioPython 的SeqIO模块解析而不是纯文本匹配。纯grep只适合快速浏览不适合批量提取。3.4 GFF/GTF九列基本格式与坐标习惯GFFGeneral Feature Format是基因组注释的交换格式GFF3 是当前主流版本。每一行表示一个注释特性共九列序列名、来源、类型、起始、终止、得分、链、相位、属性。GFF3 和 GTF 的区别集中在第九列——GTF 用键值对如gene_id XXX;GFF3 用分号分隔的标签值对如IDgene1;Namexxx。坐标体系是这里最容易翻车的地方GFF 的起始和终止坐标是1-based 闭合区间也就是说一个从第 100 个碱基到第 200 个碱基的基因写成100 200包括这两个端点长度是 101。而 SAM/BAM 格式用的是 0-based 半开区间BED 也一样。做格式转换时这个差异是 bug 高发点后面避坑章节统一讲。# 提取 GFF3 中所有注释为 gene 的行只看前九列 awk -F\t $3gene {print $1, $4, $5, $7, $9} annotation.gff3提示awk用-F\t指定制表符分隔条件$3gene筛选第三列类型为 gene 的记录输出染色体、起始、终止、链和属性。需要注意的是 GFF 文件必须用制表符分隔不能用空格否则解析会错乱这也是很多人拿 Excel 改 GFF 后文件报废的原因。4. 从下载到落地NCBI、Ensembl、UCSC 的检索下载与格式转换实操光看懂格式不够得真正把数据拉下来。这一章把课件里提到的三大来源的常用下载路径走一遍。4.1 从 NCBI 下载指定 accession 的 GenBank 记录NCBI 的 Nucleotide 数据库支持直接按 accession 号访问最常见的方式是走 E-utilities 接口。以课件里的酵母基因记录 U49845 为例# 用 E-utilities 下载 GenBank 格式记录 curl https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?dbnucleotideidU49845rettypegbretmodetext -o U49845.gb # 下载 FASTA 格式 curl https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?dbnucleotideidU49845rettypefastaretmodetext -o U49845.fa提示dbnucleotide指定数据库为核酸序列库id是 accession 号rettype控制返回格式——gb是 GenBankfasta是 FASTAgbwithparts用于包含分段记录的完整 GenBank。retmodetext表示返回纯文本而不是 XML。批量下载多条序列时可以用逗号分隔多个 ID也可以把 ID 列表写入文件配合batch模式使用。E-utilities 有频率限制未加 API key 时每秒最多 3 次请求批量下载建议先注册 NCBI API key 并在请求中加上api_key参数。4.2 从 Ensembl 下载 GFF3 与配套 FASTAEnsembl 的下载路径比 NCBI 更规整所有物种的注释文件都在一个固定目录下以人类为例# 下载人类基因注释 GFF3Release 最新版本 wget https://ftp.ensembl.org/pub/current_gff3/homo_sapiens/Homo_sapiens.GRCh38. current.gff3.gz # 下载配套的 cDNA 序列 FASTA wget https://ftp.ensembl.org/pub/current_fasta/homo_sapiens/cdna/ Homo_sapiens.GRCh38.cdna.all.fa.gz提示下载后先检查current目录下实际的 release 版本号不同版本对应不同的参考基因组版本。GFF3 和 FASTA 必须使用相同版本否则会出现转录本 ID 对不上序列的情况。解压后建议用gzip -t验证完整性再随机抽查几条转录本 ID 是否在 FASTA 中存在。Ensembl 还提供 BioMart 接口做基因注释批量导出比直接下载整个 GFF3 更灵活。在 BioMart 页面选择数据集、过滤条件和属性可以只导出目标基因列表的染色体位置、基因名、转录本数量等信息适合不想要全基因组注释文件的场景。4.3 UCSC Table Browser 下载与字段设置UCSC 的 Table Browser 适用于从 UCSC 下载特定区域的注释数据。核心操作分四步第一步选择 clade如 Mammal、genome如 Human和 assembly如 GRCh38第二步选择 group 和 track比如 RefSeq Genes第三步在 region 中输入 chr1:100000-200000 之类的区间或者直接选 whole genome第四步在 output format 里选 BED 或 GFF并决定是否包含 header 行。# 用 UCSC 的 REST API 下载特定区域序列 curl https://api.genome.ucsc.edu/getData/sequence?genomehg38;chromchr1;start100000;end200000 -o region.fa提示UCSC API 的坐标是 0-based 半开区间start100000;end200000实际覆盖的是 chr1 第 100001 到第 200000 个碱基1-based 计数。和 GFF 的 1-based 闭合区间对照时GFF 里的chr1 100000 200000在 UCSC API 里应写成start99999;end200000。这个转换做错一次整段序列就偏移一个碱基如果是编码区直接导致移码。4.4 一个最小可用的格式互转脚本课件里没有给代码但实际工作中格式互转是最高频的需求。这里给一个 FASTA 转 GFF 的最小 Python 脚本用于给一段序列快速生成粗略的 gene 注释行。#!/usr/bin/env python3 将单条 FASTA 序列转为最小 GFF3 文件仅 gene 特征。 import sys def fasta_to_gff(fasta_file, gff_file, sourcecustom): seq_id None seq_len 0 with open(fasta_file, r) as f, open(gff_file, w) as g: for line in f: line line.strip() if not line: continue if line.startswith(): # 输出上一条序列的 gene 行 if seq_id and seq_len 0: # GFF3 坐标是 1-based 闭合区间 g.write(f{seq_id}\t{source}\tgene\t1\t{seq_len}\t.\t\t.\tID{seq_id}\n) seq_id line[1:].split()[0] seq_len 0 else: seq_len len(line) if seq_id and seq_len 0: g.write(f{seq_id}\t{source}\tgene\t1\t{seq_len}\t.\t\t.\tID{seq_id}\n) if __name__ __main__: fasta_to_gff(sys.argv[1], sys.argv[2])提示脚本的核心逻辑是逐行扫描 FASTA遇到结束前一条序列并输出 GFF 行最后在文件结束处补上最后一条序列的注释。seq_id取描述行第一个字段seq_len累加所有非行的碱基长度。输出时 GFF 的起始写 1、终止写总长度、链写。这个脚本只生成 gene 级注释不识别 CDS 和 exon真实项目中建议用gffread、AGAT或BioPython做完整转换。5. 避坑LOCUS 行、坐标体系与质量值编码最容易翻车的五个场景这部分内容基本是课件里那些「看着简单用起来全是问题」的知识点的集合也是我带新人时反复强调的雷区。5.1 现象从 NCBI 下载的 GenBank 文件用 BioPython 解析报错说LOCUS行长度不足原因BioPython 的SeqIO对 LOCUS 行的解析要求相对严格不同来源的 GenBank 文件在 LOCUS 行缩进和字段间距上有差异。NCBI 官方下载的文件通常规范但从 SRA 转换出来或第三方站点下载的文件LOCUS 行可能被截断或字段缺失。 解决先用文本编辑器打开看 LOCUS 行是否完整包含六个字段。缺失严重时不要手工修直接回到 NCBI 用 E-utilities 重新下载原始记录或者用efetch的rettypegbwithparts获取完整版本。5.2 现象GFF 文件在 IGV 里打开后基因位置整体往右偏移了一个碱基原因GFF 是 1-based 闭合区间而 IGV 内部使用 0-based 半开坐标。例如 GFF 里写100 200IGV 实际显示的是从第 99 个碱基0-based到第 200 个碱基的区域。这个偏移单个基因不明显但连锁基因一起偏就会导致注释和比对峰图错位。 解决不需要手动改 GFFIGV 会自动处理这个转换。真正要注意的是自己写脚本时不要混用两种坐标系——从一个格式转另一个格式时起点要减 1GFF 转 BED终点保持不变。我一般会在脚本里写死这个规则并加注释防止自己后来忘了。5.3 现象FASTQ 质量值用fastqc检查报出大量bad quality但实际数据质量尚可原因质量值编码判断错误。fastqc会自动猜测编码类型但如果文件里恰好只有极端的质量值字符可能导致误判为 Phred64。另一个常见原因是测序平台本身用了老式编码如 Solexa 的 Phred64但下机数据没有做编码转换。 解决优先用fastqc输出的Encoding字段确认编码方式。如果从!到5之间的字符占比很高基本是 Phred33如果字符集中在B到h之间考虑 Phred64。批量处理前先随机抽取 1000 条记录做编码检查不要对整批数据直接跑质控。5.4 现象下载 Ensembl GFF3 和 FASTA 之后工具报错说转录本 ID 对不上原因GFF3 里的转录本 ID 如ENST00000456328.2但 FASTA 下载的是.cdna.all.fa.gz版本号可能不一致。常见情况是 GFF3 下载了当前 releaseFASTA 下载的是旧 release 缓存导致部分转录本 ID 在新版本中已经合并或删除。 解决下载时强制检查 release 版本GFF3 和 FASTA 必须来自同一个 release 目录。Ensembl 的 FTP 上current是软链接但下载工具如wget有时不会解析软链接建议先curl看一下current指向的实际版本号再拼接完整路径下载。5.5 现象KEGG 通路富集结果里出现大量「Unknown」条目原因基因 ID 类型没有统一。KEGG 富集分析要求输入的基因 ID 是 Entrez Gene ID 或 KEGG ID如果输入的是 Ensembl Transcript ID映射过程中大部分条目匹配不上。这个问题的根源是没搞清楚不同数据库 ID 之间的对应关系而这个对应关系才是注释环节真正的「体力活」。 解决先用biomart或gprofiler做 ID 转换再跑富集分析。转换时注意一对多的情况——一个 Ensembl 基因可能对应多个 Entrez ID建议保留多个映射结果而不是随机取第一个。6. 把课件变成自己的工作台一张速查表与批量提取 CDS 的脚本课件最后讲到的内容可以进一步沉淀成工具。我自己的习惯是每次讲完这个课件就给学员发两张表一张是数据库入口速查表一张是格式选择速查表。数据库速查表只写三列——库名、访问入口、典型用途格式速查表写四列——格式、后缀、核心结构、常见坑。这两张表的价值在于把课件里的知识点变成了「不需要二次思考就能查」的东西尤其对刚入门的人来说比反复翻 PPT 效率高得多。配合速查表还有一个高频场景值得做一个小工具从 GenBank 批量提取 CDS 序列。课件里提到的/translation限定符是提取蛋白序列的捷径但核酸序列需要按 FEATURES 坐标从 ORIGIN 部分截取。下面这个 Python 脚本用 BioPython 实现批量提取。#!/usr/bin/env python3 从 GenBank 文件批量提取 CDS 核酸序列与蛋白序列。 from Bio import SeqIO import sys def extract_cds(genbank_file, output_nuc, output_prot): nuc_handle open(output_nuc, w) prot_handle open(output_prot, w) for record in SeqIO.parse(genbank_file, genbank): for feature in record.features: if feature.type ! CDS: continue # 取 locus_tag 或 protein_id 作为序列名 seq_id feature.qualifiers.get(protein_id, [unknown])[0] # 从原始序列按坐标切片Python 切片是 0-based 半开区间 if join in str(feature.location): nuc_seq feature.location.extract(record.seq) else: start feature.location.start end feature.location.end nuc_seq record.seq[start:end] # 蛋白序列优先取 translation 限定符缺失时自行翻译 prot_seq feature.qualifiers.get(translation, [])[0] if not prot_seq: prot_seq str(nuc_seq.translate()) nuc_handle.write(f{seq_id}\n{nuc_seq}\n) prot_handle.write(f{seq_id}\n{prot_seq}\n) nuc_handle.close() prot_handle.close() if __name__ __main__: extract_cds(sys.argv[1], sys.argv[2], sys.argv[3])提示脚本用 BioPython 解析 GenBank 特征表只取type CDS的特征。核酸序列提取用的是feature.location.extract()方法这个方法已经处理了正负链和join拼接的情况——比手写坐标切片安全得多尤其是遇到跨外显子的 CDS。蛋白序列优先用translation限定符这是 GenBank 记录里注释好的成品序列如果该字段缺失再调用nuc_seq.translate()兜底。输出文件名通过命令行参数传入。这套东西用顺手之后其实等于把课件的知识全部变成了长期管用的工具。我自己每次拿到新的参考基因组注释都会先强制走一遍流程确认数据库版本、检查格式坐标系、抽查 CDS 提取结果、建立本地索引。从那以后我几乎没再遇到坐标偏移或者 ID 对不上的情况写论文跑项目都省心不少。希望帮到你——认真过一遍这个课件值得。本文还有配套的精品资源点击获取
返回列表