精华内容
下载资源
问答
  • Python读取fasta文件

    千次阅读 2019-04-23 22:13:27
    AMPs.fasta文件格式: >AP00001 |antibacterial |anticancer/tumor |antifungal GLWSKIKEVGKEAAKAAAKAAGKAALGAVSEAV >AP00004 |antibacterial |antifungal ...

    AMPs.fasta文件格式:

    >AP00001 |antibacterial |anticancer/tumor |antifungal
    GLWSKIKEVGKEAAKAAAKAAGKAALGAVSEAV
    
    >AP00004 |antibacterial |antifungal
    NLCERASLTWTGNCGNTGHCDTQCRNWESAKHGACHKRGNWKCFCYFDC
    
    >AP00005 |antibacterial
    VFIDILDKVENAIHNAAQVGIGFAKPFEKLINPK
    
    >AP00006 |antibacterial
    GNNRPVYIPQPRPPHPRI
    
    >AP00008 |antibacterial
    RLCRIVVIRVCR
    
    >AP00020 |antibacterial |anticancer/tumor |antifungal
    GLFDIVKKIAGHIAGSI
    
    >AP00026 |antibacterial |anticancer/tumor |antifungal |anti-HIV |antiviral
    FKCRRWQWRMKKLGAPSITCVRRAF
    
    >AP00027 |antibacterial
    ITPATPFTPAIITEITAAVIA
    
    >AP00035 |antibacterial |anticancer/tumor
    KSSAYSLQMGATAIKQVKKLFKKWGW
    
    >AP00050 |antibacterial
    GIGASILSAGKSALKGLAKGLAEHFAN

     

    代码任务:读取序列并统计序列长度分布

    import seaborn as sns
    import pandas as pd
    import matplotlib.pyplot as plt
    
    sns.set_style('dark')
    
    train_pos_data_path = r"./AMPs.fasta"
    train_neg_data_path = r"./nonAMPs.fasta"
    test_pos_data_path = r'./amp920.fasta'
    test_neg_data_path = r'./nonamp920.fasta'
    
    def read_from_file_with_enter(filename):
        fr = open(filename,'r')
        sample = ""
        samples = []
        for line in fr:
            if line.startswith('>'):
                sample = ""
                continue
            if line.startswith('\n'):
                samples.append(sample)
                continue
            sample += line[:-1]
        return samples
    
    def read_from_file(filename):
        fr = open(filename, 'r')
        sample = ""
        samples = []
        for line in fr:
            if line.startswith('>'):
                if sample != "":
                    samples.append(sample)
                    sample=""
                continue
            sample+=line[:-1]
        return samples
    
    def statistics_length(samples):
        lengths = []
        for sample in samples:
            lengths.append(len(sample))
        return pd.DataFrame(lengths,columns=['Length'])
    
    
    train_pos_data = read_from_file_with_enter(train_pos_data_path)
    train_neg_data = read_from_file(train_neg_data_path)
    test_pos_data = read_from_file(test_pos_data_path)
    test_neg_data = read_from_file(test_neg_data_path)
    
    
    train_pos_len = statistics_length(train_pos_data)
    train_neg_len = statistics_length(train_neg_data)
    
    plt.figure(figsize=(6,4))
    g = sns.distplot(a=train_pos_len.Length,label="Train_Pos",kde=False,bins=20)
    g = sns.distplot(a=train_neg_len.Length,label="Train_Neg",kde=False,bins=20)
    plt.ylabel("Number")
    plt.title("Length of sample")
    plt.legend(loc='best')
    
    plt.figure(figsize=(6,4))
    g = sns.kdeplot(data=train_pos_len.Length,shade=True)
    g = sns.kdeplot(data=train_neg_len.Length,shade=True)
    plt.title("Kernel density estimation")
    
    plt.show(g)

     

    效果展示

    直方图:

     

    核密度估计图:

    展开全文
  • 利用Python读取fasta文件并进行一系列操作(上) 概述 语言:python3.8 模块:pysam collections 可选:jupyter 整体思路:将fasta格式的基因原始数据处理为方便读写的txt格式并进行操作 步骤: 获取自己的fasta...

    利用Python读取fasta文件并进行一系列操作(上)


    概述


    语言:python3.8
    模块:pysam collections
    可选:jupyter
    整体思路:将fasta格式的基因原始数据处理为方便读写的txt格式并进行操作

    步骤:


    1. 获取自己的fasta文件(这里我将从NCBI上下载人类的ABO基因参考序列的fasta文件为例)
      在这里插入图片描述
    2. 利用pysam模块的FastaFile函数读取fasta,之后即可获取fasta的基本信息:filename 文件名,references 染色体编号(因为这里我下载的是ABO基因的fasta,所以不是染色体编号),nreferences 染色体数,lengths 每条染色体长度(这个非常重要!)
    import pysam as sam
    # 读取fasta
    fasta = sam.FastaFile('ABO_sequence.fasta')
    # 获取染色体编号
    fasta.references
    # 获取文件名
    fasta.filename
    # 获取染色体数
    fasta.nreferences
    # 获取每条染色体长
    fasta.lengths
    

    返回结果如下:

    ['NG_006669.2']
    
    b'ABO_sequence.fasta'
    
    1
    
    [42144]
    
    1. 通过获取输出的长度,我们知道了这个ABO基因的长度值为42144,我们再利用fetch函数获取ABO基因碱基序列(注意索引是从0开始)
    data = fasta.fetch('NG_006669.2', 0, 42144)
    

    输出结果如下图:
    在这里插入图片描述

    1. 我们的最终目的是要输出txt,所以就需要将我们的碱基序列写到txt上
    with open('ABO_seq.txt', 'w') as f:
        f.write(data)
    
    1. 我们可以利用我们输出的txt来做点事情,例如看看这个序列中4种碱基数量,先读取文件,再处理为列表,利用collections模块的Counter函数来输出结果
    from collections import Counter
    
    with open('ABO_seq.txt', 'r') as f:
        seq = f.read()
        seq.split()
        print(Counter(seq))
    

    结果展示:

    在这里插入图片描述
    在这里插入图片描述

    下面总结一下将fasta输出为txt的代码如下:


    import pysam as sam
    # 读取fasta
    fasta = sam.FastaFile('ABO_sequence.fasta')
    # 获取每条染色体长
    fasta.lengths
    # 获取指定的碱基序列
    data = fasta.fetch('NG_006669.2', 0, 42144)
    # 将序列输出为txt文件
    with open('ABO_seq.txt', 'w') as f:
        f.write(data)
    

    更多操作将在下篇解锁!

    展开全文
  • 利用Python读取fasta文件并进行一系列操作(二) 概览: 本节目标:通过上一节所输出的txt输出ABO蛋白以及ABO基因的外显子fasta文件 语言: python3.8 模块:biopython ssl 可选:jupyter 整体思路:通过ncbi获取...

    利用Python读取fasta文件并进行一系列操作(二)


    概览:


    本节目标:通过上一节所输出的txt输出ABO蛋白以及ABO基因的外显子fasta文件
    语言: python3.8
    模块:biopython ssl
    可选:jupyter
    整体思路:通过ncbi获取ABO基因外显子位置(爬虫),并读取txt,根据位置信息获取外显子序列,再通过外显子
    序列输出mRNA序列以及蛋白序列
    前排提示:本教程不管生物,有知识盲区自己补

    步骤:


    1. 设置ssl,要不debug后总会有bug
    2. 从这步开始调用biopython,设置邮箱(不要瞎填!!瞎填还不如不填!!!)
    3. 设置有关搜索ncbi的函数以及参数,获取人类ABO基因的mRNA相关资料所对应id
    4. 获取人类ABO基因的mRNA相关资料
    5. 将获取的所有资料进行处理,得到外显子位置
    6. 获取外显子序列,并将所有外显子数据写入新创建的fasta文件中

    代码:


    from Bio import Entrez
    from Bio.Seq import Seq
    from Bio.SeqRecord import SeqRecord
    from Bio import SeqIO
    import ssl
    
    # 设置ssl,要不debug后总会有bug
    ssl._create_default_https_context = ssl._create_unverified_context
    # 设置邮箱(不要瞎填!!瞎填还不如不填!!!)
    Entrez.email = 'yhl030410@163.com'
    # 设置有关搜索ncbi的函数以及参数,获取人类ABO基因的mRNA相关资料所对应id
    handle = Entrez.esearch(db = 'nucleotide',term = 'Homo sapiens ABO, alpha 1-3-N-acetylgalactosaminyltransferase and alpha 1-3-galactosyltransferase (ABO), mRNA')
    record = Entrez.read(handle)
    # 获取人类ABO基因的mRNA相关资料
    id = record["IdList"][0]
    handle = Entrez.efetch(db = 'nucleotide', id = id, rettype = 'gb', retmode = 'text')
    data = handle.read()
    # 将获取的所有资料进行处理,得到外显子位置
    exons_l = []
    for i in data.split('\n'):
        if 'exon' in i:
            exons_l.append(i.split())
    exons_num = {}
    for i in range(len(exons_l)):
        exons_num[f'exon{i+1}'] = exons_l[i][-1].split('..')
    # 获取外显子序列,并将所有外显子数据写入新创建的fasta文件中
    iter_exon = exons_num.keys()
    with open('ABO_seq.txt', 'r') as f:
        seq = f.read()
        data = []
        for i in iter_exon:
            start = int(exons_num[i][0]) - 1
            end = int(exons_num[i][-1]) - 1
            exon_seq = seq[start:end]
            seqs = SeqRecord(Seq(exon_seq), id = i, description = 'Homo sapiens ABO exons')
            data.append(seqs)
        SeqIO.write(data, 'ABO_exons.fa', 'fasta')
    

    获取蛋白序列继续往下看:

    接上面步骤:


    1. 将获取的外显子序列用Seq函数转为biopython类型
    2. 将外显子倒序并转录
    3. 将转录出的mRNA序列翻译为蛋白质序列

    代码:


    #获取mRNA序列
        seq = []
        for i in iter_exon:
            start = int(exons_num[i][0]) - 1
            end = int(exons_num[i][-1]) - 1
            exon_seq = seq[start:end]
            seq.append(exon_seq)
        seq = Seq(seq.join())
        mrna = seq.reverse_complement().transcribe()
        protein = mrna.translate()
    
    展开全文
  • 利用Python读取fasta文件并进行一系列操作(三) 概述: 本节目标:计算智人与猩猩ABO基因的相对熵 语言:python3.8 模块:pysam, scipy 整体思路:先计算出序列中“AG”“CT”“AC”“AT”“GC”“GT”六种组合...

    利用Python读取fasta文件并进行一系列操作(三)


    概述:


    本节目标:计算智人与猩猩ABO基因的相对熵
    语言:python3.8
    模块:pysam scipy
    整体思路:先计算出序列中AG CT AC AT GC GT六种组合序列所占比,再计算相对熵

    步骤:


    1. 利用pysam模块分别读取智人ABO基因所有序列和猩猩ABO基因所有序列
    import pysam as sam
    from scipy import stats
    hfasta = sam.FastaFile('homo_ABO.fasta')
    cfasta = sam.FastaFile('chimp_ABO.fasta')
    c_all = cfasta.fetch('NC_036888.1:c105542142-105516107', 0, 26036)
    h_all = hfasta.fetch('NG_006669.2', 0, 42144)
    
    1. 利用split函数切序列并用len函数返回结果-1,分别得出六种组合序列数量,再求占比
    cAG = len(c_all.split('AG')) - 1
    cCT = len(c_all.split('CT')) - 1
    cAC = len(c_all.split('AC')) - 1
    cAT = len(c_all.split('AT')) - 1
    cGC = len(c_all.split('GC')) - 1
    cCT = len(c_all.split('CT')) - 1
    p = cCT/sum([cAG,cCT,cAC,cAT,cGC,cCT])
    h_all = hfasta.fetch('NG_006669.2', 0, 42144)
    hAG = len(h_all.split('AG')) - 1
    hCT = len(h_all.split('CT')) - 1
    hAC = len(h_all.split('AC')) - 1
    hAT = len(h_all.split('AT')) - 1
    hGC = len(h_all.split('GC')) - 1
    hCT = len(h_all.split('CT')) - 1
    q = hCT/sum([hAG,hCT,hAC,hAT,hGC,hCT])
    
    1. 利用scipy模块的entropy函数求相对熵
    result = stats.entropy([p,q])
    print(result)
    

    结果:


    0.6930868205144018
    
    展开全文
  • 使用python读取和分析fasta文件

    千次阅读 2019-12-05 20:28:49
    分享一些处理fasta文件python函数
  • [笔记]pythonFASTA文件的处理

    万次阅读 多人点赞 2018-04-27 19:17:06
    这学期选了生信的选修课—perl/python在生物信息学中的应用把结课作业的代码整理出来主要是pythonFASTA文件读取和数据处理FASTA文件数据处理FASTA文件读取:只含一个基因序列将FASTA文件的基因序列读取到一个...
  • 使用Python计算fasta文件的序列长度

    千次阅读 2018-11-30 15:29:38
    使用Python计算fasta文件的序列长度 在这里插入代码片使用Python计算fasta文件的序列长度 #!/usr/bin/python #-- coding:utf-8 -- import sys f = open(sys.argv[1],‘r’) out = open(sys.argv[2],‘w’) def chr_...
  • 使用Pythonfasta格式的序列进行基本信息统计预期设计输出文件中包括fasta文件名,序列长度,GC含量以及ATCG各自的含量。Python脚本编辑使用的文件test.fastastat.py输入 sys模块#!/usr/bin/env pythonimport sys从...
  • 我正在尝试读取FASTA文件,然后查找特定的motif(string)并打印出序列和次数. AFASTA file只是一系列序列(字符串),以标题行开头,标题或新序列的开头是“>”.在标题之后的一个新行中是字母序列.我没有完成代码但到...
  • 于是决定只挑基因组前十条染色体拿来练习(所以需要从基因组文件里选取序列,尝试自己用python写脚本处理)。python之前因为实在不聪明也就刚学到字典部分,借鉴了网上的一些内容整理了一下,写了几种情况关于利用...
  • 在如果是,请将accessorID更改为:accessorID = accessorIDWithArrow[1:5]一些方法可以让这更像Python:使用集合而不是字典,使用strip()而不是切片来删除换行符,并使用生成器表达式来构建集合^{pr2}$使用True和{}...
  • 最近,用Python脚本提取,在基因号已知,位置已知条件下,相对应位置的基因序列时发现,这样很简单但是很实用的脚本,在网上却比较难找。而且,能被找到的脚本,相对于具有初级编程能力的人而言,有点难。本人写了相...
  • Python蛋白质序列文件(.fasta)拆分

    千次阅读 2020-10-04 16:48:12
    主要记录在python应用中遇到的一些小问题,小编平时做的主要是DNA结合蛋白的识别与预测,需要借助多种软件进行特征提取,如PSI——BLAST,PSIpred等等,小编这次要解决的问题就是小编拿到的原始序列数据为一个fasta...
  • python使用SeqIO #pip install biopython from Bio import SeqIO with open("./clean_total_genomic.fna",'r') as fq: for record in SeqIO.parse(fq,'fasta'): print(record.seq) 自己写文件处理 ...
  • 1、需要将下面的fasta文件进行一个整理,将序列单行输出 方法一 f1 = open('test1.fa','r').readlines()#需要整理的文件 f2 = open('2.fasta','w')#整理之后的文件 for i in f1: if i.startswith('>'): ...
  • 设计需求统计fasta文件中多条序列信息,设计目标效果:图片.png将结果输入到csv格式的表格中,因为csv格式表格用,分隔数据。脚本使用argparse模块,提示输入数据。import argparseparser = argparse.ArgumentParser...
  • 最近用Bio.SeqIO模块进行读取fasta文件到字典中的时候发现一个问题,如果你的fasta文件>开头的那一行header中含有空格的话,该行内容以键存到字典里,这个header会被从第一个空格的地方截断,比如原本的文件是...
  • out.fas提取,但是这次遇到了一个大文件,用grep就太费时了,然后又试了一下TBtools的提取序列功能,发现时间也很长,所以就写了个python。提取将近100万条reads耗时也就需要10s左右 #!/usr/bin/python # -*- ...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 552
精华内容 220
关键字:

python读取fasta文件

python 订阅