一、生物序列下载

NCBI:美国国立生物技术信息中心(National Center for Biotechnology Information)

1、生物序列一般是去NCBI下载,它有很多数据库,其中Gene是一般下载基因序列的库,Genome是下载基因组序列的库,搜索栏里键入oct4并回车。

2、能看到该基因在不同物种和实验中所测得的相同基因序列,选择智人的POU5F1基因,POU5F1Oct4基因的别名,本质一样。

3、点击FASTA链接。

4、可以看到测序技术得到的DNA序列。

5、点击Send to,选中FileFormatFASTA,点击Create File就可以下载

二、Python实现处理基因序列的基本算法

目录结构:

bio
 ├─main.py
 └─sequence.fasta

1、获取规范化数据:

打开序列文件,观察信息。Fasta格式的文件里每条序列的都有以>开头的序列描述,也就是每个>及下一个>之间的碱基就是一条序列。

def get_fasta(f_name):
    '''获取序列规范化数据'''
    fasta = {}
    with open('sequence.fasta','r') as f:
        for line in f:
            if line.startswith('>'):
                name = line[1:].rstrip()  # 去除描述字段里的\n和>
                fasta[name] = ''
                continue
            # 去除序列字段里的\n并是字符规范为大写字符
            fasta[name] += line.rstrip().upper()
    return fasta

规范化的数据才具有生物学意义。

2、核苷酸计数,碱基偏好性:

数值可以查看碱基偏好性。例如:一定类型的小型RNA会有特定的碱基偏好性,它第一个碱基偏好U,可以用于评价数据质量。如果miRNA第一碱基不偏好U,说明数据或分析过程出错。

def nt_count(seq):
    '''核苷酸计数'''
    ntCounts = []
    for nt in ['A', 'C', 'G', 'T']:
        ntCounts.appen(seq.count(nt))
    return ntCounts

3、GC含量:

$(A+T)/(G+C)$随着DNA的种类不同而异。GC含量越高,DNA密度也越高,同时温度和酸碱不易使之变形,利用这一特性可进行DNA分离或测定。同时物种的GC含量有着特异性,由此可判断测序后的数据是否合格

from __future__ import division

def gc_content(seq):
    '''GC含量'''
    total = len(seq)
    gcCount = seq.count('G') + seq.count('C')
    gcContent = format(float(gcCount/total*100), '.6f')
    return gcContent

4、DNA转录为RNA

def dna_trans_rna(dnaSeq):
    '''DNA翻译为RNA'''
    rnaSeq = dnaSeq.replace('T', 'U')
    return rnaSeq

5、RNA翻译为蛋白质

def rna_trans_protein(rnaSeq):
    codonTable = {
        'AUA':'I', 'AUC':'I', 'AUU':'I', 'AUG':'M',
        'ACA':'T', 'ACC':'T', 'ACG':'T', 'ACU':'T',
        'AAC':'N', 'AAU':'N', 'AAA':'K', 'AAG':'K',
        'AGC':'S', 'AGU':'S', 'AGA':'R', 'AGG':'R',
        'CUA':'L', 'CUC':'L', 'CUG':'L', 'CUU':'L',
        'CCA':'P', 'CCC':'P', 'CCG':'P', 'CCU':'P',
        'CAC':'H', 'CAU':'H', 'CAA':'Q', 'CAG':'Q',
        'CGA':'R', 'CGC':'R', 'CGG':'R', 'CGU':'R',
        'GUA':'V', 'GUC':'V', 'GUG':'V', 'GUU':'V',
        'GCA':'A', 'GCC':'A', 'GCG':'A', 'GCU':'A',
        'GAC':'D', 'GAU':'D', 'GAA':'E', 'GAG':'E',
        'GGA':'G', 'GGC':'G', 'GGG':'G', 'GGU':'G',
        'UCA':'S', 'UCC':'S', 'UCG':'S', 'UCU':'S',
        'UUC':'F', 'UUU':'F', 'UUA':'L', 'UUG':'L',
        'UAC':'Y', 'UAU':'Y', 'UAA':'', 'UAG':'',
        'UGC':'C', 'UGU':'C', 'UGA':'', 'UGG':'W',
    }
    proteinSeq = ""
    for codonStart in range(0, len(rnaSeq), 3):
        codon = rnaSeq[codonStart:codonStart + 3]
        if codon in codonTable:
            proteinSeq += codonTable[codon]
    return proteinSeq

