生物信息资源与工具 ==================== 生物信息资源与工具是现代生物信息学研究的重要基础。这些资源包括各种数据库、检索工具、数据标准和交换格式，它们共同构成了生物信息学研究和数据分析的基础设施。理解和掌握这些资源的使用方法，对于开展生物信息学研究至关重要。本章将系统介绍国际主要数据中心、核心数据库分类、数据检索与提交方法、生物本体与术语库，以及数据标准与交换格式，帮助读者建立起完整的生物信息资源知识体系。 关键术语解释 ------------ 在深入学习本章内容之前，读者需要了解一些关键术语的含义。 **数据库（Database）**：数据库是按照特定结构组织和存储数据的集合。在生物信息学中，数据库用于存储基因序列、蛋白质结构、表达谱等各种类型的生物数据。例如，GenBank是一个存储核酸序列的数据库，UniProt是一个存储蛋白质序列和功能的数据库。 **序列（Sequence）**：序列是指生物大分子中单体单元的排列顺序。核酸序列由四种碱基（A、T、C、G）组成，蛋白质序列由二十种氨基酸组成。序列是生物信息学研究的基本对象，包含了生物体的遗传信息。 **注释（Annotation）**：注释是对序列数据添加功能、结构、位置等描述性信息的过程。例如，标注一个DNA序列中哪些区域是基因、编码什么蛋白质、有什么功能等。高质量的注释可以帮助研究者理解序列的生物学意义。 **比对（Alignment）**：比对是将两个或多个序列进行比较，找出它们之间的相似性和差异性的过程。序列比对是生物信息学中最基础的分析方法之一，用于判断序列之间的同源性、预测基因功能等。 **本体（Ontology）**：本体是一种用于描述概念及其相互关系的规范化术语体系。在生物信息学中，本体用于标准化地描述基因功能、表型、疾病等概念，使得不同数据库之间的数据可以相互理解和关联。 国际主要数据中心 ---------------- 国际主要数据中心分布在全球各地，它们共同构成了生物信息数据的全球共享网络。这些数据中心不仅存储和管理着海量的生物信息数据，还提供了丰富的数据检索和分析工具，是生物信息学研究不可或缺的重要基础设施。当你需要查找某个基因的序列信息、下载公开的测序数据、检索相关文献时，都需要访问这些数据中心。熟悉这些数据中心的资源和使用方法，是开展生物信息学研究的基本功。 NCBI（美国国家生物技术信息中心）是全球最大的生物信息数据中心之一，位于美国马里兰州贝塞斯达。NCBI隶属于美国国立卫生研究院（NIH），于1988年成立，负责开发和维护多种重要的生物信息数据库和计算工具。NCBI管理的数据库包括PubMed文献数据库、GenBank核酸序列数据库、Protein蛋白质数据库、SNP单核苷酸多态性数据库、GEO基因表达数据库等，这些数据库是生物信息学研究中使用最频繁的数据资源。NCBI还提供了BLAST序列比对工具、Entrez检索系统、NCBI Genome Workbench等重要的分析工具。NCBI的数据资源对全球科研人员免费开放，是生物信息学领域最重要的数据来源之一。 .. code-block:: python # 使用Biopython访问NCBI数据库 from Bio import Entrez # 设置邮箱（NCBI要求） Entrez.email = "your_email@example.com" # 搜索PubMed文献 handle = Entrez.esearch(db="pubmed", term="TP53 AND cancer", retmax=10) record = Entrez.read(handle) handle.close() # 获取文献ID列表 pmid_list = record["IdList"] print(f"找到 {record['Count']} 篇相关文献") # 获取文献详细信息 handle = Entrez.efetch(db="pubmed", id=pmid_list, rettype="abstract") abstracts = handle.read() handle.close() EBI（欧洲生物信息学研究所）位于英国剑桥，是欧洲重要的生物信息数据中心。EBI成立于1994年，负责维护和开发欧洲的生物信息基础设施。EBI管理的数据库包括Ensembl基因组数据库、UniProt蛋白质数据库、EMBL-EBI核酸档案库（ENA）、ArrayExpress表达数据库、PDBe蛋白质结构数据库等。EBI还提供了BLAST@EBI、InterPro蛋白质功能分析、EMBOSS序列分析工具包等服务。EBI与NCBI、日本DDBJ并列为全球三大核酸数据中心，它们通过INSDC（国际核苷酸序列数据库协作）协议进行数据共享和同步。 .. code-block:: python # 使用REST API访问Ensembl数据库 import requests # 获取人类TP53基因信息 server = "https://rest.ensembl.org" endpoint = "/lookup/symbol/homo_sapiens/TP53" response = requests.get(server + endpoint, headers={"Content-Type": "application/json"}) if response.ok: gene_info = response.json() print(f"基因ID: {gene_info['id']}") print(f"基因描述: {gene_info['description']}") print(f"染色体位置: {gene_info['seq_region_name']}:{gene_info['start']}-{gene_info['end']}") print(f"基因方向: {gene_info['strand']}") DDBJ（日本DNA数据库）成立于1986年，位于日本静冈县。DDBJ是日本唯一的核苷酸序列数据库，也是INSDC的创始成员之一。DDBJ主要收集和发布来自日本研究机构的核酸序列数据，同时也接收全球范围的提交。DDBJ提供了与NCBI GenBank和EBI ENA的数据同步功能，研究者向任何一个数据库提交的数据都会被同步到其他两个数据库。DDBJ还开发了ARNIPS、SYSTINGS等工具用于序列数据的分析和管理。 CNCB/NGDC（中国国家生物信息中心）是中国科学院下属的生物信息学研究机构，负责建设和管理中国的生物信息基础设施。NGDC建立了多个重要的数据库，包括中国科学院蛋白质科学研究所（Institute of Bioinformatics, Chinese Academy of Sciences）等。NGDC还推出了中国国家基因组数据中心云平台，为中国科研人员提供数据存储和分析服务。中国的生物信息资源建设近年来发展迅速，正在成为全球生物信息数据共享的重要力量。 GISAID（全球流感数据共享倡议）成立于2008年，是专门用于共享流感病毒基因组数据的国际数据库。GISAID在新冠疫情期间发挥了关键作用，成为全球新冠病毒基因组数据的主要共享平台。GISAID采用独特的开放获取模式，在保护数据提供者权益的同时，促进了全球病毒数据的快速共享。 BOLD（生命条形码数据系统）是用于存储和管理DNA条形码数据的数据库系统。DNA条形码技术利用特定基因片段来鉴定物种，BOLD收集了来自全球的DNA条形码数据，为物种鉴定和生物多样性研究提供了重要资源。 PDB（蛋白质数据银行）是全球唯一的蛋白质三维结构数据库。PDB成立于1971年，最初在美国布鲁克海文国家实验室建立，现在由wwPDB（全球蛋白质数据银行）联盟共同维护。wwPDB联盟包括美国RCSB PDB、欧洲PDBe、日本PDBj和韩国PDB等成员。PDB存储了通过X射线晶体学、核磁共振和冷冻电镜等方法解析的蛋白质结构数据，是结构生物学研究的核心资源。 .. code-block:: python # 使用Biopython解析PDB文件 from Bio.PDB import PDBParser # 创建PDB解析器 parser = PDBParser() # 解析PDB文件（假设已有文件） structure = parser.get_structure("protein", "protein.pdb") # 遍历结构层级 for model in structure: for chain in model: print(f"链ID: {chain.id}") for residue in chain: # 获取残基信息 resname = residue.get_resname() res_id = residue.get_id() print(f" 残基: {resname} {res_id[1]}") 国家基因组科学数据中心是中国科学技术部下属的国家级基因组数据资源平台，隶属于中国科学院北京基因组研究所（国家生物信息中心）。该中心整合了中国的基因组数据资源，提供数据存储、共享和分析服务。 国家生物信息中心云平台是中国科学院北京基因组研究所推出的云计算服务平台，为生物信息学研究提供在线计算和分析工具，降低了生物信息分析的门槛。 核心数据库分类 -------------- 核心数据库按照数据类型可以进行更细致的分类，这种分类方式有助于研究者根据具体的研究需求快速定位所需的数据库资源。不同类型的数据库侧重于不同层面的生物信息数据，了解这些数据库的特点和用途是开展生物信息学研究的基础技能。掌握数据库分类可以帮助你快速找到所需的数据类型，例如，当你需要查找某个基因的DNA序列时，应该去核酸序列数据库；当你需要了解蛋白质的功能注释时，应该去蛋白质数据库；当你需要分析基因表达数据时，应该去表达数据库。正确选择数据库可以大大提高研究效率。 核酸序列数据库 ~~~~~~~~~~~~~~ 核酸序列数据库存储和管理着全球最重要的生物遗传信息。这些数据库不仅保存了核酸序列数据，还包含丰富的注释信息，如基因名称、功能描述、编码蛋白、变异信息等。当你克隆了一个新基因、设计PCR引物、分析基因变异、或者进行进化分析时，都需要从这些数据库获取参考序列。RefSeq数据库因其高质量的注释，常被用作分析的标准参考序列。 GenBank是NCBI建立的最重要的核酸序列数据库之一，是INSDC的核心组成部分。GenBank收录了来自全球所有公开可获得的核酸序列数据，每个序列条目都包含详细的注释信息。GenBank的数据每日与欧洲ENA和日本DDBJ进行同步更新，确保全球数据的完整性。GenBank采用Locus Definition格式存储数据，每条记录包含序列本身、注释信息和参考文献等部分。 .. code-block:: python # 解析GenBank格式文件 from Bio import SeqIO # 读取GenBank文件 for record in SeqIO.parse("sequence.gb", "genbank"): print(f"序列ID: {record.id}") print(f"序列描述: {record.description}") print(f"序列长度: {len(record.seq)} bp") # 获取注释信息 annotations = record.annotations if "organism" in annotations: print(f"物种: {annotations['organism']}") # 获取特征信息 for feature in record.features: if feature.type == "gene": gene_name = feature.qualifiers.get("gene", ["未知"])[0] print(f"基因: {gene_name}, 位置: {feature.location}") ENA（欧洲核苷酸档案）是EBI维护的核酸序列数据库，是欧洲最主要的核酸数据存档机构。ENA采用归档优先的策略，特别强调数据的长期保存和可追溯性。ENA提供了Wave工具用于序列数据的可视化和分析。 DDBJ作为日本唯一的核酸序列数据库，与GenBank和ENA保持数据同步。DDBJ特别重视亚洲地区的数据收集，是亚洲生物信息研究的重要数据来源。 RefSeq（参考序列数据库）是NCBI提供的经过人工审校的高质量参考序列数据库。与GenBank的自动注释不同，RefSeq的每条记录都经过专家审核，具有更高的准确性和可靠性。RefSeq为每种生物的基因组提供了完整的参考序列，包括染色体、线粒体、质粒等不同类型的序列，以及对应的基因、转录本和蛋白质注释。RefSeq使用统一的ID命名规范，如NM_001126表示mRNA，NR_001568表示非编码RNA，NP_001119表示蛋白质。 NCBI Nucleotide是NCBI核酸数据库的总入口，整合了GenBank、RefSeq、dbEST、dbSTS等多个子数据库的核酸序列信息。用户可以通过统一的界面搜索和获取不同来源的核酸序列数据。 SRA（序列读档档案）是NCBI存储高通量测序原始读数数据的数据库。SRA保存了来自各种测序项目的原始测序数据，包括Illumina、Ion Torrent、PacBio、Nanopore等不同平台的数据。SRA的数据量增长迅速，已成为全球最大的测序原始数据存储库。研究者可以从SRA下载公开的测序数据进行重新分析。 .. code-block:: bash # 使用SRA Toolkit下载SRA数据 # 安装: conda install -c bioconda sra-tools # 查询SRA数据库 esearch -db sra -query "SRP123456" | efetch -format runinfo # 下载SRA文件 prefetch SRR123456 # 转换为FASTQ格式 fastq-dump --split-files --gzip SRR123456 # 使用fasterq-dump（更快） fasterq-dump --split-files --threads 8 SRR123456 TSA（转录组鸟枪法组装）是NCBI存储转录组组装结果的数据库。TSA收录了来自不同物种、不同组织的转录组测序数据组装结果，为研究基因表达和可变剪接提供了重要资源。 EST（表达序列标签）是NCBI较早建立的数据库，存储了来自cDNA文库的单次测序序列。EST数据虽然相对老旧，但对于研究某些物种的基因表达仍然具有重要价值。 蛋白质序列数据库 ~~~~~~~~~~~~~~~~~~ 蛋白质序列数据库存储蛋白质氨基酸序列和功能注释信息。这些数据库不仅提供序列数据，还包含蛋白质功能、结构域、翻译后修饰、功能分类等丰富的注释信息。当你需要了解某个蛋白质的功能、结构域组成、翻译后修饰位点、或者进行蛋白质功能预测时，都需要查询蛋白质数据库。UniProt因其全面的功能注释，是蛋白质功能研究的首选数据库。 UniProt是全球最重要的蛋白质序列和功能数据库，由Swiss-Prot、TrEMBL和PIR三个部分整合而成。Swiss-Prot是经过人工审校的高质量蛋白质注释数据库，TrEMBL是自动注释的蛋白质数据，PIR是蛋白质信息资源数据库。UniProt为每条蛋白质记录提供详细的注释，包括蛋白质名称、基因名称、功能描述、亚细胞定位、翻译后修饰、变异信息等。UniProt还提供了UniProtKB/Swiss-Prot、UniProtKB/TrEMBL、UniRef（聚类）和UniParc（归档）等不同类型的数据服务。 .. code-block:: python # 使用UniProt REST API获取蛋白质信息 import requests # 获取TP53蛋白质信息 uniprot_id = "P04637" url = f"https://rest.uniprot.org/uniprotkb/{uniprot_id}.json" response = requests.get(url) if response.ok: data = response.json() print(f"蛋白质名称: {data['proteinDescription']['recommendedName']['fullName']['value']}") print(f"基因名称: {data['genes'][0]['geneName']['value']}") print(f"物种: {data['organism']['scientificName']}") # 获取功能注释 for comment in data['comments']: if comment['commentType'] == 'FUNCTION': print(f"功能: {comment['texts'][0]['value']}") NCBI Protein是NCBI的蛋白质序列数据库，整合了RefSeq Protein、GenPept、PDB等来源的蛋白质数据。NCBI Protein与GenBank核酸数据库保持密切关联，提供了核酸到蛋白质的映射关系。 RefSeq Protein是RefSeq项目的蛋白质部分，提供了经过人工审校的高质量蛋白质序列和注释。RefSeq蛋白质的ID以 ``NP_``（编码蛋白）、``XP_``（预测蛋白）等开头，与对应的mRNA（``NM_``、``XM_``）保持关联。 Pfam是蛋白质家族域数据库，使用隐马尔可夫模型（HMM）描述蛋白质家族和结构域。Pfam将蛋白质分为不同的家族，每个家族由多个序列比对构建的模型代表。Pfam广泛用于蛋白质功能注释和进化分析，是蛋白质家族分类的重要资源。 .. code-block:: python # 使用Pfam扫描蛋白质结构域 # 需要安装hmmer: conda install -c bioconda hmmer import subprocess # 准备蛋白质FASTA文件 protein_fasta = """ >TP53_human MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGP DEAPRMPEAAPPVAPAPAAPTPAAPAPAPSWPLSSSVPSQKTYQGSYGFRLGFLHSGTAK SVTCTYSPALNKMFCQLAKTCPVQLWVDSTPPPGTRVRAMAIYKQSQHMTEVVRRCPHHE RCSDSDGLAPPQHLIRVEGNLRVEYLDDRNTFRHSVVVPYEPPEVGSDCTTIHYNYMCNS SCMGGMNRRPILTIITLEDSSGNLLGRNSFEVRVCACPGRDRRTEEENLRKKGEPHHELP PGSTKRALPNNTSSSPQPKKKPLDGEYFTLQIRGRERFEMFRELNEALELKDAQAGKEPG GSRAHSSHLKSKKGQSTSRHKKLMFKTEGPDSD """ with open("protein.fasta", "w") as f: f.write(protein_fasta) # 运行hmmscan搜索Pfam结构域 # subprocess.run(["hmmscan", "--domtblout", "pfam_results.txt", "Pfam-A.hmm", "protein.fasta"]) InterPro是蛋白质特征整合数据库，综合了多个蛋白质家族和结构域数据库的信息，包括Pfam、SMART、PROSITE、PRINTS、Coils等。InterPro提供统一的蛋白质注释和分类服务，是蛋白质功能预测的常用工具。 SMART（简单模块架构研究工具）是蛋白质结构域数据库，专注于信号转导蛋白和蛋白质结构域的分析。SMART提供了结构域的注释和可视化功能。 PROSITE是蛋白质模式与谱数据库，收集了蛋白质序列中的功能位点、模式和各种生物学特征。PROSITE使用正则表达式描述蛋白质特征，用于蛋白质功能注释。 结构数据库 ~~~~~~~~~~ 结构数据库存储蛋白质和核酸的三维结构信息，这些结构数据对于理解生物大分子的功能和相互作用具有重要价值。 PDB（蛋白质数据银行）是全球唯一的蛋白质三维结构数据库，存储了通过X射线晶体学、核磁共振和冷冻电镜等方法解析的生物大分子结构。PDB采用PDB格式存储结构数据，每个结构条目包含原子坐标、晶体学参数、结构因子等信息。PDB还提供结构可视化、序列比对、功能注释等服务。 .. code-block:: python # 使用Biopython分析蛋白质结构 from Bio.PDB import PDBParser, DSSP import numpy as np # 解析PDB文件 parser = PDBParser() structure = parser.get_structure("1TUP", "1tup.pdb") # 计算两个原子之间的距离 def calc_distance(atom1, atom2): vector = atom1.coord - atom2.coord return np.sqrt(np.sum(vector * vector)) # 获取特定残基的原子 model = structure[0] chain = model['A'] residue = chain[(' ', 100, ' ')] # 获取CA原子 ca_atom = residue['CA'] print(f"CA原子坐标: {ca_atom.coord}") # 计算两个残基CA原子之间的距离 residue2 = chain[(' ', 101, ' ')] ca_atom2 = residue2['CA'] distance = calc_distance(ca_atom, ca_atom2) print(f"残基间距离: {distance:.2f} Å") PDBsum是PDBe维护的结构摘要数据库，提供PDB结构的详细分析和可视化。PDBsum包括结构概述、二级结构、配体信息、互作网络等内容，是结构分析的重要辅助工具。 SCOP（蛋白质结构分类）和CATH是两种主要的蛋白质结构分类系统。SCOP按照蛋白质的折叠类型、超家族、家族等层次进行分类，CATH按照类别（Class）、架构（Architecture）、拓扑（Topology）、同源性（Homologous superfamily）的层次进行分类。这两个系统帮助研究者理解蛋白质结构的进化关系。 AlphaFold DB是DeepMind开发的AlphaFold蛋白质结构预测系统的数据库，提供超过2亿个蛋白质结构的预测结果。AlphaFold DB的发布极大丰富了蛋白质结构信息，为结构生物学研究提供了重要补充。 .. code-block:: python # 从AlphaFold DB下载预测结构 import requests # AlphaFold DB API uniprot_id = "P04637" # TP53 url = f"https://alphafold.ebi.ac.uk/api/prediction/{uniprot_id}" response = requests.get(url) if response.ok: data = response.json() # 获取PDB文件下载链接 pdb_url = data[0]['pdbUrl'] print(f"AlphaFold预测结构下载链接: {pdb_url}") # 下载PDB文件 pdb_response = requests.get(pdb_url) with open("TP53_alphafold.pdb", "w") as f: f.write(pdb_response.text) EMDB（电子显微镜数据库）是存储冷冻电镜三维重构结果的数据库。随着冷冻电镜技术的发展，EMDB的数据量快速增长，成为结构生物学的重要数据来源。 BMRB（生物磁共振数据库）存储通过核磁共振技术确定的生物分子结构数据。BMRB与PDB保持数据同步，为NMR结构研究提供专门的数据服务。 基因组数据库 ~~~~~~~~~~~~~~~~~~ 基因组数据库存储完整基因组的序列和注释信息，是功能基因组学研究的核心资源。当你需要了解某个基因在基因组上的位置、基因结构、转录本信息、或者进行基因组比较分析时，都需要访问基因组数据库。Ensembl和UCSC Genome Browser因其丰富的注释轨道和强大的可视化功能，是基因组研究中最常用的工具。 Ensembl是EBI和Wellcome Sanger Institute联合开发的基因组数据库，提供高质量的人类和小鼠基因组注释。Ensembl还提供Ensembl Genomes子数据库，覆盖植物、真菌、原生生物、细菌等非动物物种。Ensembl的注释包括基因结构、转录本、变异、表达等多种类型，并提供丰富的浏览器和检索工具。 .. code-block:: python # 使用Ensembl REST API获取基因组信息 import requests server = "https://rest.ensembl.org" # 获取基因序列 gene_id = "ENSG00000141510" # TP53 endpoint = f"/sequence/id/{gene_id}" response = requests.get(server + endpoint, headers={"Content-Type": "application/json"}) if response.ok: data = response.json() print(f"基因ID: {data['id']}") print(f"序列长度: {len(data['seq'])} bp") print(f"序列前100bp: {data['seq'][:100]}") # 获取基因的转录本信息 endpoint = f"/lookup/id/{gene_id}?expand=1" response = requests.get(server + endpoint, headers={"Content-Type": "application/json"}) if response.ok: gene_data = response.json() print(f"\n基因名称: {gene_data['display_name']}") print(f"转录本数量: {len(gene_data['Transcript'])}") for transcript in gene_data['Transcript'][:3]: print(f" 转录本: {transcript['id']}, 长度: {transcript['length']} bp") UCSC Genome Browser是加州大学圣克鲁斯分校开发的基因组浏览器，提供人类、小鼠、大鼠等多个物种的基因组数据。UCSC Genome Browser以其强大的可视化功能和丰富的数据轨道著称，是基因组数据可视化的标准工具。 NCBI Genome是NCBI的基因组数据库门户，整合了多种生物的基因组序列和注释信息。NCBI Genome提供基因组装配、基因注释、变异数据等多种类型的数据服务。 GENCODE是专门提供人类和小鼠基因组基因注释的项目，GENCODE注释被认为是人类基因组注释的黄金标准。GENCODE由EBI和Wellcome Sanger Institute联合维护，注释质量经过严格审核。 1000 Genomes Project是重要的人类遗传变异资源，提供了来自不同人群的数千个个体的全基因组测序数据。该项目建立了人类遗传变异的基准图谱，是群体遗传学研究的基石。 WormBase、FlyBase、SGD、ZFIN、MGI分别是线虫、果蝇、酵母、斑马鱼、小鼠的模型生物基因组数据库。这些数据库为相应的模式生物研究提供专门的基因组资源和工具。 TAIR是拟南芥的基因组数据库，提供拟南芥基因组的序列、注释、功能等信息，是植物生物学研究的重要资源。 Phytozome是植物基因组数据库，整合了多种植物的基因组数据，支持跨物种的基因组比较分析。 转录调控数据库 ~~~~~~~~~~~~~~ 转录调控数据库存储转录因子结合位点、表观遗传修饰、基因调控网络等调控信息，是理解基因表达调控机制的重要资源。 ENCODE项目是旨在解析人类基因组功能元件的大型项目，产生了大量的人类基因组功能注释数据，包括转录因子结合位点、增强子、染色质状态等信息。ENCODE数据已成为理解基因调控的重要参考。 .. code-block:: python # 分析ENCODE ChIP-seq数据 # 使用pandas处理BED格式的peak文件 import pandas as pd # 读取ENCODE peak文件（BED格式） peak_file = "encode_peak.bed" columns = ['chrom', 'start', 'end', 'name', 'score', 'strand', 'signalValue', 'pValue', 'qValue', 'peak'] peaks = pd.read_csv(peak_file, sep='\t', header=None, names=columns) print(f"Peak总数: {len(peaks)}") print(f"Peak平均长度: {(peaks['end'] - peaks['start']).mean():.0f} bp") # 统计各染色体的peak数量 chrom_counts = peaks['chrom'].value_counts() print("\n各染色体Peak数量:") print(chrom_counts.head(10)) Roadmap Epigenomics项目绘制了人类表观遗传修饰的图谱，包括组蛋白修饰、DNA甲基化、染色质可及性等多种表观遗传标记。 JASPAR是转录因子结合谱数据库，收集了转录因子的位置权重矩阵（PWM）和结合位点信息。JASPAR数据广泛用于转录因子结合位点预测和ChIP-seq数据分析。 .. code-block:: python # 使用JASPAR数据库预测转录因子结合位点 from Bio import motifs from Bio.Seq import Seq # 从JASPAR获取TF motif # 这里使用示例数据 motif_matrix = { 'A': [3, 0, 0, 0, 0, 0], 'C': [0, 5, 0, 0, 0, 0], 'G': [2, 0, 5, 0, 5, 0], 'T': [0, 0, 0, 5, 0, 5] } motif = motifs.Motif(counts=motif_matrix) motif.background = {'A': 0.25, 'C': 0.25, 'G': 0.25, 'T': 0.25} # 搜索序列中的结合位点 test_seq = Seq("ACGTACGTGTCAGTACGT") for pos, score in motif.search(test_seq, threshold=5.0): print(f"位置: {pos}, 得分: {score:.2f}, 序列: {test_seq[pos:pos+motif.length]}") TRANSFAC是商业化的转录因子数据库，包含转录因子、结合位点、基因调控网络等信息。TRANSFAC提供专业版和免费版，是转录因子分析的重要工具。 ChIP-Atlas和CistromeDB是ChIP-seq和ATAC-seq数据的整合数据库，提供转录因子和表观遗传标记的结合位点信息。这些数据库允许研究者查询特定转录因子在不同细胞类型中的结合模式。 GTRD（转录调控区域数据库）整合了转录因子ChIP-seq数据，提供统一的数据检索和分析接口。 RegulonDB是大肠杆菌的基因调控数据库，是细菌调控网络研究的经典资源。 AnimalTFDB和KnockTF是动物转录因子和敲除影响数据库，提供转录因子分类和敲除实验数据。 代谢通路数据库 ~~~~~~~~~~~~~~ 代谢通路数据库存储生物代谢通路和分子网络信息，帮助理解细胞内的生物化学过程。 KEGG（京都基因与基因组百科全书）是日本京都大学开发的代谢通路数据库，是生物通路分析最常用的资源之一。KEGG数据库包括代谢通路、基因功能、疾病、药物等多个子系统，通过KEGG ID可以将基因、蛋白、代谢物、通路联系起来。KEGG Mapper工具允许研究者绘制自己的代谢通路图。 .. code-block:: python # 使用KEGG REST API获取通路信息 import requests # 获取KEGG通路信息 pathway_id = "hsa04110" # 细胞周期通路 url = f"http://rest.kegg.jp/get/{pathway_id}" response = requests.get(url) if response.ok: print(response.text[:500]) # 获取通路中的基因列表 url = f"http://rest.kegg.jp/link/genes/{pathway_id}" response = requests.get(url) if response.ok: genes = [] for line in response.text.strip().split('\n'): parts = line.split('\t') if len(parts) == 2: genes.append(parts[1]) print(f"\n通路中的基因数量: {len(genes)}") print(f"前10个基因: {genes[:10]}") Reactome是人工审校的生物通路数据库，提供人类和其他物种的代谢和信号通路注释。Reactome的注释经过专家审核，质量较高，并提供可视化和分析工具。 BioCyc是包含MetaCyc、EcoCyc等物种特异性通路数据库的集合。MetaCyc是实验验证的代谢通路数据库，EcoCyc是大肠杆菌的专门数据库。BioCyc提供物种特异性通路分析和预测功能。 WikiPathways是社区驱动的开放获取 pathway 数据库，允许研究者创建和维护通路图谱。 PANTHER（蛋白质分析通过进化关系）提供基于进化关系的蛋白质分类和通路注释。 表达数据库 ~~~~~~~~~~ 表达数据库存储基因表达谱数据，是转录组学研究的重要资源。 GEO（基因表达综合数据库）是NCBI最大的基因表达数据库，存储了来自各种实验的基因表达数据。GEO采用平台、系列、样本的三层结构组织数据，支持多种实验类型包括芯片表达谱、RNA-seq、ChIP-seq等。 .. code-block:: python # 使用GEOparse分析GEO数据 # 安装: pip install GEOparse import GEOparse # 下载并解析GEO数据集 # gse = GEOparse.get_GEO(geo="GSE12345", destdir="./") # 示例：处理已下载的数据 # gse = GEOparse.get_GEO(filepath="GSE12345_family.soft.gz") # 获取样本信息 # for gsm_name, gsm in gse.gsms.items(): # print(f"样本: {gsm_name}") # print(f"标题: {gsm.metadata.get('title', ['N/A'])[0]}") # 简化的表达矩阵处理示例 import pandas as pd # 模拟表达矩阵 expression_data = { 'Gene': ['TP53', 'BRCA1', 'MYC', 'EGFR'], 'Sample1': [10.5, 8.3, 15.2, 7.8], 'Sample2': [12.1, 9.1, 14.8, 8.2], 'Sample3': [8.9, 7.5, 16.1, 9.5] } df = pd.DataFrame(expression_data) df.set_index('Gene', inplace=True) print("基因表达矩阵:") print(df) # 计算基因的平均表达 df['Mean'] = df.mean(axis=1) print("\n平均表达值:") print(df['Mean']) ArrayExpress是EBI维护的表达数据库，与GEO类似，存储基因表达实验数据。ArrayExpress提供与实验相关的元数据信息，便于数据检索和再分析。 Expression Atlas是EBI提供的基因表达查询数据库，允许研究者查看基因在不同条件下的表达模式。 GTEx（基因型-组织表达）项目收集了来自多个组织的人类正常表达数据，是理解基因组织特异性表达的重要资源。 CCLE（癌细胞系百科全书）提供癌细胞系的基因表达、突变、药物敏感性等数据，是癌症研究的重要资源。 TCGA（癌症基因组图谱）虽然主要是突变数据库，但也包含丰富的表达数据，是癌症多组学研究的基石。 HPA（人类蛋白质图谱）提供蛋白质在各种组织中的表达信息，包括mRNA和蛋白质水平的表达数据。 互作数据库 ~~~~~~~~~~ 互作数据库存储蛋白质与蛋白质、蛋白质与核酸之间的相互作用信息，是构建分子互作网络和理解信号通路的重要资源。 STRING是最大的蛋白质互作数据库之一，整合了实验验证的互作数据、数据库记录和预测的互作关系。STRING提供互作网络可视化和功能富集分析。 .. code-block:: python # 使用STRING API获取蛋白质互作网络 import requests import json # STRING API基础URL string_api_url = "https://string-db.org/api" # 获取TP53的互作蛋白 params = { "identifiers": "TP53", "species": 9606, # 人类 "limit": 20, "caller_identity": "my_app" } # 获取互作网络 output_format = "json" method = "network" url = f"{string_api_url}/{output_format}/{method}" response = requests.get(url, params=params) if response.ok: interactions = response.json() print(f"找到 {len(interactions)} 个互作伙伴") for interaction in interactions[:10]: partner = interaction['preferredName_B'] score = interaction['score'] print(f"互作蛋白: {partner}, 置信度: {score:.3f}") BioGRID是实验验证的蛋白质和遗传互作数据库，收录了来自多种物种的互作实验数据。 HPRD（人类蛋白质参考数据库）是人类蛋白质的专门数据库，提供人类蛋白质的注释和互作信息。 IntAct和MINT是分子互作数据库，提供实验验证的蛋白质互作数据。 CORUM是哺乳动物蛋白复合物数据库，存储已知的蛋白质复合物信息。 DrugBank是药物-靶点互作数据库，包含药物与蛋白质靶点的相互作用信息，是药物研究的重要资源。 物种特异性数据库 ~~~~~~~~~~~~~~~~~~~~ 物种特异性数据库针对特定物种提供专门的基因组和功能注释，是模式生物研究的重要资源。 dbSNP是NCBI的单核苷酸多态性数据库，收录了人类和其他物种的SNP信息。dbSNP的ID以rs开头，是遗传变异研究的核心资源。 .. code-block:: python # 使用NCBI API查询dbSNP from Bio import Entrez Entrez.email = "your_email@example.com" # 查询特定SNP snp_id = "rs1042522" # TP53基因的SNP handle = Entrez.efetch(db="snp", id=snp_id, rettype="xml") record = Entrez.read(handle) handle.close() # 解析SNP信息 # print(record) dbVar和dbGaP是NCBI的结构变异和基因型-表型关联数据库。 ClinVar是临床意义变异数据库，收集了与疾病相关的遗传变异信息，是临床遗传学诊断的重要参考。 .. code-block:: python # 查询ClinVar数据库 import requests # ClinVar API gene = "TP53" url = f"https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi" params = { "db": "clinvar", "term": f"{gene}[gene] AND pathogenic[clinical_significance]", "retmax": 10, "rettype": "json" } response = requests.get(url, params=params) if response.ok: print(response.text[:500]) GWAS Catalog收录了全基因组关联研究（GWAS）发现的表型相关变异位点。 OMIM是人类孟德尔遗传在线数据库，收录了遗传病和基因的信息。 ExAC和gnomAD是群体频率数据库，提供了大规模人群的遗传变异频率信息，是变异解读的重要参考。 .. code-block:: python # 使用gnomAD API查询变异频率 import requests # gnomAD API v3 variant = "7-140753336-A-G" # 示例变异 url = f"https://gnomad.broadinstitute.org/api/variant/{variant}" # 注意：实际使用需要更详细的API调用 # 这里展示基本结构 非编码RNA数据库 ~~~~~~~~~~~~~~~~ 非编码RNA数据库存储各种非编码RNA的序列和功能信息。 miRBase是microRNA数据库，收录了已报道的microRNA序列和注释。 .. code-block:: python # 分析microRNA靶基因 # 使用miRBase获取miRNA信息 import requests # miRBase API mirna_id = "hsa-miR-21" url = f"https://mirbase.org/api/mirna/{mirna_id}" response = requests.get(url) if response.ok: data = response.json() print(f"miRNA ID: {data.get('name', 'N/A')}") print(f"序列: {data.get('sequence', 'N/A')}") NONCODE和lncRNAdb是长链非编码RNA数据库。 Rfam是RNA家族数据库，包含非编码RNA的家族和特征信息。 RNAcentral是整合多种RNA数据的统一数据库。 疾病与药物基因组学数据库 ~~~~~~~~~~~~~~~~~~~~~~~~~~ 疾病与药物基因组学数据库将遗传变异、药物反应和疾病联系起来，是精准医学研究的重要资源。 CTD（比较毒理基因组学数据库）整合了化学物质、基因和疾病之间的关联信息。 PharmGKB是药物基因组学知识库，收录了药物代谢和响应相关的遗传变异。 TTD（Therapeutic Target Database）是治疗靶点数据库。 ChEMBL、PubChem BioAssay、BindingDB是药物化合物和生物活性数据库。 IUPHAR/BPS Guide to PHARMACOLOGY是药理学指南数据库。 数据检索与提交 -------------- 数据检索与提交是生物信息学研究的基本技能。研究者不仅需要从各种数据库获取数据，有时也需要向数据库提交自己的研究成果。掌握这些技能对于开展生物信息学研究至关重要。 Entrez系统检索 ~~~~~~~~~~~~~~ Entrez系统是NCBI提供的统一检索平台，可以跨多个数据库进行信息检索。Entrez支持全局搜索，用户可以在一个界面中搜索PubMed、GenBank、Protein、GEO等数据库。Entrez还提供高级搜索语法，用户可以通过字段限定和布尔运算精确查找所需信息。E-utilities是Entrez的编程接口，包括ESearch（搜索）、ESummary（获取摘要）、EFetch（获取详细信息）、ELink（获取链接）等工具，支持自动化数据检索和处理。 .. code-block:: python # 使用E-utilities进行高级检索 from Bio import Entrez Entrez.email = "your_email@example.com" # 高级搜索：搜索特定基因在特定物种中的序列 search_term = "TP53[Gene] AND Homo sapiens[Organism] AND RefSeq[Filter]" handle = Entrez.esearch(db="nucleotide", term=search_term, retmax=20) record = Entrez.read(handle) handle.close() id_list = record["IdList"] print(f"找到 {record['Count']} 条记录") # 批量获取序列信息 handle = Entrez.efetch(db="nucleotide", id=id_list, rettype="gb", retmode="xml") records = Entrez.read(handle) handle.close() for record in records[:5]: print(f"序列ID: {record['GBSeq_primary-accession']}") print(f"描述: {record['GBSeq_definition'][:80]}...") print(f"长度: {record['GBSeq_length']} bp\n") .. code-block:: python # 使用E-utilities获取关联数据 from Bio import Entrez Entrez.email = "your_email@example.com" # 从Gene数据库获取基因信息 gene_id = "7157" # TP53 # 获取基因记录 handle = Entrez.efetch(db="gene", id=gene_id, rettype="xml") gene_record = Entrez.read(handle) handle.close() # 获取关联的PubMed文献 handle = Entrez.elink(dbfrom="gene", db="pubmed", id=gene_id) links = Entrez.read(handle) handle.close() if links[0]['LinkSetDb']: pmid_list = [link['Id'] for link in links[0]['LinkSetDb'][0]['Link']] print(f"关联的PubMed文献数: {len(pmid_list)}") # 获取文献摘要 if pmid_list: handle = Entrez.efetch(db="pubmed", id=pmid_list[:5], rettype="abstract") abstracts = handle.read() handle.close() print(abstracts[:500]) BLAST序列搜索 ~~~~~~~~~~~~~ BLAST是生物信息学最常用的序列比对工具，用于在数据库中搜索与查询序列相似的序列。BLAST使用启发式算法，在速度和灵敏度之间取得平衡。基本BLAST包括BLASTn（核酸比核酸）、BLASTp（蛋白比蛋白）、BLASTx（核酸翻译后比蛋白）、tBLASTn（蛋白比核酸翻译）、tBLASTx（核酸翻译比核酸翻译）等多种模式。PSI-BLAST通过多轮迭代搜索可以发现远缘相关的序列。DELTA-BLAST使用蛋白质结构域信息增强比对灵敏度。BLAST的输出格式包括文本、XML、JSON等多种形式，便于后续分析。BLAST是序列分析的基础工具，几乎所有涉及序列分析的研究都会用到。当你获得一个新的序列（如克隆的基因、测序结果、预测的蛋白质），需要了解它的功能时，第一步通常就是进行BLAST搜索，找出数据库中与之相似的已知序列，从而推断其可能的功能。BLAST也常用于验证序列的正确性、发现同源基因、分析序列变异等场景。 .. code-block:: python # 使用Biopython进行在线BLAST from Bio.Blast import NCBIWWW, NCBIXML # 准备查询序列 query_seq = """ MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGP DEAPRMPEAAPPVAPAPAAPTPAAPAPAPSWPLSSSVPSQKTYQGSYGFRLGFLHSGTAK SVTCTYSPALNKMFCQLAKTCPVQLWVDSTPPPGTRVRAMAIYKQSQHMTEVVRRCPHHE RCSDSDGLAPPQHLIRVEGNLRVEYLDDRNTFRHSVVVPYEPPEVGSDCTTIHYNYMCNS SCMGGMNRRPILTIITLEDSSGNLLGRNSFEVRVCACPGRDRRTEEENLRKKGEPHHELP PGSTKRALPNNTSSSPQPKKKPLDGEYFTLQIRGRERFEMFRELNEALELKDAQAGKEPG GSRAHSSHLKSKKGQSTSRHKKLMFKTEGPDSD """ # 运行BLASTp搜索 result_handle = NCBIWWW.qblast("blastp", "nr", query_seq, hitlist_size=10) # 解析结果 blast_records = NCBIXML.parse(result_handle) for record in blast_records: print(f"查询序列: {record.query[:50]}...") print(f"数据库: {record.database}") print(f"比对数量: {len(record.alignments)}\n") for alignment in record.alignments[:5]: print(f"比对序列: {alignment.hit_id}") print(f"描述: {alignment.hit_def[:80]}...") print(f"E值: {alignment.hsps[0].expect:.2e}") print(f"一致性: {alignment.hsps[0].identities}/{alignment.hsps[0].align_length}") print() .. code-block:: bash # 本地BLAST使用示例 # 安装: conda install -c bioconda blast # 创建BLAST数据库 makeblastdb -in reference_proteins.fasta -dbtype prot -out my_database # 运行BLASTp blastp -query query.fasta -db my_database -out results.txt -evalue 0.001 -outfmt 6 # 输出格式说明（-outfmt 6为表格格式）: # qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore # 运行BLASTn blastn -query nucleotide_query.fasta -db nt_database -out blastn_results.txt -evalue 0.001 # 运行BLASTx（核酸翻译后比对蛋白质数据库） blastx -query mrna.fasta -db protein_db -out blastx_results.txt -evalue 0.001 # 运行tBLASTn（蛋白质比对核酸翻译数据库） tblastn -query protein.fasta -db nucleotide_db -out tblastn_results.txt -evalue 0.001 数据提交工具 ~~~~~~~~~~~~ 数据提交工具用于向各种数据库提交序列和实验数据。BankIt是向GenBank提交数据的Web工具。Sequin是NCBI提供的桌面提交工具。Webin是ENA的数据提交门户。SRA submission用于提交高通量测序原始数据。BioProject和BioSample是NCBI用于组织项目元数据和样本信息的系统。 .. code-block:: python # 准备GenBank提交数据 from Bio import SeqIO from Bio.Seq import Seq from Bio.SeqRecord import SeqRecord from Bio.SeqFeature import SeqFeature, FeatureLocation # 创建序列记录 seq = Seq("ATGGCCGTTACCGAAGGTGAAACCGTTGAAATGACCGGTGAACTGACCGGTGAACTG") record = SeqRecord( seq, id="MYGENE001", name="MyGene", description="Example gene for submission", annotations={ "molecule_type": "DNA", "organism": "Homo sapiens", "source": "human tissue sample", } ) # 添加基因特征 gene_feature = SeqFeature( FeatureLocation(0, 57), type="gene", qualifiers={ "gene": "example_gene", "locus_tag": "EG001" } ) record.features.append(gene_feature) # 添加CDS特征 cds_feature = SeqFeature( FeatureLocation(0, 57), type="CDS", qualifiers={ "gene": "example_gene", "codon_start": 1, "transl_table": 1, "product": "example protein", "translation": str(seq.translate()) } ) record.features.append(cds_feature) # 写入GenBank格式文件 with open("submission.gb", "w") as f: SeqIO.write(record, f, "genbank") print("GenBank提交文件已创建: submission.gb") 生物本体与术语库 ------------------ 生物本体与术语库为生物信息数据提供统一的语义框架，使得不同数据库之间可以进行语义互操作。本体库定义了概念之间的关系，如包含（is_a）、部分（part_of）等，使得计算机可以理解数据的含义。生物本体是功能富集分析的基础。当你进行RNA-seq或芯片数据分析后得到一组差异表达基因，想要了解这些基因在功能上有什么共同特点时，就需要使用GO进行功能富集分析。本体也用于标准化不同数据库之间的术语，使得数据可以相互关联和比较。 Gene Ontology（GO）是使用最广泛的生物本体，提供了描述基因功能的标准化词汇。GO本体包括分子功能、生物过程和细胞组分三个本体，每个概念都有唯一的GO ID。GO注释将基因与GO术语关联起来，用于基因功能富集分析。 .. code-block:: python # GO富集分析示例 from scipy import stats import pandas as pd # 模拟GO注释数据 go_annotations = { 'Gene': ['TP53', 'BRCA1', 'MYC', 'EGFR', 'KRAS', 'PIK3CA'], 'GO_terms': [ ['GO:0006915', 'GO:0006355', 'GO:0003700'], # TP53 ['GO:0006281', 'GO:0006915', 'GO:0005515'], # BRCA1 ['GO:0006355', 'GO:0003700', 'GO:0008283'], # MYC ['GO:0007165', 'GO:0005524', 'GO:0004672'], # EGFR ['GO:0007264', 'GO:0005525', 'GO:0003924'], # KRAS ['GO:0004674', 'GO:0005524', 'GO:0045944'] # PIK3CA ] } df = pd.DataFrame(go_annotations) # 统计GO term出现频率 all_go_terms = [] for terms in df['GO_terms']: all_go_terms.extend(terms) go_counts = pd.Series(all_go_terms).value_counts() print("GO term出现频率:") print(go_counts) # GO term描述（示例） go_descriptions = { 'GO:0006915': 'apoptotic process', 'GO:0006355': 'regulation of transcription, DNA-templated', 'GO:0003700': 'DNA-binding transcription factor activity', 'GO:0006281': 'DNA repair', 'GO:0005515': 'protein binding', 'GO:0008283': 'cell population proliferation', 'GO:0007165': 'signal transduction', 'GO:0005524': 'ATP binding', 'GO:0004672': 'protein kinase activity' } print("\nGO term描述:") for go_id in go_counts.index[:5]: desc = go_descriptions.get(go_id, 'Unknown') print(f"{go_id}: {desc} (出现 {go_counts[go_id]} 次)") 疾病本体（DO）提供疾病概念的标准化描述。表型本体（HPO）描述人类表型特征。细胞类型本体（CL）标准化细胞类型命名。解剖学本体（UBERON）提供解剖学术语的统一描述。化学实体生物学本体（ChEBI）描述化学分子和功能基团。序列本体（SO）描述核酸和蛋白质序列特征。分子相互作用本体（MI）描述分子互作类型。NCBI Taxonomy提供生物分类学命名。Uberon是跨物种解剖本体。Experimental Factor Ontology（EFO）描述实验因素。 数据标准与交换格式 ------------------ 数据标准与交换格式是生物信息学数据管理和分析的基础。统一的数据格式使得不同软件和数据库之间可以交换和整合数据。当你使用不同的分析软件处理数据、与协作者共享数据、或者将数据提交到公共数据库时，都需要遵循相应的数据标准格式。了解这些格式的结构，可以帮助你正确解析和处理数据文件，避免因格式问题导致的分析错误。 FASTA和FASTQ是存储序列信息的文本格式。FASTA格式用于存储序列，可以是核酸或蛋白质。FASTQ在FASTA基础上增加了质量值信息，用于存储测序读数。 .. code-block:: python # FASTA格式示例 fasta_format = """ >gene1 Homo sapiens TP53 ATGGAGGAGCCGCAGTCAGATCCTAGCGTCGAGCCCCCTCTGAGTCAGGAAACATTTTCAG ACCTATGGAAACTACTTCCTGAAAACAACGTTCTGTCCCCCTTGCCGTCCCAAGCAATGGA TGATTTGATGCTGTCCCCGGACGATATTGAACAATGGTTCACTGAAGACCCAGGTCCAGAT """ # FASTQ格式示例 fastq_format = """ @SEQ_ID GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT + !''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65 """ # 解析FASTA文件 from Bio import SeqIO # 写入示例FASTA文件 with open("example.fasta", "w") as f: f.write(fasta_format.strip()) # 读取FASTA文件 for record in SeqIO.parse("example.fasta", "fasta"): print(f"序列ID: {record.id}") print(f"描述: {record.description}") print(f"序列长度: {len(record.seq)}") print(f"序列前50bp: {record.seq[:50]}") .. code-block:: python # 解析FASTQ文件 from Bio import SeqIO # 写入示例FASTQ文件 fastq_content = """@SEQ_001 GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT + !''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65 @SEQ_002 ATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCG + IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII """ with open("example.fastq", "w") as f: f.write(fastq_content) # 读取FASTQ文件 for record in SeqIO.parse("example.fastq", "fastq"): print(f"序列ID: {record.id}") print(f"序列: {record.seq}") print(f"质量值: {record.letter_annotations['phred_quality'][:20]}") print() SAM/BAM/CRAM是存储序列比对结果的标准格式。SAM是文本格式，BAM是二进制压缩格式，CRAM是更高效的压缩格式。这些格式是高通量测序数据分析的事实标准。 .. code-block:: python # SAM格式示例 sam_format = """ @HD\tVN:1.6\tSO:coordinate @SQ\tSN:chr1\tLN:248956422 @RG\tID:sample1\tSM:sample1 read001\t99\tchr1\t1000\t60\t50M\t=\t1500\t550\tATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCG\tIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII\tRG:Z:sample1 read001\t147\tchr1\t1500\t60\t50M\t=\t1000\t-550\tATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCG\tIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII\tRG:Z:sample1 """ # 使用pysam处理BAM文件 # 安装: pip install pysam # import pysam # 打开BAM文件 # bamfile = pysam.AlignmentFile("aligned.bam", "rb") # 遍历比对结果 # for read in bamfile.fetch('chr1', 1000, 2000): # print(f"Read ID: {read.query_name}") # print(f"位置: {read.reference_start}") # print(f"比对质量: {read.mapping_quality}") # print(f"序列: {read.query_sequence}") # bamfile.close() # SAM格式字段说明 print("SAM格式字段说明:") print("1. QNAME: 序列名称") print("2. FLAG: 比对标志位") print("3. RNAME: 参考序列名称") print("4. POS: 比对起始位置") print("5. MAPQ: 比对质量") print("6. CIGAR: CIGAR字符串") print("7. RNEXT: 配对序列的参考序列名") print("8. PNEXT: 配对序列的位置") print("9. TLEN: 模板长度") print("10. SEQ: 序列") print("11. QUAL: 质量值") VCF/BCF/GVCF是存储遗传变异的格式。VCF是文本格式，BCF是二进制格式，GVCF是包含所有位点信息的VCF变体。 .. code-block:: python # VCF格式示例 vcf_format = """##fileformat=VCFv4.2 ##INFO= ##INFO= ##FORMAT= ##FORMAT= #CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tSample1\tSample2 chr1\t1000\trs123\tA\tG\t100\tPASS\tDP=100;AF=0.5\tGT:DP\t0/1:50\t0/0:50 chr1\t2000\trs456\tC\tT\t200\tPASS\tDP=150;AF=0.3\tGT:DP\t0/1:75\t1/1:75 chr2\t500\t.\tG\tA\t50\tq10\tDP=80;AF=0.1\tGT:DP\t0/0:40\t0/1:40 """ # 使用pysam或cyvcf2处理VCF文件 # 安装: pip install cyvcf2 # from cyvcf2 import VCF # 读取VCF文件 # vcf = VCF("variants.vcf") # 遍历变异 # for variant in vcf: # print(f"染色体: {variant.CHROM}") # print(f"位置: {variant.POS}") # print(f"参考碱基: {variant.REF}") # print(f"变异碱基: {variant.ALT}") # print(f"质量: {variant.QUAL}") # print(f"信息: {variant.INFO}") # vcf.close() print("VCF格式说明:") print("CHROM: 染色体名称") print("POS: 变异位置（1-based）") print("ID: 变异ID（如rs号）") print("REF: 参考碱基") print("ALT: 变异碱基") print("QUAL: 变异质量") print("FILTER: 过滤标志") print("INFO: 附加信息") print("FORMAT: 基因型格式") print("样本列: 各样本的基因型数据") GFF/GTF/GFF3是存储基因组注释的格式。GFF3是目前最常用的版本。 .. code-block:: python # GFF3格式示例 gff3_format = """##gff-version 3 ##sequence-region chr1 1 248956422 chr1\tNCBI\tgene\t1000\t5000\t.\t+\t.\tID=gene001;Name=TP53;biotype=protein_coding chr1\tNCBI\tmRNA\t1000\t5000\t.\t+\t.\tID=transcript001;Parent=gene001;Name=TP53-001 chr1\tNCBI\texon\t1000\t1500\t.\t+\t.\tID=exon001;Parent=transcript001 chr1\tNCBI\texon\t2000\t2500\t.\t+\t.\tID=exon002;Parent=transcript001 chr1\tNCBI\texon\t4000\t5000\t.\t+\t.\tID=exon003;Parent=transcript001 chr1\tNCBI\tCDS\t1200\t1500\t.\t+\t0\tID=cds001;Parent=transcript001 chr1\tNCBI\tCDS\t2000\t2500\t.\t+\t1\tID=cds001;Parent=transcript001 chr1\tNCBI\tCDS\t4000\t4800\t.\t+\t2\tID=cds001;Parent=transcript001 """ # 解析GFF3文件 # 使用gffutils库 # 安装: pip install gffutils # import gffutils # 创建数据库 # db = gffutils.create_db("annotations.gff3", dbfn="annotations.db", force=True, keep_order=True) # 查询基因 # for gene in db.features_of_type('gene'): # print(f"基因ID: {gene.id}") # print(f"基因名称: {gene.attributes.get('Name', ['N/A'])[0]}") # print(f"位置: {gene.chrom}:{gene.start}-{gene.end}") # print(f"链: {gene.strand}") print("GFF3格式字段说明:") print("seqid: 序列ID（染色体）") print("source: 注释来源") print("type: 特征类型（gene, mRNA, exon, CDS等）") print("start: 起始位置（1-based）") print("end: 结束位置") print("score: 得分") print("strand: 链方向（+, -, .）") print("phase: 编码相位（0, 1, 2）") print("attributes: 属性键值对") BED/bigBed是存储基因组区间的格式。bigBed是二进制索引格式，适合大数据。 .. code-block:: python # BED格式示例 bed_format = """chr1\t1000\t2000\tgene001\t0\t+ chr1\t3000\t4000\tgene002\t0\t- chr2\t5000\t6000\tgene003\t0\t+ chr2\t7000\t8000\tgene004\t0\t- """ # BED格式字段说明 print("BED格式字段说明:") print("1. chrom: 染色体名称") print("2. chromStart: 起始位置（0-based）") print("3. chromEnd: 结束位置（半开区间）") print("4. name: 特征名称（可选）") print("5. score: 得分（可选）") print("6. strand: 链方向（可选）") # 使用pybedtools处理BED文件 # 安装: pip install pybedtools # import pybedtools # 读取BED文件 # bed = pybedtools.BedTool("regions.bed") # 区间操作 # intersect = bed.intersect("other_regions.bed") # merge = bed.merge() # sort = bed.sort() WIG/bigWig是存储连续数值的格式，如测序覆盖度。 .. code-block:: python # WIG格式示例 wig_format = """track type=wiggle_0 name="coverage" description="Sequencing coverage" fixedStep chrom=chr1 start=1000 step=1 10 15 20 18 22 variableStep chrom=chr1 span=10 2000 50 2010 55 2020 60 """ print("WIG格式说明:") print("fixedStep: 固定步长格式") print("variableStep: 可变步长格式") print("chrom: 染色体名称") print("start: 起始位置") print("step: 步长") print("span: 跨度") PDB/mmcIF/PDBe是存储蛋白质三维结构的格式。 .. code-block:: python # PDB格式示例 pdb_format = """HEADER PROTEIN 01-JAN-00 1ABC TITLE EXAMPLE PROTEIN STRUCTURE COMPND MOL_ID: 1; COMPND 2 MOLECULE: EXAMPLE PROTEIN; ATOM 1 N MET A 1 10.000 20.000 30.000 1.00 50.00 N ATOM 2 CA MET A 1 11.000 21.000 31.000 1.00 50.00 C ATOM 3 C MET A 1 12.000 22.000 32.000 1.00 50.00 C ATOM 4 O MET A 1 13.000 23.000 33.000 1.00 50.00 O END """ print("PDB格式说明:") print("ATOM: 原子记录") print(" - 序号: 原子序号") print(" - 原子名: 原子类型（N, CA, C, O等）") print(" - 残基名: 氨基酸三字母代码") print(" - 链ID: 蛋白质链标识") print(" - 残基序号: 残基编号") print(" - X, Y, Z: 三维坐标") print(" - 占有率: 原子占有率") print(" - 温度因子: B因子") print(" - 元素: 元素符号") mzML/mzXML是质谱数据的标准格式。 HDF5/H5AD是存储单细胞测序数据的格式。 .. code-block:: python # 使用h5py处理HDF5文件 # 安装: pip install h5py import h5py import numpy as np # 创建HDF5文件示例 with h5py.File("single_cell_data.h5", "w") as f: # 创建组 grp = f.create_group("expression") # 存储表达矩阵 expression_matrix = np.random.rand(100, 2000) # 100个细胞，2000个基因 grp.create_dataset("matrix", data=expression_matrix) # 存储基因名称 gene_names = [f"Gene_{i}" for i in range(2000)] grp.create_dataset("gene_names", data=np.array(gene_names, dtype='S')) # 存储细胞条形码 cell_barcodes = [f"Cell_{i}" for i in range(100)] grp.create_dataset("cell_barcodes", data=np.array(cell_barcodes, dtype='S')) print("HDF5文件结构:") print("expression/") print(" ├── matrix (100 x 2000 表达矩阵)") print(" ├── gene_names (2000个基因名称)") print(" └── cell_barcodes (100个细胞条形码)") MAGE-TAB和ISA-Tab是表达实验的元数据格式。 RO-Crate和BioCompute Objects是工作流和数据可重复性的标准。 通过本章的学习，读者应该对生物信息资源与工具有了全面的了解。这些数据库和工具是开展生物信息学研究的基础设施，熟练掌握它们的使用方法对于进行高质量的研究工作至关重要。在实际研究中，研究者需要根据具体的研究问题选择合适的数据库和工具，并遵循数据标准进行数据管理和共享。