python腳本提取葉綠體基因組的大小單拷貝區、反向重複區
- 2020 年 3 月 3 日
- 筆記
葉綠體基因組結構保守，包含四部分結構：大單拷貝區、小單拷貝區、兩個反向重複區。葉綠體基因組類的文章通常會計算這四個區域的變異位點。那麼第一步便是從完整的葉綠體基因組的序列中分別將這四個區域提取出來，然後比對計算。 本篇文章記錄提取這四個區域用到的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(">%sn%sn"%(rec.id,str(LSC))) fwIR.write(">%sn%sn"%(rec.id,str(IR))) fwSSC.write(">%sn%sn"%(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(">%sn%sn%sn"%(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
如果需要以上腳本，在我的公眾號留言就可以了！！
