前言

Crispr基因编辑正越来越广泛的应用在各个方面，包括科研，医疗等等，针对经过筛选的药物靶向基因设计gRNA，使其由原始的基因序列突变为终止密码子，从而无法表达蛋白，进而治疗疾病或者抵抗药物

算法

Gene 沉默效率计算方法 - 图1

初步过滤：

  1. 将12ktrap中，N1-N20内C>T成功的trap提取出来(所有编辑后的reads，不做trap的去重)，定义为dataset#1
  2. 将dataset#1分成2部分，一部分是gRNA在基因的正义链（+）上，定义为dataset#1（+），另一部分在反义链（-）上，定义为dataset#1（-）；

正义链的情况（CAG/CGA/CAA > TAG/TGA/TAA）：

  1. 找出dataset#1（+）N1-N20中 C>T成功的位置，并往后推2个碱基，判断其是否为CAG或CGA或CAA（编辑前） ，满足其一即抽提出来，完成后定义为dataset#2（+），同时标记出突变位点在gRNA中的位置，譬如N=5, N=6等；
  2. 将dataset#2（+）中的gRNA（不包含我们人工添加的第一个碱基G）mapping回对应基因的CDS中，数出从ATG（翻译起始位置）的A到gRNA第一个碱基的碱基数M，随后算出 (M-1+N-1)/3, 判断结果是否为整数，如果是整数，判断该基因被stop成功，输出”TRAP ID + gRNA sequence + strand (+/-) + C>T position + C>T efficiency”;
  3. 举例说明4）中算法：如N=6, M=11, 则说明突变位点的C距离ATG位置是10+5=15个碱基，即15/3=5个氨基酸，到突变位点（如CAG）恰好是3的整数倍，该基因被stop成功；

反义链的情况(CCA > TCA (TGA) / CTA(TAG) / TTA (TAA))：

  1. 找出dataset#1（-）N1-N20中 C>T成功的位置，并判断其是否在CCA这样的motif中（编辑前），如果在，输出trap放到dataset#2（-）中，并标记出CCA motif中A所在gRNA中的位置N，譬如N=5, N=6等；
  2. 将dataset#2（-）中的gRNA提取出来（此时gRNA的序列应该是和CDS正义链互补的），数出从ATG（翻译起始位置）的A到gRNA第一个碱基的碱基数M （注意此时gRNA为CCNN20，gRNA第一个碱基在最右端），随后算出(M-N)/3是否为整数，如果是整数，判断该基因被stop成功，输出”TRAP ID + gRNA sequence + strand (+/-) + C>T position + C>T efficiency”;

需求是：

  1. 获取每条TRAP序列所对应的stop效率，效率的计算方法是针对一条TRAP来说，至少有一个位点能够stop成功，那么这条对应的read为成功一次，一次类推，既不能重复计数，也不能少计数
  2. 获取每个位点对应的stop效率，如上原理所述
  3. 获取基因的stop率：stop-gene/all-gene

我使用的代码