6、获取反向互补序列:

def reverse_comple(type, seq):
    '''获取反向互补序列'''
    seq = seq[::-1]
    dnaTable = {"A": "T", "T": "A", "C": "G", "G": "C"}
    rnaTable = {"A": "T", "U": "A", "C": "G", "G": "C"}
    res = ""
    if type == "dna":
        for ele in seq:
            if ele in seq:
                if type == "dna":
                    res += dnaTable[ele]
                else:
                    res += rnaTable[ele]
    return res

7、验证

def main():
    oct4 = get_fasta('sequence.fasta')
    for name, sequence in oct4.items():
        print ("name: ", name)
        print ("sequence: ", sequence)
        print ("nt_count: ", nt_count(sequence))
        print ("cg_content: ", cg_content(sequence))
        rna = dna_trans_rna(sequence)
        print ("rna: ", rna)
        protein = rna_trans_protein(rna)
        print ("protein: ", protein)
        print ("reverse_comple: ", reverse_comple("dna", sequence))

if __name__ == '__main__':
    main()

三、Biopython分析生物序列

Biopython

Biopython是一个使用Python开发的计算分子生物学工具。

几乎使用编程的生信都开发出了适应自用语言的包,比如BioPerlBioGoBioRubyPyCogentBioJavaSeqAn。作为Python爱好者当然选择Piopython。安装见 Bioconda生物信息环境搭建

测试安装:

$ python
>>> from Bio.Seq import Seq
>>> seq = Seq('ATCG')
>>> seq
Seq('ATCG')

基本用法

1、读取常见fastagenbank格式的序列文件

导入包:

from Bio import SeqIO

读取:

# 读取包含单个序列的Fasta格式文件
faSeq = SeqIO.read('sequence.fasta', 'fasta')
print(faSeq)

# 读取包含多个序列的fasta格式文件
for fa in SeqIO.parse('multi.fasta', 'fasta'):
    print(fa.seq)

# 一个多序列文件中的所有序列
seqs = [fa.seq for fa in SeqIO.parse("multi.fasta",  "fasta")]
print(seqs)
# 如果不想要seq对象中的字母表,可以用str()来强制类型转换
seqs = [str(fa.seq) for fa in SeqIO.parse("multi.fasta",  "fasta")]
print(seqs)

# 读取包含单个序列的 gb 格式文件
gbSeq = SeqIO.read("sequence.gb", "genbank")
print(gbSeq)

2、浏览fasta序列文件内容

from Bio import SeqIO

# 读取包含单个序列 Fasta 格式文件
faSeq = SeqIO.read("sequence.fasta", "fasta")

# 获取详细的信息
# # 提取基因ID,name
# # Fasta 文件中序列名所在行的第一个词被作为id和name
print ("id: ", faSeq.id)
print ("name: ", faSeq.name)
# # 基因Description是fasta文件格式中的第一行
print ("description: ", faSeq.description)
# # 序列
print ("seq: ", faSeq.seq)
# # 序列来源库信息(NCBI的数据库信息会包括数据库交叉引用)
print ("dbxrefs: ", faSeq.dbxrefs)
# # 全部序列的注释信息
print ("annotations: ", faSeq.annotations)
# # 序列中每个字母的注释信息
print ("letter_annotations: ", faSeq.letter_annotations)
# # 部分序列的注释信息
print ("features: ", faSeq.features)

3、浏览genbank序列文件内容

from Bio import SeqIO

# 读取包含单个序列的 gb 格式文件
gbSeq = SeqIO.read("sequence.gb", "genbank")    
print(Sseq)

