前言

CRISPR因为具有NHEJ效应，在修复过程中会造成一定长度的突变，这种突变包括小片段插入(insertion),小片段缺失(deletion),单碱基突变(substitution)。

为了有效准确的分析突变情况，我们需要对测序序列和参考序列进行全局比对，以期望发现不同种类的突变情况

为了加速双序列比对，我们采用cython模块来加速，这里就直接采用别人写好的开源程序Github: brentp / align

  1. git clone https://github.com/brentp/align.git
  2. cd align
  3. python3 setup.py build
  4. python3 setup.py install

这个模块的用法也很简单，但要注意该模块包含python和cython两种比对，当然cython的速度会更快

  1. #for python align
  2. from align.align import aligner as pyaligner
  3. #for cython align
  4. from align.calign import aligner as caligner
  5. from align.matrix import DNAFULL

序列比对有三种方法，分别是global,glogal_cfe,glocal分别针对全局比对，全局长序列比对以及局部比对。

下面直接就是主函数

  1. from align.calign import aligner
  2. from align.matrix import DNAFULL
  3. import re
  4. from collections import Counter
  5. import argparse
  7. parser = argparse.ArgumentParser(
  8.     description='This script was used to calculate indel of 37bp trap library')
  9. parser.add_argument(
  10.     '-r', '--ref', help='Please input your reference file')
  11. parser.add_argument(
  12.     '-f', '--fasta', help='Please input your filter sequence.')
  13. parser.add_argument('-e', '--edits', help='Please input your edits style. eg: insertion;deletion;sub')
  14. parser.add_argument('-m', '--methods', help='Please input your select method. eg: size;pos;ss')
  15. args = parser.parse_args()
  17. ##seq1 is ref
  20. def find_insertion(str1, str2):
  21.     indel = []
  22.     j = 0
  23.     for i in re.finditer("-+", str1):
  24.         pos = i.start() + 1 - j
  25.         size = i.end() - i.start()
  26.         ss = str2[i.start(): i.end()]
  27.         j = j + size
  28.         indel.append((pos, size, ss))
  29.     return indel
  31. def find_deletion(str1, str2):
  32.     indel = []
  33.     new_s1 = str1.replace("-", "")
  34.     new_s2 = "".join([str2[i] for i in range(len(str1)) if str1[i] != "-"])
  35.     for i in re.finditer("-+", new_s2):
  36.         pos = i.start() + 1
  37.         size = i.end() - i.start()
  38.         ss = new_s1[i.start(): i.end()]
  39.         indel.append((pos, size, ss))
  40.     return indel
  43. def find_sub(str1, str2):
  44.     sub = []
  45.     new_s1 = str1.replace("-", "")
  46.     new_s2 = "".join([str2[i] for i in range(len(str1)) if str1[i] != "-"])
  47.     for i in range(len(new_s1)):
  48.         if new_s1[i] != new_s2[i] and new_s2[i] != "-":
  49.             pos = i + 1
  50.             size = 1
  51.             ss = new_s2[i]
  52.             sub.append((pos, size, ss))
  53.         else:
  54.             continue
  55.     return sub
  58. if args.ref:
  59.     lb_lst = []
  60.     rf_lst = []
  61.     with open(args.ref, "r") as f:
  62.         for line in f:
  63.             label, ref_seq = line.split(",")
  64.             lb_lst.append(label.strip())
  65.             rf_lst.append(ref_seq.strip("\n"))
  66. else:
  67.     print("No ref input!")
  68.     exit(0)
  70. if args.fasta:
  71.     fasta = {}
  73.     with open(args.fasta, 'r') as f:
  74.         for line in f:
  75.             if line.startswith(">"):
  76.                 ss_lst = []
  77.                 name = line.strip("\n").replace(">", "")
  78.             elif line != "\n":
  79.                 ss_lst.append(line.strip("\n"))
  80.                 fasta[name] = ss_lst
  81.             else:
  82.                 continue
  83. else:
  84.     print("No sequence input!")
  85.     exit(0)
  87. haha_lst = []
  88. # calculate insertion for every site
  89. for lb, rf in zip(lb_lst, rf_lst):
  90.     for val in fasta[lb]:
  91.         ag = aligner(val, rf, method="global_cfe", matrix=DNAFULL)
  92.         bidui = ag.pop()
  93.         ctrl = bidui.seq2.decode()
  94.         case = bidui.seq1.decode()
  95.         if args.edits == "insertion":
  96.             haha_lst.extend(find_insertion(ctrl, case))
  97.         elif args.edits == "deletion":
  98.             haha_lst.extend(find_deletion(ctrl, case))
  99.         elif args.edits == "sub":
  100.             haha_lst.extend(find_sub(ctrl, case))
  101.         else:
  102.             print("please input correct call methods!eg:insertion;deletion;sub")
  104. if len(haha_lst) > 0:
  105.     if args.methods == "pos":
  106.         new_lst = [haha[0] for haha in haha_lst if len(haha) == 3]
  107.         fin = Counter(new_lst)
  108.         print("{} in every site is:".format(args.edits))
  109.         for i in range(37):
  110.             if fin.get((i+1), 0):
  111.                 print("{}\t{}".format((i + 1), fin[i + 1]))
  112.             else:
  113.                 print("{}\t{}".format((i+1), 0))
  114.     elif args.methods == "size":
  115.         new_lst = [haha[1] for haha in haha_lst if len(haha) == 3]
  116.         fin = Counter(new_lst)
  117.         print("{} in every size is:".format(args.edits))
  118.         for i in range(1, 10):
  119.             if fin.get((i), 0):
  120.                 print("{}\t{}".format((i), fin[i]))
  121.             else:
  122.                 print("{}\t{}".format((i), 0))
  123.     elif args.methods == "ss":
  124.         new_lst = [haha[2] for haha in haha_lst if len(haha) == 3]
  125.         fin = Counter(new_lst)
  126.         print("{} in every sequence is:".format(args.edits))
  127.         for key, vals in fin.items():
  128.             print("{}\t{}".format(key, vals))
  129. else:
  130.     print("No {}".format(args.edits))