代码一：此代码用于输出所有通过的stop reads并计算每条TRAP的stop效率

  1. # -*- coding=utf-8 -*-
  2. import re
  5. def complement(seq):
  6.     return seq.translate(str.maketrans('ACGT', 'TGCA'))
  7. # Obtain reverse complementary sequence
  10. def revcomp(seq):
  11.     return complement(seq)[::-1]
  12. # Input n site, all CDs sequences and trap37bp sequences, and check whether the site can stop
  13. # 1 for success, 0 for failure
  16. def stop_test_for_Z(n, cds_lst, trap):
  17.     new_cds_lst = [cds for cds in cds_lst if trap in cds]
  18.     score_lst = []
  19.     if len(new_cds_lst) > 0:
  20.         for haha in new_cds_lst:
  21.             m = haha.find(trap) + 1 + 10
  22.             score = (m - 1 + n - 1)
  23.             score_lst.append(score)
  24.         if len([new_score for new_score in score_lst if new_score % 3 == 0]) > 0:
  25.             return 1
  26.         else:
  27.             return 0
  28.     else:
  29.         return 0
  32. def stop_test_for_F(n, cds_lst, trap):
  33.     trap = revcomp(trap)
  34.     new_cds_lst = [cds for cds in cds_lst if trap in cds]
  35.     score_lst = []
  36.     if len(new_cds_lst) > 0:
  37.         for haha in new_cds_lst:
  38.             m = haha.find(trap) + 1 + 10 + 23
  39.             score = (m - n)
  40.             score_lst.append(score)
  41.         if len([new_score for new_score in score_lst if new_score % 3 == 0]) > 0:
  42.             return 1
  43.         else:
  44.             return 0
  45.     else:
  46.         return 0
  47. # Read all CDS files and return a dictionary. The key is gene and the value is CDS sequence
  50. def read_cds(fs):
  51.     fasta = {}
  52.     with open(fs, 'r') as f:
  53.         for line in f:
  54.             if line.startswith(">"):
  55.                 seq_lst = []
  56.                 name = line.strip("\n").replace(">", "")
  57.             else:
  58.                 seq_lst.append(line.strip("\n"))
  59.                 fasta[name] = "".join(seq_lst)
  60.     new_fasta = {key: value for key,
  61.                  value in fasta.items() if "unavailable" not in value}
  62.     return new_fasta
  63. # Read all the trap sequences and return a dictionary. The key is trap label and the value is 37bp
  66. def read_seq(fs):
  67.     fasta = {}
  68.     with open(fs, 'r') as f:
  69.         for line in f:
  70.             if line.startswith(">"):
  71.                 seq_lst = []
  72.                 name = line.strip("\n").replace(">", "")
  73.             else:
  74.                 seq_lst.append(line.strip("\n"))
  75.                 fasta[name] = seq_lst
  76.     return fasta
  77. # Read all reference sequences. The key is label and the value is 37bptrap sequence
  80. def read_ref(fs):
  81.     ref = {}
  82.     with open(fs, 'r') as f:
  83.         for line in f:
  84.             if line.startswith(">"):
  85.                 name = line.strip("\n").replace(">", "")
  86.             else:
  87.                 ref[name] = line.strip("\n")[118:155]
  88.     return ref
  89. # Read the corresponding relationship between gene and label and return to the dictionary
  92. def read_duiying(fs):
  93.     duiying = {}
  94.     with open(fs, 'r') as f:
  95.         for line in f:
  96.             lable, gene = line.split("\t")
  97.             duiying[lable.strip("\n")] = gene.strip("\n")
  98.     return duiying
  101. def motif_Z(trap, site):
  102.     if trap[10:30][site: site + 3] in ["CAG", "CGA", "CAA"]:
  103.         return 1
  104.     else:
  105.         return 0
  108. def motif_F(trap, site):
  109.     if "CCA" in trap[10:30][site-2: site + 3]:
  110.         return 1
  111.     else:
  112.         return 0
  115. lable_map_gene = read_duiying("../duiying.txt")
  116. all_cds = read_cds("../mart_export.txt")
  117. all_ref = read_ref("ref.fa")
  118. all_trapseq = read_seq("new508.filter.fa")
  119. eff_lst = []
  120. for name, traps in all_ref.items():
  121.     gene = lable_map_gene[name]
  122.     cds_lst = [vals for keys, vals in all_cds.items() if gene in keys]
  123.     pass_lst = []
  124.     if name in list(all_trapseq.keys()):
  125.         alltrap = len(all_trapseq[name])
  126.         seq_lst = all_trapseq[name]
  127.         for m in re.finditer("C", traps[10:30]):
  128.             if motif_Z(traps, m.start()):
  129.                 if stop_test_for_Z((m.start() + 1), cds_lst, traps):
  130.                     while seq_lst:
  131.                         gg=seq_lst.pop(0)
  132.                         if gg[10:30][m.start()] == "T":
  133.                             pass_lst.append(gg)
  134.         if len(pass_lst) > 0:
  135.             with open("pass_lst_z.fa", 'a') as wocao:
  136.                 wocao.write(">{}\n".format(name))
  137.                 for nima in pass_lst:
  138.                     wocao.write("{}\n".format(nima))
  139.             if alltrap > 0:
  140.                 eff_lst.append([name, len(pass_lst) / alltrap])
  141.             else:
  142.                 print(name)
  144. for name, traps in all_ref.items():
  145.     gene = lable_map_gene[name]
  146.     cds_lst = [vals for keys, vals in all_cds.items() if gene in keys]
  147.     pass_lst = []
  148.     if name in list(all_trapseq.keys()):
  149.         alltrap = len(all_trapseq[name])
  150.         seq_lst2 = all_trapseq[name]
  151.         for m in re.finditer("C", traps[10:30]):
  152.             if motif_F(traps, m.start()):
  153.                 if stop_test_for_F(((traps[10:30][m.start()-2: m.start() + 3]).find("CCA")+1+2), cds_lst, traps):
  154.                     while seq_lst2:
  155.                         gg=seq_lst2.pop(0)
  156.                         if gg[10:30][m.start()] == "T":
  157.                             pass_lst.append(gg)
  158.         if len(pass_lst) > 0:
  159.             with open("pass_lst_f.fa", 'a') as wocao:
  160.                 wocao.write(">{}\n".format(name))
  161.                 for nima in pass_lst:
  162.                     wocao.write("{}\n".format(nima))
  163.             if alltrap > 0:
  164.                 eff_lst.append([name, len(pass_lst) / alltrap])
  165.             else:
  166.                 print(name)
  168. with open("trap_stop_eff.txt", 'a') as l:
  169.     l.write("trap-lable\ttrap-effective\n")
  170.     for lable, effctive in eff_lst:
  171.         l.write("{}\t{}\n".format(lable, effctive))