# 获取详细的信息
# # 提取基因ID,name
# # gb文件中序列名包含比fasta更加详细的序列信息,下面分别是id和name
print("id: ", gbSeq.id)
print ("name: ", gbSeq.name)
# # 基因 Description 是fasta文件格式中的第一行
print("description: ",  gbSeq.description)
# # 序列信息, 这里的序列信息是以Biopython中的seq对象存储
print("seq: ", gbSeq.seq)
# # 序列来源库信息(NCBI的数据库信息会包括数据库交叉引用)
print("dbxrefs: ", gbSeq.dbxrefs)
# # 全部序列的注释信息
print("annotations: ", gbSeq.annotations)
# # 序列中每个字母的注释信息
print("letter_annotations: ", gbSeq.letter_annotations)
# # 部分序列的注释信息,SeqFeature 对象的形式保存了features table中的所有entries(如genes和CDS等)
print("features: ", gbSeq.features)
# # 该基因的物种信息
print("organism: ", gbSeq.annotations["organism"])
# # 关于序列的注释信息,相关数据库的交叉引用号
print("comment: ", gbSeq.annotations["comment"])
# # 序列来源的物种名
print("source: ", gbSeq.annotations["source"])
# # 该基因的分类学信息
print("taxonomy: ", gbSeq.annotations["taxonomy"])
# # 该基因的整理后的注释信息
print("structured_comment: ", gbSeq.annotations["structured_comment"])
# # 该基因序列相关的关键词
print("keywords: ", gbSeq.annotations["keywords"])
# # 该基因的相关文献编号,或递交序列的注册信息
print("references: ", gbSeq.annotations["references"])
# # 该基因的入库时,给的基因编号,以及在染色体上的位点信息
print("accessions: ", gbSeq.annotations["accessions"])
# # 该基因的分子类型,一般为 DNA
print("molecule_type: ", gbSeq.annotations["molecule_type"])
# # 该基因的数据文件划分方式
print("data_file_division: ", gbSeq.annotations["data_file_division"])
# # 基因发布时间
print("date: ", gbSeq.annotations["date"])
# # 该基因的更新版本
print("sequence_version: ", gbSeq.annotations["sequence_version"])
# # 该基因的拓扑结构
print("topology: ", gbSeq.annotations["topology"])

GeneBank比fasta格式更加详细,但对序列处理来说内存占用和运行时间比这些信息更加重要。这就使fasta成为在序列分析中常用的格式。

4、新建序列文件

导入包:

from Bio.Seq import Seq

创建序列对象:

# 新建一个DNA序列对象
dnaSeq = Seq("GGATGGTTGTCTATTAACTTGTTCAAAAAAGTATCAGGAGTTGTCAAGGCAGAGAAGAGAGTGTTTGCA", IUPAC.unambiguous_dna)
# 新建一个RNA序列对象
rnaSeq = Seq("GGATGGTTGTCTATTAACTTGTTCAAAAAAGTATCAGGAGTTGTCAAGGCAGAGAAGAGAGTGTTTGCA", IUPAC.unambiguous_rna)
# 新建一个蛋白质序列对象
proteinSeq = Seq("GGATGGTTGTCTATTAACTTGTTCAAAAAAGTATCAGGAGTTGTCAAGGCAGAGAAGAGAGTGTTTGCA", IUPAC.protein)

序列对象一段字符串和其对应的编码表所定义。上述代码中,字符串内容一样,唯一不同的是第二个参数IUPAC值不一样。IUPAC (International Union of Pure and Applied Chemistry ) 是一个制定化学相关标准的组织,Biopython所使用的编码表就是由它制定的,细节可参考http://www.bioinformatics.org/sms2/iupac.html ,详细定义如下:

名称编码表:ambiguous_dna_lettersGATCRYWSMKHBVDNunambiguous_dna_lettersGATCambiguous_rna_lettersGAUCRYWSMKHBVDNunambiguous_rna_lettersGAUCproteinARNDCQEGHILKMFPSTWYV

5、修改序列文件

在生物学意义上,序列不能随便更改,即序列不可变。如果强行修改,则报错TypeError: 'Seq' object does not support item assignment

