python脚本提取叶绿体基因组的大小单拷贝区、反向重复区

用户7010445

发布于 2020-03-03 15:16:37

1.7K0

发布于 2020-03-03 15:16:37

叶绿体基因组结构保守，包含四部分结构：大单拷贝区、小单拷贝区、两个反向重复区。叶绿体基因组类的文章通常会计算这四个区域的变异位点。那么第一步便是从完整的叶绿体基因组的序列中分别将这四个区域提取出来，然后比对计算。本篇文章记录提取这四个区域用到的python脚本

第一步：利用叶绿体基因组的fasta文件得到反向重复区的位置信息

叶绿体基因组类的文章通常是我们自己做几个，然后结合已经发表的数据做分析。已经公布在NCBI的叶绿体基因组中通常没有反向重复区的信息。这个时候就需要我们自己重新注释。注释用到的是在线工具GeSeq https://chlorobox.mpimp-golm.mpg.de/geseq.html

我这里用4个苹果属的叶绿体基因组做演示，序列号分别是：

NC_031163
NC_035625
NC_035671
NC_036368

将4个fasta文件上传到GeSeq，点击提交就可以了

image.png

很快就可以运行完，下载标注的文件用于后续分析

这个文件里包含里两个反向重复区的位置信息

image.png

提取脚本

import os
import sys
from Bio import SeqIO

inputFile = sys.argv[1]

fwLSC = open("LSC_region.fasta",'w')
fwIR = open("IR_region.fasta",'w')
fwSSC = open("SSC_region.fasta",'w')

for rec in SeqIO.parse(inputFile,'gb'):
    IR_pos = []
    for feature in rec.features:
        if feature.type == "repeat_region":
            for part in feature.location.parts:
                IR_pos.append(part.start)
                IR_pos.append(part.end)
    if len(rec.seq) == max(IR_pos):
        print(rec.id," 可以进行下一步分析！")
        IR_pos.sort()
        LSC = rec.seq[0:(IR_pos[0]-1)]
        IR = rec.seq[(IR_pos[0]-1):IR_pos[1]]
        SSC = rec.seq[IR_pos[1]:IR_pos[2]]
        fwLSC.write(">%s\n%s\n"%(rec.id,str(LSC)))
        fwIR.write(">%s\n%s\n"%(rec.id,str(IR)))
        fwSSC.write(">%s\n%s\n"%(rec.id,str(SSC)))
    else:
        fwLSC.close()
        fwIR.close()
        fwSSC.close()
        os.remove("LSC_region.fasta")
        os.remove("IR_region.fasta")
        os.remove("SSC_region.fasta")
        print(rec.id," 需要调整IR区域的相对位置！")
        break
if "LSC_region.fasta" in os.listdir():
    print("结果文件分别是：","LSC_region.fasta ","SSC_region.fasta ","IR_region.fasta ")
else:
    print("调整后重新注释再来提取！")

运行脚本

 python .\extract_LSC_SSC_IRs_from_cp_genome.py .\GeSeqJob-20200205-103624_GLOBAL_multi-GenBank.gbff

会提示我

NC_031163.  可以进行下一步分析！
NC_036368.  需要调整IR区域的相对位置！
调整后重新注释再来提取！

这是因为这条序列的反向重复区位置和通常的不一样

image.png

因为叶绿体基因组是环状的，放到文件里存储你可以选择任意一个碱基作为开始的第一个，叶绿体基因组通常是大单拷贝区的第一个碱基作为起始，但是这条序列不符合普遍情况，我们需要将序列起始的31个碱基放到序列的尾部

用到的脚本

import sys
from Bio import SeqIO

inputFile = sys.argv[1]
pos = int(sys.argv[2])

for rec in SeqIO.parse(inputFile,'fasta'):
    fw = open(rec.id+"_1.fasta",'w')
    seq_1 = rec.seq[0:pos]
    seq_2 = rec.seq[pos:]
    fw.write(">%s\n%s\n%s\n"%(rec.id,str(seq_2),str(seq_1)))

fw.close()

使用方法

python .\adjust_IR_pos.py .\NC_036368.fasta 31

然后利用输出文件NC_036368.1_1.fasta重新去注释注释完以后再来运行第一个脚本

python .\extract_LSC_SSC_IRs_from_cp_genome.py .\GeSeqJob-20200205-121603_GLOBAL_multi-GenBank.gbff

得到结果

NC_031163.  可以进行下一步分析！
NC_035625.  可以进行下一步分析！
NC_035671.  可以进行下一步分析！
NC_036368.  可以进行下一步分析！
结果文件分别是： LSC_region.fasta  SSC_region.fasta  IR_region.fasta

如果需要以上脚本，在我的公众号留言就可以了！！

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2020-02-28，如有侵权请联系 cloudcommunity@tencent.com 删除

腾讯云开发者社区

本文分享自小明的数据分析笔记本微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

腾讯云开发者社区

登录后参与评论

0 条评论

热度