代码二：该代码用于计算需求2

  1. # -*- coding=utf-8 -*-
  2. import re
  5. def complement(seq):
  6.     return seq.translate(str.maketrans('ACGT', 'TGCA'))
  7. # 获取反向互补序列
  10. def revcomp(seq):
  11.     return complement(seq)[::-1]
  12. # 输入n位点，对应的所有CDS序列以及trap37bp序列，检查该位点能否stop
  13. # 成功返回1，失败返回0
  16. def stop_test_for_Z(n, cds_lst, trap):
  17.     new_cds_lst = [cds for cds in cds_lst if trap in cds]
  18.     score_lst = []
  19.     if len(new_cds_lst) > 0:
  20.         for haha in new_cds_lst:
  21.             m = haha.find(trap) + 1 + 10
  22.             score = (m - 1 + n - 1)
  23.             score_lst.append(score)
  24.         if len([new_score for new_score in score_lst if score % 3 == 0]) > 0:
  25.             return 1
  26.         else:
  27.             return 0
  28.     else:
  29.         return 0
  32. def stop_test_for_F(n, cds_lst, trap):
  33.     trap = revcomp(trap)
  34.     new_cds_lst = [cds for cds in cds_lst if trap in cds]
  35.     score_lst = []
  36.     if len(new_cds_lst) > 0:
  37.         for haha in new_cds_lst:
  38.             m = haha.find(trap) + 1 + 10 + 23
  39.             score = (m - n)
  40.             score_lst.append(score)
  41.         if len([new_score for new_score in score_lst if score % 3 == 0]) > 0:
  42.             return 1
  43.         else:
  44.             return 0
  45.     else:
  46.         return 0
  47. # 读取所有的cds文件返回一个字典，键为gene，值为cds序列
  50. def read_cds(fs):
  51.     fasta = {}
  52.     with open(fs, 'r') as f:
  53.         for line in f:
  54.             if line.startswith(">"):
  55.                 seq_lst = []
  56.                 name = line.strip("\n").replace(">", "")
  57.             else:
  58.                 seq_lst.append(line.strip("\n"))
  59.                 fasta[name] = "".join(seq_lst)
  60.     new_fasta = {key: value for key,
  61.                  value in fasta.items() if "unavailable" not in value}
  62.     return new_fasta
  63. # 读取所有的trap序列返回一个字典，键为trap-lable，值为37bp列表
  66. def read_seq(fs):
  67.     fasta = {}
  68.     with open(fs, 'r') as f:
  69.         for line in f:
  70.             if line.startswith(">"):
  71.                 seq_lst = []
  72.                 name = line.strip("\n").replace(">", "")
  73.             else:
  74.                 seq_lst.append(line.strip("\n"))
  75.                 fasta[name] = seq_lst
  76.     return fasta
  77. # 读取所有的参考序列，键为lable,值为37bptrap序列
  80. def read_ref(fs):
  81.     ref = {}
  82.     with open(fs, 'r') as f:
  83.         for line in f:
  84.             if line.startswith(">"):
  85.                 name = line.strip("\n").replace(">", "")
  86.             else:
  87.                 ref[name] = line.strip("\n")[118:155]
  88.     return ref
  89. # 读取基因和lable的对应关系，返回字典
  92. def read_duiying(fs):
  93.     duiying = {}
  94.     with open(fs, 'r') as f:
  95.         for line in f:
  96.             lable, gene = line.split("\t")
  97.             duiying[lable.strip("\n")] = gene.strip("\n")
  98.     return duiying
  101. def motif_Z(trap, site):
  102.     if trap[10:30][site: site + 3] in ["CAG", "CGA", "CAA"]:
  103.         return 1
  104.     else:
  105.         return 0
  108. def motif_F(trap, site):
  109.     if "CCA" in trap[10:30][site-2: site + 3]:
  110.         return 1
  111.     else:
  112.         return 0
  115. lable_map_gene = read_duiying("../duiying.txt")
  116. all_cds = read_cds("../mart_export.txt")
  117. all_ref = read_ref("ref.fa")
  118. all_trapseq = read_seq("new508.filter.fa")
  119. with open("pass_lst_f.txt", 'a') as wocao:
  120.     wocao.write(
  121.         "TRAP ID\tGene\tgRNA sequence\tstrand\tC>T position\tC>T efficiency\n")
  122. with open("pass_lst_z.txt", 'a') as wocao:
  123.     wocao.write(
  124.         "TRAP ID\tGene\tgRNA sequence\tstrand\tC>T position\tC>T efficiency\n")
  125. for name, traps in all_ref.items():
  126.     gene = lable_map_gene[name]
  127.     cds_lst = [vals for keys, vals in all_cds.items() if gene in keys]
  128.     pass_lst = []
  129.     for m in re.finditer("C", traps[10:30]):
  130.         if motif_Z(traps, m.start()):
  131.             if stop_test_for_Z((m.start() + 1), cds_lst, traps) and name in list(all_trapseq.keys()):
  132.                 i = 0
  133.                 for gg in all_trapseq[name]:
  134.                     if gg[10:30][m.start()] == "T":
  135.                         i = i + 1
  136.                 eff = i/len(all_trapseq[name])
  137.                 pass_lst.append(
  138.                     [name, traps[10:33], "+", (m.start() + 1), eff])
  139.     if len(pass_lst) > 0:
  140.         with open("pass_lst_z.txt", 'a') as wocao:
  141.             for nima in pass_lst:
  142.                 [trap_id, trap_seq, strand, ctposition, efftion] = nima
  143.                 wocao.write("{}\t{}\t{}\t{}\t{}\t{}\n".format(
  144.                     trap_id, gene, trap_seq, strand, ctposition, efftion))
  146. for name, traps in all_ref.items():
  147.     gene = lable_map_gene[name]
  148.     cds_lst = [vals for keys, vals in all_cds.items() if gene in keys]
  149.     pass_lst = []
  150.     for m in re.finditer("C", traps[10:30]):
  151.         if motif_F(traps, m.start()):
  152.             if stop_test_for_F(((traps[10:30][m.start()-2: m.start() + 3]).find("CCA")+1+2), cds_lst, traps) and name in list(all_trapseq.keys()):
  153.                 j = 0
  154.                 for gg in all_trapseq[name]:
  155.                     if gg[10:30][m.start()] == "T":
  156.                         j = j + 1
  157.                 eff = j/len(all_trapseq[name])
  158.                 pass_lst.append(
  159.                     [name, traps[10:33], "-", (m.start() + 1), eff])
  160.     if len(pass_lst) > 0:
  161.         with open("pass_lst_f.txt", 'a') as wocao:
  162.             for nima in pass_lst:
  163.                 [trap_id, trap_seq, strand, ctposition, efftion] = nima
  164.                 wocao.write("{}\t{}\t{}\t{}\t{}\t{}\n".format(
  165.                     trap_id, gene, trap_seq, strand, ctposition, efftion))

该博客仅记录工作流程，以作为备忘录，请勿过分借鉴，也不要转载，感谢！