dna_seq[0] = "G"

如果执意修改也可以,但不建议这么做

dna_seq_mutable = dna_seq.tomutable()
dna_seq_mutable[0] = "G"

6、操作序列文件

导入包:

from Bio.Seq import Seq
# 新建一个DNA序列对象
dnaSeq = Seq("GGATGGTTGTCTATTAACTTGTTCAAAAAAGTATCAGGAGTTGTCAAGGCAGAGAAGAGAGTGTTTGCA", IUPAC.unambiguous_dna)
# 序列信息
print("Sequence: ", dnaSeq)
# 序列长度
print("Length : ", len(dnaSeq))
# 单个核苷酸计数
print("G Counts: ", dnaSeq.count("G"))
# 获取反向序列
print("reverse: ", dnaSeq[::-1])
# 获取反向互补序列
print("Reverse complement: ", dnaSeq.complement())
# 获取蛋白质的反向互补序列,这里显然是报错的,因为蛋白序列没有这一属性
print("Protein reverse complement: ", proteinSeq.complement())

7、用Biopython将DNA转录为RNA

# 转录
# # 如果序列为编码链,那么直接转换
print("rna: ", dna_seq.transcribe())
# # 如果序列为模板链,就需要先转为编码链
transcribe_seq = dna_seq.reverse_complement().transcribe()
print("rna: ", transcribe_seq)

8、用Biopython将RNA翻译为蛋白质

# 翻译
print("protein: ", transcribe_seq.translate())
# # 如果翻译的是线粒体密码子,那么在参数中需要输入,其他参考 https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi?mode=c or ftp://ftp.ncbi.nlm.nih.gov/entrez/misc/data/gc.prt
print("protein: ", transcribe_seq.translate(table="Vertebrate Mitochondrial"))
# # 在现实生物世界中,一般在遇到终止密码子之后的序列不用翻译
print("protein: ", transcribe_seq.translate(table="Vertebrate Mitochondrial", to_stop=True))
# # 如果DNA序列为编码序列,可以直接翻译,DNA序列不是3的倍数时,报错
print("protein: ", dna_seq.translate())
# # 在细菌世界中,在细菌遗传密码中 GTG 是个有效的起始密码子,注意第一个密码子(正常情况下 GTG编码缬氨酸, 但是如果作为起始密码子,则翻译成甲硫氨酸)
bacterial_dna = Seq("GTGAAAAAGATGCAATCTATCGTACTCGCACTTTCCCTGGTTCTGGTCGCTCCCATGGCATAA", generic_dna)
print("protein: ", bacterial_dna.translate(table="Bacterial", to_stop=True))
print("protein: ", bacterial_dna.translate(table="Bacterial", cds=True))

9、其他

# 寻找TATA框
# # TATA框约在多数真核生物基因转录起始点上游约-30bp(-25~-32bp)处,基本上由A-T碱基对组成,是决定基因转录始的选择,为RNA聚合酶的结合处之一
print("TA Counts: ", dna_seq.count("TA"))

# GC含量
print ("GC Contenten", 100 * float(dna_seq.count("G") + dna_seq.count("C")) / len(dna_seq))

# 得到promoter序列
# 在寻找基因的promoter时(一般promoter的位点不确定),但是可以通过将起始位点左右2kb基因视为promoter
# 这里训练切取,将切取设起始位点为前10bp
print("Promoter seq: ",dna_seq[:10])

四、访问NCBI Entrez数据库

Entrez

Entrez 在线资源检索器是一组服务器端程序,为国家生物技术信息中心(NCBI)的Entrez查询和数据库系统提供稳定的接口。使用固定的URL语法,将一组标准输入参数转换为各种NCBI软件组件搜索和检索所请求数据所需的值。目前包括38个数据库,涵盖各种生物医学数据,包括核苷酸和蛋白质序列,基因记录,三维分子结构和生物医学文献。该在线资源检索器可以使用任何计算机语言(Perl,Python,Java和C ++等)将URL发送到应用程序服务器并解析响应。