用法

如果想知道fasta序列中相对于ref序列，分别在1-37bp上的插入是多少，可以键入一下参数

  1. python3 cal_indel.py -r ref.txt -f fasta.fa -e insertion -m pos
  1. #cat ref.txt
  2. 1,AGAAGTGGAGCTGCAGCTGCAGGCAGCTCCCGGATCC
  3. 2,AGTGGCTGAGTATGATCAGTGTCCAGTGTCTGGCCCA

参考序列是两列文件，第一列代表序号，第二列代表参考序列

  1. #cat fasta.fa
  2. >1
  3. AGAAGTGGAGCTGCAGCTGCAGGCAGCCTCCCGGATCC
  4. AGAAGTGGAGCTGCAGATCC
  5. AGAAGTGGAGCTGCAGATCC
  6. AGAAGTGGAGCT
  7. AGAAGTGGAGCTGCAGCTGCAGGCAGCCTCCCGGATCC
  8. AGAAGTGGAGCTGCAGCTGCAGGCAGAGCTCCTGGATCC
  9. AGAAGTGGAGCTGCAGCTGCAGGCAGAGCTCCTGGATCC
  10. AGAAGTGGAGCT
  11. AGAAGTGGAGCTGCAGCTGCAGGCTGGCCCGGATCC
  12. >2
  13. AGTGGCTGAGTATGATCAGTGTCCAGTCTGGCCCA
  14. AGTGGCTGAGTATGATCAGTGTCTGGCCCA
  15. AGTGGCTGAGTATGATCAGCCCA
  16. AGTGGCTGAGTATGATCAGTGTCCA
  17. AGTGGCTGAGTATGATCAGTGTCCAGTGTCTGGCCCA
  18. AGTGGCTGAGTATGATCAGTGTGTCTGGCCCA
  19. AGTGGCTGAATATGATCAGTGTCCAGTTGTCTGGCCCA
  20. AGTGGCTGAGTATGATCAGTGTCCAG
  21. AGTGGCTGAGTATGATCAGTGTGTCTGGCCCA
  22. AATGGCTGAGTATGATCAGTGTCCAGATCTGGCCCA
  23. AGTGGCTGAATAGTCTGGCCCA

输入序列为序列名称为参考序列对应的编号，而后面跟着以换行符分隔的小片段序列

-e参数可以分别输入insertion,deletion,sub,分别对应插入，缺失以及单碱基替换
-m参数可以分别输入size,pos,ss分别对应按照大小，位置和序列进行统计

注意：如果报global_cfe的错误，就是因为比对序列太多了，更换比对参数为global即可！