Python 自動化提取基因 CDS

wwwhjw1688com1622-8719999發表於2021-03-18

環境準備及背景介紹

環境搭建:Pycharm + Anaconda

目錄結構:

圖片

圖片

2

Python 實現

BaimoTools.py

  1#!/usr/bin/env python  2# -*- coding: utf-8 -*-  3# @Author  : Baimoc  4# @Email   : baimoc@163.com  5# @Time    : 2021/3/17 14:28  6# @File    : BaimoTools  7import os  8import time  9 10from Bio import SeqIO, SeqFeature 11 12 13class BaimoTools(): 14    def __init__(self, gb_file, fasta_file): 15        self.complete_fasta = "" 16        self.fasta_file = fasta_file 17        self.gb_file = gb_file 18        self.feature = None 19        self.record = None 20 21    def format_val(self, object=None): 22        """ 23        格式化物件值為字串 24        :param object: 物件或物件鍵值 25        :return: 26        """ 27        key = "" 28        # 判斷引數是否為字元 29        if isinstance(object, str): 30            obj = self.feature.qualifiers 31            key = object 32        else: 33            obj = object 34        # 為字元,提取 feature.qualifiers 物件關鍵字 35        if key != "" and not obj.get(key): 36            return "" 37        elif key == "": 38            val = obj 39        else: 40            val = obj[key] 41        # 轉換為字串 42        if not len(val): 43            val = "" 44        elif len(val) == 1: 45            val = val[0] 46        else: 47            if isinstance(val, SeqFeature.CompoundLocation) or isinstance(val, SeqFeature.FeatureLocation): 48                val = str(val) 49            else: 50                val = " | ".join(val) 51        return val 52 53    def extract_cds(self, cds): 54        """ 55        獲取 CDS 的 Fasta 序列 56        :param cds: 獲取指定基因的 CDS 區域,如果為空,則獲取全部 57        """ 58        records = list(SeqIO.parse(self.gb_file, "genbank")) 59        for record in records: 60            print(f"{record.id}") 61            for feature in record.features: 62                # 提取 CDS 資訊 63                try: 64                    if feature.type == "CDS": 65                        self.feature = feature 66                        self.record = record 67                        cds_gene = self.format_val('gene') 68                        if cds == "": 69                            self.complete_fasta += self.format_fasta() 70                        elif isinstance(cds, str) and cds_gene == cds: 71                            self.complete_fasta += self.format_fasta() 72                        elif isinstance(cds, list) and cds_gene in cds: 73                            self.complete_fasta += self.format_fasta() 74                 except: 75                    print(f"{record.id}: No CDS features") 76        self.write_file() 77 78    def write_file(self): 79        """ 80        寫入檔案 81        """ 82        with open(self.fasta_file, "w") as f_obj: 83            f_obj.writelines(self.complete_fasta) 84 85    def format_fasta(self, num=0): 86        """ 87        整理 Fasta 格式 88        :param num: 每行字元數,超出則換行 89        :return: Fasta 文字 90        """ 91        cds_gene = self.format_val('gene') 92        cds_location = self.format_val(self.feature.location) 93        cds_product = self.format_val('product') 94        cds_protein_id = self.format_val('protein_id') 95        cds_translation = self.format_val('translation') 96        complete_ana = f">{self.record.id} | {cds_gene} | {cds_product} | {cds_protein_id} | {str(cds_location)}\n" 97        format_seq = "" 98        if num: 99            for i, char in enumerate(cds_translation):100                format_seq += char101                if (i + 1) % num == 0:102                    format_seq += "\n"103        else:104            format_seq = cds_translation105        return complete_ana + format_seq + "\n"

圖片

3

使用示例

1

資料介紹

示例資料為新冠病毒的基因組 genbank 檔案,檔案中包含:

兩個基因組:LC553263.1 和 LC553262.1

一個基因組會有多個基因,下面是它的基因組結構:

圖片

2

提取單個基因CDS

main.py

1from BaimoTools import BaimoTools23gb_file = f"res/genbank/SARS-CoV-2.gb"4fasta_file = f"out/output_s.fasta"5baimoTools = BaimoTools(gb_file, fasta_file)6# baimoTools.extract_cds('S')

輸出檔案 output_s.fasta,分別提取到兩個基因組的 S 基因 CDS 區域:

3

提取多個基因CDS

main.py

1from BaimoTools import BaimoTools23gb_file = f"res/genbank/SARS-CoV-2.gb"4fasta_file = f"out/output_s_m_orf10.fasta"5baimoTools = BaimoTools(gb_file, fasta_file)6baimoTools.extract_cds(['S', 'M', 'ORF10'])

輸出檔案 output_s_m_orf10.fasta,分別提取到兩個基因組的 S,M,ORF10 基因 CDS 區域::

4

提取全部基因CDS

main.py

1from BaimoTools import BaimoTools23gb_file = f"res/genbank/SARS-CoV-2.gb"4fasta_file = f"out/output_s.fasta"5baimoTools = BaimoTools(gb_file, fasta_file)6# baimoTools.extract_cds("")

輸出檔案 output_all.fasta,分別提取到兩個基因組的全部基因 CDS 區域:

圖片

下一步更新其他基因特徵提取,及格式轉換功能。

本作品採用《CC 協議》,轉載必須註明作者和本文連結

相關文章