注意事项

  • 最小化请求数

    • 如果任务需要搜索和/或下载大量记录,则使用Entrez历史记录批量上载和/或检索这些记录而不是对每条记录使用单独的请求会更有效
    • 可以使用单个EPost请求上载数千个ID
    • 可以使用一个EFetch请求下载数百个记录
  • 访问限制

    • 为了不使服务器过载,NCBI建议用户每秒发布不超过三个URL请求
    • 将大型作业限制在工作日的周末或东部时间晚上9:00到凌晨5:00之间
  • 设置邮箱

    • 使用email参数,这样如果遇到什么问题,NCBI可以通过邮件联系到你
    • 邮件的参数从2010年6月1日是强制的参数,所以每次必须告诉 NCBI 是谁在访问
  • URL字符处理

    • 所有参数使用小写字符
    • 参数没有必需的顺序,通常会忽略空值或不适当的参数
    • 避免在URL中放置空格,尤其是在查询中。如果需要空格,请使用加号(+)代替空格
    • 其他特殊字符(例如引号(“)或用于引用历史记录服务器上的查询键的#符号)应由其URL编码表示(%22表示”;%23表示#)

基本操作

1、一般参数设置

# 设置 email 参数,为了方便 NCBI 的工作人员可以联系到你
# 邮件的参数从2010年6月1日是强制的参数,所以每次必须告诉 NCBI 是谁在访问
Entrez.email = "example@163.com"
# 如果你是通过其他脚本调用,可以设定 tool 的名字,默认为 `Biopython`
Entrez.tool = "exampleScript"
# 可选参数,使用代理,一般在无法正常访问时设置
os.environ["http_proxy"] = "http://proxyhost.example.com:8080"

2、查看概况

2.1 查看目前 NCBI 所有数据库

# 查看数据库概况
# 获取 Entrez 所有数据库的句柄
hd_info = Entrez.einfo()
# 获取所有数据库列表
read_info = Entrez.read(hd_info)
for db in read_info['DbList']:
    print(db)

2.2 查看单个数据库概况

# 获取 Entrez 的 gene 数据库句柄
hd_info_gene = Entrez.einfo(db="gene")
read_info_gene = Entrez.read(hd_info_gene)
# 数据库名
print("DbName: ", read_info_gene["DbInfo"]["DbName"])
# 在 NCBI 首页顶部下拉菜单栏中的命名
print("MenuName: ", read_info_gene["DbInfo"]["MenuName"])
# 数据库描述
print("Description: ", read_info_gene["DbInfo"]["Description"])
# 数据库收录总数
print("Count: ", read_info_gene["DbInfo"]["Count"])
# 最新更新时间
print("LastUpdate: ", read_info_gene["DbInfo"]["LastUpdate"])
# Gene 数据库中可用的搜索关键字列表
print("FieldList: ", read_info_gene["DbInfo"]["FieldList"])
# 我们把它遍历下
for field in read_info_gene["DbInfo"]["FieldList"]:
    print("%(Name)s\t %(FullName)s\t %(Description)s" % field)

3、查询

3.1 全局搜索 | EGQuery

EGQuery:https://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.ESummary
这里只关注搜索关键字在数据库中所有的个数,而不关注它的具体内容
在实际使用中我们可以通过在这里得到的数字来确定下载策略

# 全局搜索
hd_egquery = Entrez.egquery(term="oct4")
read_egquery = Entrez.read(hd_egquery)
print(read_egquery)
for ele in read_egquery["eGQueryResult"]:
    print(ele["DbName"], ele["Count"], ele["Status"])

3.2 查询单个数据库中的基因 | ESearch

关于 ESearch 的官方文档 https://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.ESearch

# 在数据库搜索基因
# 搜索 Xenopus laevis 物种中名为 oct4 的基因
handle = Entrez.esearch(db="gene", term="oct4[Gene] AND Xenopus laevis[ORGN]")
read_gene = Entrez.read(handle)
print(read_gene)

3.3 查询基因详细描述信息 | Esummary

ESummary :https://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.ESummary

# 获取摘要
# 通过 id 来获取 item 的详细信息
hd_esummary = Entrez.esummary(db="gene", id="397784")
read_esummary = Entrez.read(hd_esummary)
# 获取该基因的详细描述
for key, value in read_esummary['DocumentSummarySet']['DocumentSummary'][0].items():
    print(key, value)

3.4 查询交叉引用条目 | Elink

Elink:https://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.ELink

# 搜索交叉引用条目
# 接下来我们看看id为5460的基因相关的文献资料
read_elink = Entrez.read(Entrez.elink(dbfrom="gene", db="pubmed", id="5460"))
print ("LinkSetDb: ", read_elink[0]["LinkSetDb"])
# 查看所有相关的目标库
for lsd in read_elink[0]["LinkSetDb"]:
    print (lsd["DbTo"], lsd["LinkName"], len(lsd["Link"]))
# 查看相关的所有文献Id
for link in read_elink[0]["LinkSetDb"][0]["Link"]:
    print (link["Id"])

3.5 其他操作

3.5.1 上传id列表到服务器 | EPost

EPost :https://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.EPost
为什么要上传列表到服务器?
你要上传的 id 的列表会以 url 的形式上传到服务器,这里有一个问题,如果 id 很多,就会导致url很长。但是在 HTTP 的协议中,上传一般以 GET 形式,这种方式会限制 url 的长度,也就是说如果用户上传的 URL 太长就会只能局限在一定的长度内,而不能完整的上传到服务器。 为了解决这个问题,只能使用 POST 方式上传,它没有限制文本长度,随后以 HTTP 头文件的形式上传服务器,并以历史记录的形式存储在服务器

# 上传历史记录
# EPost:https://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.EPost
id_list = ['379522', '397784', '398336']
hd_epost = Entrez.epost("gene", id=",".join(id_list))
read_epost = Entrez.read(hd_epost)
print("web_env: ", read_epost["WebEnv"])
print("query_key: ", read_epost["QueryKey"])

3.5.2 拼写建议与纠正 | Espell

这个模块用来纠正输入的查询词条

# 拼写建议
hd_espell = Entrez.espell(term="steem cell")
read_espell = Entrez.read(hd_espell)
print("Query: ", read_espell["Query"])
print("CorrectedQuery: ", read_espell["CorrectedQuery"])

3.5.3 下载 | EFetch

EFetch:https://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.EFetch所有参数组合:https://www.ncbi.nlm.nih.gov/books/NBK25499/table/chapter4.T._valid_values_of__retmode_and/?report=objectonly

# 下载
hd_efetch_gb = Entrez.efetch(db='nucleotide', id="397784", rettype='gb', retmode='text')
hd_efetch_fa = Entrez.efetch(db='nucleotide', id="397784", rettype='fasta')
print (hd_efetch_gb.read())
print (hd_efetch_fa.read())
with open("res/397784.txt", "w") as file:
    hd_efetch_ml = Entrez.efetch(db='pubmed', id="397784", rettype='medline', retmode='text')
    file.write(hd_efetch_ml.read())

with open("res/397784.txt") as file:
    read_medline = Medline.read(file)
    print("PMID", read_medline["PMID"])
    print("TI", read_medline["TI"])

3.5.4 解析大文件| parse

一般在 NCBI 中的资源会有较大的内存占用,
这里的parse使用迭代器的方式,而不是像列表全部加载,因此了避免了大文件读取时占满内存

  • Linux 系统下准备工作

下载实例文件:ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/ASN_BINARY/Mammalia/Homo_sapiens.ags.gz
下载格式转换工具:ftp://ftp.ncbi.nlm.nih.gov/toolbox/ncbi_tools/converters/by_program/gene2xml/linux64.gene2xml.gz

  • 在终端依次运行下列命令
mkdir ncbi
cd ncbi
mkdir ags
mkdir tool
cd tool
wget ftp://ftp.ncbi.nlm.nih.gov/toolbox/ncbi_tools/converters/by_program/gene2xml/linux64.gene2xml.gz
gunzip linux64.gene2xml.gz
mv linux64.gene2xml gene2xml
cd ../ags
wget ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/ASN_BINARY/Mammalia/Homo_sapiens.ags.gz
gunzip Homo_sapiens.ags.gz
../tool/gene2xml -b T -i Homo_sapiens.ags -o Homo_sapiens.xml
  • 下载目录结构如下,这里的Homo_sapiens.xml 大约有15G(2018.09)
ncbi
 ├─ags
 |  ├─Homo_sapiens.ags
 |  └─Homo_sapiens.xml
 └─tool
    └─gene2xml
  • 使用 Biopython 解析
# 解析大文件
hd_parse = open("Homo_sapiens.xml")
res_parse = Entrez.parse(hd_parse)
for record in res_parse:
     status = record['Entrezgene_track-info']['Gene-track']['Gene-track_status']
     if status.attributes['value']=='discontinued':
         continue
     geneid = record['Entrezgene_track-info']['Gene-track']['Gene-track_geneid']
     genename = record['Entrezgene_gene']['Gene-ref']['Gene-ref_locus']
     print(geneid, genename)

五、批量下载基因与文献

Entrez极其强大,Biopython将它几乎所有操作都封装为方法,更方便的使用Entrez。使用Entrez实现批量的序列与文献下载。

批量获取氨基酸序列数据

1、利用 Nucleotide 数据库来查询所有 oct4 基因的序列数据,为了展示基础的流程,这里采用逐条下载的方式

from Bio import Entrez,SeqIO

# 参数设置
Entrez.email = "example@163.com"
Entrez.tool  = "exampleScript"

# 查询 oct4 基因的在 Nucleotide 中的总数
hd_egquery = Entrez.egquery(term="oct4")
read_egquery = Entrez.read(hd_egquery)
total = 0
for ele in read_egquery["eGQueryResult"]:
    if ele["MenuName"] == "Nucleotide":
        total = ele["Count"]

# 得到查询 id 列表
hd_esearch = Entrez.esearch(db="nucleotide", term="oct4", retmax=total)
read_esearch = Entrez.read(hd_esearch)
# 这里我们只取前两个序列
ids = read_esearch["IdList"][:2]

# 用得到的 id 列表去下载每一条 fasta 文件,并合并,以便后续分析使用(比如进化树构建)
hd_efetch_fa = Entrez.efetch(db='nucleotide', id=ids, rettype='fasta')
read_efetch_fa = hd_efetch_fa.read()
with open("oct4.fasta","w") as file:
    file.write(read_efetch_fa)
# 同理你可以得到 xml 格式的序列信息
hd_efetch_xml = Entrez.efetch(db="nucleotide", id=ids, retmode="xml")
read_efetch_xml = Entrez.read(hd_efetch_xml)
print(read_efetch_xml)
hd_efetch_gb = Entrez.efetch(db="nuccore", id=ids, rettype="gb", retmode="text")
# 这里读取的是文本文件,保存为本地数据
read_efetch_gb = hd_efetch_gb.read()
with open("res/oct4.gb","w") as file:
    file.write(read_efetch_gb)

# 如果需要提取其中一些信息,可以按照以下步骤, 这里需要注意需要重新得到 efetch 句柄
hd_efetch_gb = Entrez.efetch(db="nuccore", id=ids, rettype="gb", retmode="text")
parse_efetch_gb = SeqIO.parse(hd_efetch_gb, "gb")
# 这里可以保存为 xls 或者 csv 格式
for ele in parse_efetch_gb:
    print(ele.name, ele.annotations['molecule_type'], ele.seq)

1.1 用历史记录特性提高效率

还记得上一篇教程中提到的历史记录吗?
利用这个特性,不仅可以减轻 Entrez 服务器的负载,更可以同时获取多条数据,节省大量时间精力

from Bio import Entrez

# 参数设置
Entrez.email = "example@163.com"
Entrez.tool  = "exampleScript"

hd_search = Entrez.esearch(db="nucleotide", term="oct4", usehistory="y")
read_search = Entrez.read(hd_search)
webenv = read_search["WebEnv"]
query_key = read_search["QueryKey"]

# 使用历史记录特性来进行搜索。
# Entrez 将会提前进行缓冲,提高查询效率
step = 5
total = 10
with open("res_env_oct4.fasta", "w") as res_file:
    for start in range(0, total, step):
        end = min(total, start+step)
        print("Download record %i to %i" % (start+1, end))
        hd_fetch = Entrez.efetch(db="nucleotide", rettype="fasta", retmode="text", retstart=start, retmax=step, webenv=webenv, query_key=query_key)
        records = hd_fetch.read()
        res_file.write(records)

批量获取参考文献

1、利用PubMed数据库来查询所有关于小鼠的文献资料,这里采用逐条下载的方式

from Bio import Entrez
from Bio import Medline

# 参数设置
Entrez.email = "example@163.com"
Entrez.tool  = "exampleScript"

# 用 esearch 在 pubmed 库中搜索关键字为 "mouse" 的文章
# RetMax 这个参数为每次返回的最大个数,因此如果把Count的值赋给RetMax就会获取全部的mouse的文章,这里为实例设置为100
hd_esearch = Entrez.esearch(db="pubmed", term="mouse", RetMax="100")
read_esearch = Entrez.read(hd_esearch)
idlist = read_esearch["IdList"]
print("Total: ", read_esearch["Count"])
# 用 efetch下载
hd_efetch = Entrez.efetch(db="pubmed", id=idlist, rettype="medline", retmode="text", )
# 用 Medline 来解析
parse_medline = Medline.parse(hd_efetch)
with open("mouse_pubmed.xls", "w") as file:
    file.write("title\tauthors\tsource\tPubMed\n")
    for i, ele in enumerate(list(parse_medline)):
        line = ele['TI'] + "\t" + ",".join(ele['AU']) + "\t" + ele['SO'] + "\t" + ele['PMID'] + "\n"
        file.write(line)
        print(i, line)

2、提高上面脚本的效率,这里查询近一年的关于 Sus scrofa 的综述

from Bio import Entrez
# 参数设置
Entrez.email = "example@163.com"
Entrez.tool  = "exampleScript"

# 搜索
hd_esearch = Entrez.esearch(db="pubmed", term="Sus scrofa", reldate=365, ptyp="Review", usehistory="y")
read_esearch = Entrez.read(hd_esearch)
total = int(read_esearch["Count"])
webenv = read_esearch["WebEnv"]
query_key = read_esearch["QueryKey"]

# 这里演示设定total为 10
total = 10
step = 5
print("Result items: ", total)
with open("recent_review_sus_scrofa.txt", "w") as file:
    for start in range(0, total, step):
        print("Download record %i to %i" % (start + 1, int(start+step)))
        hd_efetch = Entrez.efetch(db="pubmed", retstart=start, retmax=step, webenv=webenv, query_key=query_key, rettype="medline", retmode="text")
        file.write(hd_efetch.read())

获取物种谱系

NCBI 提供了很多生物相关数据库,用法几乎差不多,可以根据自身研究或者感兴趣的方向自行选择。
下面的例子是利用NCBI中的分类库 Taxonomy 来查询我们人类在分类学中的位置。

# 查看物种谱系
from Bio import Entrez

# 参数设置
Entrez.email = "example@163.com"
Entrez.tool  = "exampleScript"

# 在 Taxonomy 库中搜索 Homo sapiens
hd_esearch = Entrez.esearch(db="Taxonomy", term="Homo sapiens")
read_esearch = Entrez.read(hd_esearch)
id = read_esearch["IdList"][0]
hd_eftech = Entrez.efetch(db="Taxonomy", id=id, retmode="xml")
read_eftech = Entrez.read(hd_eftech)
print(read_eftech[0]["Lineage"])

参考链接:
白墨 - 知乎
Biopython 教程与手册 — Biopython-cn 0.1 文档
National Center for Biotechnology Information

Last modification:March 16th, 2020 at 01:00 pm
如果觉得我的文章对你有用,请随意赞赏