双序列比对后统计插入与删除（indels）的数量

• 叶绿体基因组类的文章偶尔会有这类分析，但是在方法部分的描述自己还没有看得明白。比如 Comparative chloroplast genomics between the invasive weed Mikania micrantha and its indigenous congener Mikania cordata: Structure variation, identification of highly divergent regions, divergence time estimation, and phylogenetic analysis 这篇文章中方法部分提到 Base substitution events and insertion/deletion (indel) were identified through DnaSP v5.0. 结果部分提到 Using the M.cordata cp genome as reference, we detected 58 indels between M.micrantha and M.cordata. 自己的做法是使用mafft对两个序列进行比对，然后统计序列中出现 “-” 的位置进行统计； 本来想尝试用文章中的数据来验证自己的想法是否正确，但是文章中提到的叶绿体基因组序列暂时还不能下载。根据自己的想法自己的数据使用DnaSP软件统计了插入删除的数量，结果统计到了500多。正常同一个属的叶绿体基因组差别会有这么大吗？想到了两个原因：一是自己最初组装的就有问题；二是分析插入删除的方法是错误的。如果是问题一的话自己的麻烦可就大了，但是如何评估叶绿体基因组组装结果的准确性这方面自己还不太明白（自己的叶绿体基因组是利用全基因组重测序数据使用chloroExtractor软件组装来的，回头仔细看一下这款软件对应的文章）；
• 在假设组装结果没有问题并且自己关于插入与删除分析方法正确的前提下，使用DnaSP软件只能够统计到插入删除的数量，不能够给出插入与删除位置的序列信息，或者说自己还不知道如何利用DnaSP软件得到上述信息；自己尝试写了一个简单的python脚本得到插入与删除位置的序列信息，简单记录

#导入re和SeqIO模块，re模块用来写正则表达式，SeqIO模块用来解析fasta文件
import re
from Bio import SeqIO
#分别将比对后的两条序列存入到列表中
seq = []
for record in SeqIO.parse("aligned.fasta","fasta"):
  seq.append(str(record.seq))
#使用正则表达式匹配序列中出现“-”的位置，然后输出到文件中
motif = "-+"
regexp = re.compile(motif)
#统计第一条序列中是 “-” 的位置
match1 = regexp.finditer(seq[0])
fw1 = open("indels_1.txt","w")
i = 0
for s in match1:
  i = i + 1
  location = list(s.span())
  fw1.write(str(i)+"\n")
  fw1.write(seq[0][location[0]:location[1]]+"\n")
  fw1.write(seq[1][location[0]:location[1]]+"\n")
fw1.close()
#统计第二条序列中是 “-” 的位置，输出到文件中
match2 = regexp.finditer(seq[1])
fw2 = open("indels_2.txt","w")
i = 0
for s in match2:
  i = i + 1
  location = list(s.span())
  fw2.write(str(i)+"\n")
  fw2.write(seq[0][location[0]:location[1]]+"\n")
  fw2.write(seq[1][location[0]:location[1]]+"\n")
fw2.close()

• 还剩一个问题：在基于序列和方法都没有问题的前提下， 如何获得插入与删除在原序列中的位置（位于基因内或者基因间隔区，位于内含子还是外显子），这个得仔细想一下。

未完待续。。。。。。