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

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

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

欢迎大家关注我的公众号
小明的数据分析笔记本

小明的数据分析笔记本

中国加油！武汉加油！