OSCR

Integrative analysis of drug-gene signatures in human pluripotent stem cells reveals prazosin as a novel SQSTM1 regulator for ALS therapeutics.

Code ↔ Paper

2 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 2 matches · 1 of them tie a paragraph to a whole file, not to given lines: a weak match, whose lines are not tinted
  1. [1] § STAR★Methods › Method details › Generation of heterozygous and homozygous SQSTM1 knockout PSC lines ↔ crispor.py, lines 2721–2779 · score 0.70 · spCas9, length polymorphism, mutation induced, enzyme, RFLP, CRISPOR
  2. [2] § STAR★Methods › Method details › Prazosin treatment of zebrafish with sqstm1 knockdown ↔ js/jquery-ui.min.js, the whole file · a weak match · score 0.52 · duration, Touch, escape, submitted, LI, fast

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 4,761 lines · 195 KB · other · 1 match

  1. #!/data/www/crispor/venv/bin/python3
  2. # if you do not want the hardcoded PATH above, delete this line and the one above to use the default Python3 interpreter
  3. #!/usr/bin/env python3
  4. # I know that this line looks unprofessional to you, but modifying the PATH on a shared Apache webserver is not obvious.
  5. # the tefor crispr tool
  6. # can be run as a CGI or from the command line
  7. # python std library
  8. import subprocess, tempfile, optparse, logging, atexit, glob, shutil, signal, pdb
  9. import http.cookies, time, sys, cgi, re, random, platform, os, pipes
  10. import hashlib, base64, string, logging, operator, urllib.request, urllib.parse, urllib.error, time
  11. import traceback, json, pwd, gzip, zlib
  12. from io import StringIO
  13. from collections import defaultdict, namedtuple
  14. from datetime import datetime
  15. from itertools import product
  16. from os.path import abspath, basename, dirname, isdir, isfile, join, relpath
  17. try:
  18. # prefer the pip package, it's more up-to-date than the native package
  19. import pysqlite3 as sqlite3
  20. SQLITEERROR=pysqlite3.dbapi2.OperationalError
  21. except:
  22. import sqlite3
  23. SQLITEERROR=sqlite3.OperationalError
  24. try:
  25. from collections import OrderedDict
  26. except ImportError:
  27. from ordereddict import OrderedDict # python2.6 users: run 'sudo pip install ordereddict'
  28. # for matplotlib, improves "import" performance
  29. os.environ["MPLCONFIGDIR"] = "/tmp/matplotlib-cache"
  30. # try to load external dependencies
  31. # we're going into great lengths to create a readable error message
  32. needModules = set(["pytabix", "twobitreader", "pandas", "matplotlib", "scipy"])
  33. try:
  34. import tabix # if not found, install with 'pip install pytabix'
  35. needModules.remove("pytabix")
  36. except:
  37. pass
  38. try:
  39. import twobitreader # if not found, install with 'pip install twobitreader'
  40. needModules.remove("twobitreader")
  41. except:
  42. pass
  43. try:
  44. import pandas # required by doench2016 score. install with 'pip install pandas'
  45. needModules.remove("pandas")
  46. import scipy # required by doench2016 score. install with 'pip install scipy'
  47. needModules.remove("scipy")
  48. import matplotlib # required by doench2016 score. install with 'pip install matplotlib'
  49. needModules.remove("matplotlib")
  50. import numpy # required by doench2016 score. install with 'pip install numpy'
  51. needModules.remove("numpy")
  52. except:
  53. pass
  54. if len(needModules)!=0:
  55. print("Content-type: text/html\n")
  56. print(("Python interpreter path: %s<p>" % sys.executable))
  57. print(("These python modules were not found: %s<p>" % ",".join(needModules)))
  58. print("To install all requirements in one line, run: sudo pip install biopython numpy scikit-learn==0.16.1 pandas twobitreader<p>")
  59. sys.exit(0)
  60. # our own eff scoring library
  61. import crisporEffScores
  62. # don't report print as an error
  63. # pylint: disable=E1601
  64. # optional module for Excel export as native .xls files
  65. # install with 'apt-get install python-xlwt' or 'pip install xlwt'
  66. xlwtLoaded = True
  67. try:
  68. import xlwt
  69. except:
  70. sys.stderr.write("crispor.py - warning - the python xlwt module is not available\n")
  71. xlwtLoaded = False
  72. # optional module for mysql support
  73. #try:
  74. #import MySQLdb
  75. #mysqldbLoaded = True
  76. #except:
  77. #mysqldbLoaded = False
  78. # version of crispor
  79. versionStr = "5.2"
  80. # contact email
  81. contactEmail='[email hidden]'
  82. # url to this server
  83. ctBaseUrl = "http://crispor-max.tefor.net/temp/customTracks"
  84. # write debug output to stdout
  85. DEBUG = False
  86. #DEBUG = True
  87. # use bowtie for off-target search?
  88. useBowtie = False
  89. # calculate the efficienc scores?
  90. doEffScoring = True
  91. # system-wide temporary directory
  92. #TEMPDIR = os.environ.get("TMPDIR", "/var/tmp")
  93. TEMPDIR = "/var/tmp"
  94. # a hack for cluster jobs at UCSC:
  95. # - default to ramdisk
  96. if isdir("/scratch/tmp"):
  97. TEMPDIR = "/dev/shm/"
  98. # skipAlign is useful if your input sequence is not in the genome at all
  99. # - don't do bwasw
  100. # - this will trigger auto-ontarget: any perfect match is the on-target
  101. # - do not calculate efficiency scores
  102. skipAlign = False
  103. # prefix in html statements before the directories "image/", "style/" and "js/"
  104. HTMLPREFIX = ""
  105. # alternative directory on local disk where image/, style/ and js/ are located
  106. HTMLDIR = "/usr/local/apache/htdocs/crispor/"
  107. # directory of crispor.py
  108. baseDir = dirname(__file__)
  109. # filename of this script, usually crispor.py
  110. myName = basename(__file__)
  111. # the segments.bed files use abbreviated genomic region names
  112. segTypeConv = {"ex":"exon", "in":"intron", "ig":"intergenic"}
  113. # directory for processed batches of offtargets ("cache" of bwa results)
  114. batchDir = join(baseDir,"temp")
  115. # sqlite3 db with gzipped old json batch files, to avoid hitting the ext4 inode limits
  116. batchArchive = "/data/crisporJobArchive.db"
  117. # the file where the sqlite job queue is stored
  118. #JOBQUEUEDB = join(TEMPDIR, "crisporJobs.db") # TEMPDIR is mapped away for security reasons under Redhat/Centos for CGIs
  119. JOBQUEUEDB = "/data/www/temp/crisporJobs.db"
  120. # alternatively: connection info for mysql
  121. jobQueueMysqlConn = {"socket":None, "host":None, "user": None, "password" : None}
  122. # directory for platform-independent scripts (e.g. Heng Li's perl SAM parser)
  123. scriptDir = join(baseDir, "scripts")
  124. # directory for helper binaries (e.g. BWA)
  125. # system() is one of 'Linux', 'Darwin', 'Windows', machine() is one of 'x86_64', 'arm64', 'aarch64'
  126. os_name = platform.system() # 'Linux', 'Darwin', 'Windows'
  127. arch = platform.machine() # 'x86_64', 'arm64', 'aarch64'
  128. binDir = abspath(join(baseDir, "bin", platform.system()+"-"+platform.machine()))
  129. # directory for genomes
  130. genomesDir = join(baseDir, "genomes")
  131. DEFAULTORG = 'hg19'
  132. DEFAULTSEQ = 'cttcctttgtccccaatctgggcgcgcgccggcgccccctggcggcctaaggactcggcgcgccggaagtggccagggcgggggcgacctcggctcacagcgcgcccggctattctcgcagctcaccatgGATGATGATATCGCCGCGCTCGTCGTCGACAACGGCTCCGGCATGTGCAAGGCCGGCTTCGCGGGCGACGATGCCCCCCGGGCCGTCTTCCCCTCCATCGTGGGGCGCC'
  133. # used if hg19 is not available
  134. ALTORG = 'sacCer3'
  135. ALTSEQ = 'ATTCTACTTTTCAACAATAATACATAAACatattggcttgtggtagCAACACTATCATGGTATCACTAACGTAAAAGTTCCTCAATATTGCAATTTGCTTGAACGGATGCTATTTCAGAATATTTCGTACTTACACAGGCCATACATTAGAATAATATGTCACATCACTGTCGTAACACTCT'
  136. pamDesc = [ ('NGG','20bp-NGG - Sp Cas9, SpCas9-HF1, eSpCas9 1.1'),
  137. ('NNG','20bp-NNG - Cas9 S. canis'),
  138. ('NGN','20bp-NGN - SpG'),
  139. ('NNGT','20bp-NNGT - Cas9 S. canis - high efficiency PAM, recommended'),
  140. ('NAA','20bp-NAA - iSpyMacCas9'),
  141. ('TTN', 'TTN-23bp - hfCas12Max - as recommended by Synthego'), # Casey Jowdy by email
  142. ('TNN','TNN-23bp - hfCas12Max, broader PAM, as recommended by Synthego'), #
  143. ('NGG-22', 'NGG-22bp - eSpOT-ON (ePsCas9), as recommended by Synthego'),
  144. ('NNGRRT','21bp-NNG(A/G)(A/G)T - Cas9 S. Aureus'),
  145. ('NNGRRT-20','20bp-NNG(A/G)(A/G)T - Cas9 S. Aureus with 20bp-guides'),
  146. ('NGK','20bp-NG(G/T) - xCas9, recommended PAM, see notes'),
  147. #('NGN','20bp-NGN or GA(A/T) - xCas9 (low efficiency, not recommended)'),
  148. #('NGG-BE1','20bp-NGG - BaseEditor1, modifies C->T'),
  149. ('NNNRRT','21bp-NNN(A/G)(A/G)T - KKH SaCas9'),
  150. ('NNNRRT-20','20bp-NNN(A/G)(A/G)T - KKH SaCas9 with 20bp-guides'),
  151. ('NGA','20bp-NGA - Cas9 S. Pyogenes mutant VQR'),
  152. ('NNNNCC','24bp-NNNNCC - Nme2Cas9'),
  153. ('NGCG','20bp-NGCG - Cas9 S. Pyogenes mutant VRER'),
  154. ('NNAGAA','20bp-NNAGAA - Cas9 S. Thermophilus'),
  155. ('NGGNG','20bp-NGGNG - Cas9 S. Thermophilus'),
  156. ('NNNNGMTT','20bp-NNNNG(A/C)TT - Cas9 N. Meningitidis'),
  157. ('NNNNACA','20bp-NNNNACA - Cas9 Campylobacter jejuni, original PAM'),
  158. ('NNNNRYAC','22bp-NNNNRYAC - Cas9 Campylobacter jejuni, revised PAM'),
  159. ('NNNVRYAC','22bp-NNNVRYAC - Cas9 Campylobacter jejuni, opt. efficiency'),
  160. ('TTCN','TTCN-20bp - CasX'),
  161. ('TTTV','TTT(A/C/G)-23bp - Cas12a (Cpf1) - recommended, 23bp guides'),
  162. ('TTTV-21','TTT(A/C/G)-21bp - Cas12a (Cpf1) - 21bp guides recommended by IDT'),
  163. ('TTTN','TTTN-23bp - Cas12a (Cpf1) - low efficiency'),
  164. ('ATTN','ATTN-23bp - BhCas12b v4'),
  165. ('NGTN','NGTN-23bp - ShCAST/AcCAST, Strecker et al, Science 2019'),
  166. ('TYCV','T(C/T)C(A/C/G)-23bp - TYCV As-Cpf1 K607R'),
  167. ('TATV','TAT(A/C/G)-23bp - TATV As-Cpf1 K548V'),
  168. ('TTTA','TTTA-23bp - TTTA LbCpf1'),
  169. ('TCTA','TCTA-23bp - TCTA LbCpf1'),
  170. ('TCCA','TCCA-23bp - TCCA LbCpf1'),
  171. ('CCCA','CCCA-23bp - CCCA LbCpf1'),
  172. ('GGTT','GGTT-23bp - CCCA LbCpf1'),
  173. ('YTTV','YTTV-20bp - MAD7 Nuclease, Lui, Schiel, Maksimova et al, CRISPR J 2020'),
  174. ('TTYN','TTYN- or VTTV- or TRTV-23bp - enCas12a E174R/S542R/K548R - Kleinstiver et al Nat Biot 2019'),
  175. ('NNNNCNAA','20bp-NNNNCNAA - Thermo Cas9 - Walker et al, Metab Eng Comm 2020'),
  176. ('NNN','20bp-NNN - SpRY, Walton et al Science 2020'), # https://science.sciencemag.org/content/368/6488/290.abstract
  177. ('NRN','20bp-NRN - SpRY (high efficiency PAM)'),
  178. ('NYN','20bp-NYN - SpRY (low efficiency PAM)'),
  179. #('VTTV','(A/C)TT(A/C)-23bp - enCas12a S542R - Kleinstiver et al Nat Biot 2019'),
  180. #('TRTV','T(A/G)T(A/C)-23bp - enCas12a K548R - Kleinstiver et al Nat Biot 2019'),
  181. ]
  182. DEFAULTPAM = 'NGG'
  183. # the default base editor modification window
  184. DEFAULTBEWIN = "1-7"
  185. # for some PAMs, there are alternative main PAMs. These are also shown on the main sequence panel
  186. multiPams = {
  187. #"NGN" : ["GAW"],
  188. "TTYN" : ["VTTV", "TRTV"]
  189. }
  190. # these PAMs are not specific. Allow only short sequences for them.
  191. slowPams = ["TTYN", "NNG"]
  192. # allow only very short sequences for these
  193. verySlowPams = ["NNN", "NRN", "NYN"]
  194. # for some PAMs, we allow other alternative motifs when searching for offtargets
  195. # MIT and eCrisp do that, they use the motif NGG + NAG, we add one more, based on the
  196. # on the guideSeq results in Tsai et al, Nat Biot 2014
  197. # The NGA -> NGG rule was described by Kleinstiver...Young 2015 "Improved Cas9 Specificity..."
  198. # NNGTRRT rule for S. aureus is in the new protocol "SaCas9 User manual"
  199. # ! the length of the alternate PAM has to be the same as the original PAM!
  200. offtargetPams = {
  201. "NGG" : ["NAG","NGA"],
  202. #"NGN" : ["GAW"],
  203. "NGK" : ["GAW"],
  204. "NGA" : ["NGG"],
  205. "NNGRRT" : ["NNGRRN"],
  206. "TTTV" : ["TTTN"],
  207. 'ATTN' : ["TTTN", "GTTN"],
  208. "TTYN" : ["VTTV", "TRTV"]
  209. }
  210. # maximum size of an input sequence
  211. MAXSEQLEN = 2300
  212. # maximum input size when specifying "no genome"
  213. MAXSEQLEN_NOGENOME = 25000
  214. # maximum input size when using xCas9 or sCanis
  215. MAXSEQLEN2 = 600
  216. # maximum input size for NNN SpRY or similar PAMs
  217. MAXSEQLEN3 = 150
  218. # BWA: allow up to X mismatches
  219. maxMMs=4
  220. # maximum number of occurences in the genome to get flagged as repeats.
  221. # This is used in bwa samse, when converting the same file
  222. # and for warnings in the table output.
  223. MAXOCC = 60000
  224. # the BWA queue size is 2M by default. We derive the queue size from MAXOCC
  225. MFAC = 2000000/MAXOCC
  226. # the length of the guide sequence, set by setupPamInfo
  227. GUIDELEN=None
  228. # length of the PAM sequence
  229. PAMLEN=None
  230. # the name of the base editor, if any. This is the flag to activate
  231. # baseEditor mode in the UI
  232. baseEditor = None
  233. # input sequences are extended by X basepairs so we can calculate the efficiency scores
  234. # and can better design primers
  235. FLANKLEN=100
  236. # the name of the currently processed batch, assigned only once
  237. # in readBatchParams and only for json-type batches
  238. batchName = ""
  239. # are we doing a Cpf1 run?
  240. # this variable changes almost all processing and
  241. # has to be set on program start, as soon as we know
  242. # the PAM we're running on
  243. pamIsFirst=None
  244. saCas9Mode=False
  245. # Highly-sensitive mode (not for CLI mode):
  246. # MAXOCC is increased in processSubmission() and in the html UI if only one
  247. # guide seq is run
  248. # Also, the number of allowed mismatches is increased to 5 instead of 4
  249. #HIGH_MAXOCC=600000
  250. #HIGH_maxMMs=5
  251. # minimum off-target score of standard off-targets (those that end with NGG)
  252. # This should probably be based on the CFD score these days
  253. # But for now, I'll let the user do the filtering
  254. MINSCORE = 0.0
  255. # minimum off-target score for alternative PAM off-targets
  256. # There is not a lot of data to support this cutoff, but it seems
  257. # reasonable to have at least some cutoff, as otherwise we would show
  258. # NAG and NGA like NGG and the data shows clearly that the alternative
  259. # PAMs are not recognized as well as the main NGG PAM.
  260. # so for now, I just filter out very degenerative ones. the best solution
  261. # would be to have a special penalty on the CFD score, but CFS does not
  262. # support non-NGG PAMs (is this actually true?)
  263. ALTPAMMINSCORE = 1.0
  264. # how much shall we extend the guide after the PAM to match restriction enzymes?
  265. pamPlusLen = 5
  266. # global flag to indicate if we're run from command line or as a CGI
  267. commandLineMode = False
  268. # names/order of efficiency scores to show in UI
  269. cas9ScoreNames = ["fusi", "crisprScan", "rs3"]
  270. allScoreNames = ["fusi", "chariRank", "ssc", "wuCrispr", "doench", "wang", "crisprScan", "ccTop", "rs3"]
  271. mutScoreNames = []
  272. spCas9MutScoreNames = ["oof", 'lindel'] # lindel is only added for spCas9
  273. otherMutScoreNames = ["oof"] # lindel is only added for spCas9
  274. cpf1ScoreNames = ["seqDeepCpf1"]
  275. saCas9ScoreNames = ["najm"]
  276. # to make the CFD more comparable to the MIT score, Nicholas Parkinson suggests to multiply it with 100.
  277. # can be switched on with the URL argument fixCfd=1
  278. doCfdFix=False
  279. # how many digits shall we show for each score? default is 0
  280. scoreDigits = {
  281. "ssc" : 1,
  282. }
  283. # List of AddGene plasmids, their long and short names:
  284. addGenePlasmids = [
  285. ("43860", ("MLM3636 (Joung lab)", "MLM3636")),
  286. ("49330", ("pAc-sgRNA-Cas9 (Liu lab)", "pAcsgRnaCas9")),
  287. ("42230", ("pX330-U6-Chimeric_BB-CBh-hSpCas9 (Zhang lab) + derivatives", "pX330")),
  288. ("52961", ("lentiCRISPR v2 (Zhang lab)", "lentiCrispr")),
  289. ("52963", ("lentiGuide-Puro (Zhang lab)", "lentiGuide-Puro")),
  290. ]
  291. addGenePlasmidsAureus = [
  292. ("61591", ("pX601-AAV-CMV::NLS-SaCas9-NLS-3xHA-bGHpA;U6::BsaI-sgRNA (Zhang lab)", "pX601")),
  293. ("61592", ("pX600-AAV-CMV::NLS-SaCas9-NLS-3xHA-bGHpA (Zhang lab)", "pX600")),
  294. ("61593", ("pX602-AAV-TBG::NLS-SaCas9-NLS-HA-OLLAS-bGHpA;U6::BsaI-sgRNA (Zhang lab)", "pX602")),
  295. ("65779", ("VVT1 (Joung lab)", "VVT1"))
  296. ]
  297. # list of AddGene primer 5' and 3' extensions, one for each AddGene plasmid
  298. # format: prefixFw, prefixRw, u6-G-suffix, restriction enzyme, link to protocol
  299. addGenePlasmidInfo = {
  300. "43860" : ("ACACC", "AAAAC", "G", "BsmBI", "https://www.addgene.org/static/data/plasmids/43/43860/43860-attachment_T35tt6ebKxov.pdf"),
  301. "49330" : ("TTC", "AAC", "", "Bsp QI", "http://bio.biologists.org/content/3/1/42#sec-9"),
  302. "42230" : ("CACC", "AAAC", "", "Bbs1", "https://www.addgene.org/static/data/plasmids/52/52961/52961-attachment_B3xTwla0bkYD.pdf"),
  303. "52961" : ("CACC", "AAAC", "", "BsmBI", "https://www.addgene.org/static/data/plasmids/52/52961/52961-attachment_B3xTwla0bkYD.pdf"),
  304. "61591" : ("CACC", "AAAC", "", "BsaI", "https://www.addgene.org/static/data/plasmids/61/61591/61591-attachment_it03kn5x5O6E.pdf"),
  305. "61592" : ("CACC", "AAAC", "", "BsaI", "https://www.addgene.org/static/data/plasmids/61/61592/61592-attachment_iAbvIKnbqNRO.pdf"),
  306. "61593" : ("CACC", "AAAC", "", "BsaI", "https://www.addgene.org/static/data/plasmids/61/61592/61592-attachment_iAbvIKnbqNRO.pdf"),
  307. "65779": ("CACC", "AAAC", "", "BsmBI (aka Esp3l)", "https://www.addgene.org/static/data/plasmids/65/65779/65779-attachment_G8oNyvV6pA78.pdf"),
  308. "52963": ("CACC", "AAAC", "", "BsmBI (aka Esp3l)", "https://www.addgene.org/static/data/plasmids/52/52963/52963-attachment_IPB7ZL_hJcbm.pdf")
  309. }
  310. # the barcodes for subpool tagging for oligo pool tables
  311. satMutBarcodes = [
  312. (0, "No Subpool barcode"),
  313. (1, "Subpool 1: CGGGTTCCGT/GCTTAGAATAGAA"),
  314. (2, "Subpool 2: GTTTATCGGGC/ACTTACTGTACC"),
  315. (3, "Subpool 3: ACCGATGTTGAC/CTCGTAATAGC"),
  316. (4, "Subpool 4: GAGGTCTTTCATGC/CACAACATA"),
  317. (5, "Subpool 5: TATCCCGTGAAGCT/TTCGGTTAA"),
  318. (6, "Subpool 6: TAGTAGTTCAGACGC/ATGTACCC"),
  319. (7, "Subpool 7: GGATGCATGATCTAG/CATCAAGC"),
  320. (8, "Subpool 8: ATGAGGACGAATCT/CACCTAAAG"),
  321. (9, "Subpool 9: GGTAGGCACG/TAAACTTAGAACC"),
  322. (10, "Subpool 10: AGTCATGATTCAG/GTTGCAAGTCTAG"),
  323. ]
  324. # Restriction enzyme supplier codes
  325. rebaseSuppliers = {
  326. "B":"Life Technologies",
  327. "C":"Minotech",
  328. "E":"Agilent",
  329. "I":"SibEnzyme",
  330. "J":"Nippon Gene",
  331. "K":"Takara",
  332. "M":"Roche",
  333. "N":"NEB",
  334. "O":"Toyobo",
  335. "Q":"Molecular Biology Resources",
  336. "R":"Promega",
  337. "S":"Sigma",
  338. "V":"Vivantis",
  339. "X":"EURx",
  340. "Y":"SinaClon BioScience"
  341. }
  342. # labels and descriptions of eff. scores
  343. scoreDescs = {
  344. "doench" : ("Doench '14", "Range: 0-100. Linear regression model trained on 880 guides transfected into human MOLM13/NB4/TF1 cells (three genes) and mouse cells (six genes). Delivery: lentivirus. The Fusi score can be considered an updated version this score, as their training data overlaps a lot. See <a target='_blank' href='http://www.nature.com/nbt/journal/v32/n12/full/nbt.3026.html'>Doench et al.</a>"),
  345. "wuCrispr" : ("Wu-Crispr", "Range 0-100. Aka 'Wong score'. SVM model trained on previously published data. The aim is to identify only a subset of efficient guides, many guides will have a score of 0. Takes into account RNA structure. See <a target='_blank' href='https://genomebiology.biomedcentral.com/articles/10.1186/s13059-015-0784-0'>Wong et al., Gen Biol 2015</a>"),
  346. "ssc" : ("Xu", "Range ~ -2 - +2. Aka 'SSC score'. Linear regression model trained on data from &gt;1000 genes in human KBM7/HL60 cells (Wang et al) and mouse (Koike-Yusa et al.). Delivery: lentivirus. Ranges mostly -2 to +2. See <a target='_blank' href='http://genome.cshlp.org/content/early/2015/06/10/gr.191452.115'>Xu et al.</a>"),
  347. "crisprScan" : ["Moreno-Mateos", "Also called 'CrisprScan'. Range: mostly 0-100. Linear regression model, trained on data from 1000 guides on &gt;100 genes, from zebrafish 1-cell stage embryos injected with mRNA. See <a target=_blank href='http://www.nature.com/nmeth/journal/v12/n10/full/nmeth.3543.html'>Moreno-Mateos et al.</a>. Recommended for guides transcribed <i>in-vitro</i> (T7 promoter). Click to sort by this score. Note that under 'Show all scores', you can find a Doench2016 model trained on Zebrafish scores, Azimuth in-vitro, which should be slightly better than this model for zebrafish."],
  348. "wang" : ("Wang", "Range: 0-100. SVM model trained on human cell culture data on guides from &gt;1000 genes. The Xu score can be considered an updated version of this score, as the training data overlaps a lot. Delivery: lentivirus. See <a target='_blank' href='http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3972032/'>Wang et al.</a>"),
  349. "chariRank" : ("Chari", "Range: 0-100. Support Vector Machine, converted to rank-percent, trained on data from 1235 guides targeting sequences that were also transfected with a lentivirus into human 293T cells. See <a target='_blank' href='http://www.nature.com/nmeth/journal/v12/n9/abs/nmeth.3473.html'>Chari et al.</a>"),
  350. "fusi" : ("Doench '16", "Aka the 'Fusi-Score', since V4.4 using the version 'Azimuth', scores are slightly different than before April 2018 but very similar (click 'show all' to see the old scores). Range: 0-100. Boosted Regression Tree model, trained on data produced by Doench et al (881 guides, MOLM13/NB4/TF1 cells + unpublished additional data). Delivery: lentivirus. See <a target='_blank' href='http://biorxiv.org/content/early/2015/06/26/021568'>Fusi et al. 2015</a> and <a target='_blank' href='http://www.nature.com/nbt/journal/v34/n2/full/nbt.3437.html'>Doench et al. 2016</a> and <a target=_blank href='https://crispr.ml/'>crispr.ml</a>. Recommended for guides expressed in cells (U6 promoter). Click to sort the table by this score."),
  351. "fusiOld" : ("OldDoench '16", "The original implementation of the Doench 2016 score, as received from John Doench. The scores are similar, but not exactly identical to the 'Azimuth' version of the Doench 2016 model that is currently the default on this site, since Apr 2018."),
  352. "rs3" : ("Doench-RuleSet3", "The Doench Rule Set 3 (RS3) score (-200-+200). Similar to the Doench 2014 and Doench 2016/Fusi/Azimuth score, but updated and more accurate. See <a href='https://www.nature.com/articles/s41467-022-33024-2' target=_blank>. Scores shown are multiplied with 100 for easier display. RS3 is configured here to use the Hsu-TRACR sequence."),
  353. "najm" : ("Najm 2018", "A modified version of the Doench 2016 score ('Azimuth'), by Mudra Hegde for S. aureus Cas9. Range 0-100. See <a target=_blank href='https://www.nature.com/articles/nbt.4048'>Najm et al 2018</a>."),
  354. "ccTop" : ("CCTop", "The efficiency score used by CCTop, called 'crisprRank'."),
  355. "aziInVitro" : ("Azimuth in-vitro", "The Doench 2016 model trained on the Moreno-Mateos zebrafish data. Unpublished model, gratefully provided by J. Listgarden. This should be better than Moreno-Mateos, but we have not found the time to evaluate it yet."),
  356. "housden" : ("Housden", "Range: ~ 1-10. Weight matrix model trained on data from Drosophila mRNA injections. See <a target='_blank' href='http://stke.sciencemag.org/content/8/393/rs9.long'>Housden et al.</a>"),
  357. "proxGc" : ("ProxGCCount", "Number of GCs in the last 4pb before the PAM"),
  358. "seqDeepCpf1" : ("DeepCpf1", "Range: ~ 0-100. Convolutional Neural Network trained on ~20k Cpf1 lentiviral guide results. This is the score without DNAse information, 'Seq-DeepCpf1' in the paper. See <a target='_blank' href='https://www.nature.com/articles/nbt.4061'>Kim et al. 2018</a>"),
  359. "oof" : ("Out-of-Frame", "Range: 0-100. Out-of-Frame score, only for deletions. Predicts the percentage of clones that will carry out-of-frame deletions, based on the micro-homology in the sequence flanking the target site. See <a target='_blank' href='http://www.nature.com/nmeth/journal/v11/n7/full/nmeth.3015.html'>Bae et al. 2014</a>. Click the score to show the predicted deletions."),
  360. "lindel": ("Lindel", "Wei Chen Frameshift ratio (0-100). Predicts probability of a frameshift caused by any type of insertion or deletion. See <a href='https://academic.oup.com/nar/article/47/15/7989/5511473'>Wei Chen et al, Bioinf 2018</a>. Click the score to see the most likely deletions and insertions."),
  361. }
  362. # the headers for the guide and offtarget output files
  363. guideHeaders = ["guideId", "targetSeq", "mitSpecScore", "cfdSpecScore", "offtargetCount", "targetGenomeGeneLocus"]
  364. offtargetHeaders = ["guideId", "guideSeq", "offtargetSeq", "mismatchPos", "mismatchCount", "mitOfftargetScore", "cfdOfftargetScore", "chrom", "start", "end", "strand", "locusDesc"]
  365. # library descriptions
  366. libLabels = [
  367. # https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4486245/
  368. ("human_brunello" , "Human, Brunello, Doench Nat Bio 2016 (recommended)"),
  369. ("human_avana" , "Human, Avana, Doench Nat Bio 2016"),
  370. ("human_geckov2" , "Human, GeCKO V2, Sanjana Nat Meth 2014"),
  371. ("mouse_brie" , "Mouse, Brie, Doench Nat Bio 2016 (recommended)"),
  372. ("mouse_geckov2" , "Mouse, GeCKO V2, Sanjana Nat Meth 2014"),
  373. ("mouse_asiago" , "Mouse, Asiago, Doench Nat Bio 2016"),
  374. ]
  375. # a file crispor.conf in the directory of the script allows to override any global variable
  376. myDir = dirname(__file__)
  377. confPath =join(myDir, "crispor.conf")
  378. if isfile(confPath):
  379. exec(open(confPath).read())
  380. #execfile(confPath)
  381. cgiParams = None
  382. # ====== END GLOBALS ============
  383. def setupPamInfo(pam):
  384. " modify a few globals based on the current pam "
  385. global GUIDELEN
  386. global pamIsFirst
  387. global addGenePlasmids
  388. global PAMLEN
  389. global scoreNames
  390. global baseEditor
  391. global saCas9Mode
  392. global mutScoreNames
  393. global isSpg
  394. PAMLEN = len(pam)
  395. pamIsFirst = False
  396. scoreNames = cas9ScoreNames
  397. pamOpt = None
  398. if "-" in pam:
  399. pam, pamOpt = pam.split("-")
  400. if pamOpt=="BE1":
  401. baseEditor = "BE1"
  402. elif pamOpt=="spg":
  403. isSpg = True
  404. if pamIsCasX(pam):
  405. logging.debug("switching on CasX mode, guide length is 20bp")
  406. GUIDELEN = 20
  407. pamIsFirst = True
  408. scoreNames = cpf1ScoreNames
  409. if pamIsCas12max(pam):
  410. logging.debug("switching on hfCas12max mode, guide length is 20bp")
  411. GUIDELEN = 20
  412. pamIsFirst = True
  413. scoreNames = cpf1ScoreNames
  414. elif pamIsCpf1(pam):
  415. logging.debug("switching on Cpf1 mode, guide length is 23bp")
  416. GUIDELEN = 23
  417. pamIsFirst = True
  418. scoreNames = cpf1ScoreNames
  419. #if pamOpt:
  420. #GUIDELEN=int(pamOpt)
  421. elif pam=="NGTN":
  422. logging.debug("switching on Cpf1 mode for ShCAST, guide length is 23bp")
  423. GUIDELEN = 23
  424. pamIsFirst = True
  425. elif pam=="NNNNRYAC" or pam=="NNNVRYAC":
  426. GUIDELEN = 22
  427. elif pam=="NNGRRT" or pam=="NNNRRT":
  428. logging.debug("switching on S. aureus mode, guide length is 21bp")
  429. addGenePlasmids = addGenePlasmidsAureus
  430. GUIDELEN = 21
  431. #if pamOpt=="20":
  432. #GUIDELEN=20
  433. saCas9Mode = True
  434. scoreNames = saCas9ScoreNames
  435. elif pam=="NNNNCC":
  436. GUIDELEN = 24
  437. else:
  438. GUIDELEN = 20
  439. if pamOpt and pamOpt.isnumeric():
  440. GUIDELEN = int(pamOpt)
  441. if (GUIDELEN==20 or GUIDELEN==22) and pam=="NGG":
  442. mutScoreNames = spCas9MutScoreNames
  443. else:
  444. mutScoreNames = otherMutScoreNames
  445. logging.debug("Enzyme info: pam=%s, guideLen=%d, pamIsFirst=%s, saCas9Mode=%s" %
  446. (pam, GUIDELEN, pamIsFirst, saCas9Mode))
  447. return pam
  448. # ==== CLASSES =====
  449. class JobQueue:
  450. """
  451. simple job queue, using a db table as a backend
  452. jobs have different types and status. status can be updated while they run
  453. job running times are kept and old job info is kept in a separate table
  454. >>> q = JobQueue()
  455. >>> q.openSqlite()
  456. >>> q.clearJobs()
  457. >>> q.waitCount()
  458. 0
  459. >>> q.addJob("search", "abc123", "myParams")
  460. True
  461. only one job per jobId
  462. >>> q.addJob("search", "abc123", "myParams")
  463. False
  464. >>> q.waitCount()
  465. 1
  466. >>> q.getStatus("abc123")
  467. 'Waiting'
  468. >>> q.startStep("abc123", "bwa", "Alignment with BWA")
  469. >>> q.getStatus("abc123")
  470. 'Alignment with BWA'
  471. >>> jobType, jobId, paramStr = q.popJob()
  472. >>> q.waitCount()
  473. 0
  474. >>> q.jobDone("abc123")
  475. >>> q.waitCount()
  476. 0
  477. can't pop from an empty queue
  478. #>>> q.popJob()
  479. #(None, None, None)
  480. #>>> os.system("rm /tmp/tempCrisporTest.db")
  481. #0
  482. """
  483. _queueDef = (
  484. 'CREATE TABLE IF NOT EXISTS %s '
  485. '('
  486. ' jobType text,' # either "index" or "search"
  487. ' jobId text %s,' # unique identifier
  488. ' paramStr text,' # parameters for jobs, like db, options, etc.
  489. ' isRunning int DEFAULT 0,' # indicates steps have started, done jobs are moved to doneJobs table
  490. ' stepName text,' # currently step, internal step name for timings
  491. ' stepLabel text,' # current step, human-readable status of job, for UI
  492. ' lastUpdate float,' # time of last update
  493. ' stepTimes text,' # comma-sep list of whole msecs, one per step
  494. ' startTime text ' # date+time when job was put into queue
  495. ')')
  496. #def __init__(self):
  497. #" no inheritance needed here "
  498. #self.openSqlite(JOBQUEUEDB)
  499. def openSqlite(self, dbName=JOBQUEUEDB):
  500. self.dbName = dbName
  501. # isolation_level=None = autocommit mode: we manage transactions explicitly
  502. # timeout=30: wait up to 30s for locks (multiple daemons + CGI share this DB)
  503. self.conn = sqlite3.connect(dbName, timeout=30, isolation_level=None)
  504. # WAL mode allows concurrent readers + writer, essential for CGI + daemon(s)
  505. result = self.conn.execute("PRAGMA journal_mode=WAL;").fetchone()
  506. if result[0] != "wal":
  507. logging.warn("Could not set WAL mode, journal_mode is: %s (is another connection open with the old mode?)" % result[0])
  508. #self.conn.set_trace_callback(print) # for debugging: print all sql statements
  509. self._chmodJobDb()
  510. try:
  511. self.conn.execute(self._queueDef % ("queue", "PRIMARY KEY"))
  512. except SQLITEERROR as ex:
  513. errAbort("cannot open the sqlite jobs file %s: %s" % (JOBQUEUEDB, ex))
  514. def _chmodJobDb(self):
  515. # umask is not respected by sqlite, bug http://www.mail-archive.com/[email hidden]/msg59080.html
  516. try:
  517. os.chmod(JOBQUEUEDB, 0o666)
  518. except OSError:
  519. # if the file was created by other job, we can't chmod, as we're the CGI. Just silently ignore this
  520. pass
  521. def addJob(self, jobType, jobId, paramStr):
  522. " create a new job, returns False if not successful "
  523. self._chmodJobDb()
  524. sql = 'INSERT INTO queue (jobType, jobId, isRunning, lastUpdate, ' \
  525. 'stepTimes, paramStr, stepName, stepLabel, startTime) VALUES (:jobType, :jobId, :isRunning, :lastUpdate, ' \
  526. ':stepTimes, :paramStr, :stepName, :stepLabel, :startTime)'
  527. now = "%.3f" % time.time()
  528. values = {'jobType' : jobType, 'jobId' : jobId, 'isRunning' : 0, 'lastUpdate' : now, 'stepTimes':"", 'paramStr':paramStr, 'stepName':"wait",
  529. "stepLabel":"Waiting", "startTime":now}
  530. try:
  531. # in autocommit mode, this INSERT commits immediately
  532. self.conn.execute(sql, values)
  533. return True
  534. except sqlite3.IntegrityError:
  535. # job already in queue (e.g. resubmit, or daemon restart) - that's fine
  536. return True
  537. except SQLITEERROR:
  538. errAbort("Cannot open DB file %s. Please contact %s" % (self.dbName, contactEmail))
  539. def getStatus(self, jobId):
  540. " return current job status label or None if job is not in queue"
  541. sql = 'SELECT stepLabel FROM queue WHERE jobId=?'
  542. try:
  543. rows = self.conn.execute(sql, (jobId,)).fetchmany(1)
  544. if len(rows) == 0:
  545. status = None
  546. else:
  547. status = rows[0][0]
  548. except (StopIteration, IndexError):
  549. logging.debug("getStatus: job %s not found" % jobId)
  550. status = None
  551. return status
  552. def dump(self):
  553. " for debugging, write the whole queue table to stdout "
  554. sql = 'SELECT * FROM queue'
  555. for row in self.conn.execute(sql):
  556. print("\t".join([str(x) for x in row]))
  557. def jobInfo(self, jobId, isDone=False):
  558. " for debugging, return all job info as a tuple "
  559. print("job info<br>")
  560. if isDone:
  561. sql = 'SELECT * FROM doneJobs WHERE jobId=?'
  562. else:
  563. sql = 'SELECT * FROM queue WHERE jobId=?'
  564. try:
  565. row = next(self.conn.execute(sql, (jobId,)))
  566. except StopIteration:
  567. return []
  568. return row
  569. def startStep(self, jobId, newName, newLabel):
  570. " start a new step. Update lastUpdate, status and stepTime "
  571. try:
  572. self.conn.execute('BEGIN IMMEDIATE')
  573. sql = 'SELECT lastUpdate, stepTimes, stepName FROM queue WHERE jobId=?'
  574. logging.debug(sql)
  575. rows = self.conn.execute(sql, (jobId,)).fetchmany(1)
  576. if len(rows)==0:
  577. logging.error("startStep: no row for jobId %s" % jobId)
  578. self.conn.commit()
  579. return
  580. lastTime, timeStr, lastStep = rows[0]
  581. lastTime = float(lastTime)
  582. # append a string in format "stepName:milliSecs" to the timeStr
  583. now = time.time()
  584. timeDiff = "%d" % int((1000.0*(now - lastTime)))
  585. newTimeStr = timeStr+"%s=%s" % (lastStep, timeDiff)+","
  586. sql = 'UPDATE queue SET lastUpdate=?, stepName=?, stepLabel=?, stepTimes=?, isRunning=? WHERE jobId=?'
  587. self.conn.execute(sql, (now, newName, newLabel, newTimeStr, 1, jobId))
  588. self.conn.commit()
  589. except:
  590. self.conn.rollback()
  591. raise
  592. def jobDone(self, jobId):
  593. " remove the job from the queue and add it to the queue log"
  594. print("job done<br>")
  595. try:
  596. self.conn.execute('BEGIN IMMEDIATE')
  597. sql = 'SELECT * FROM queue WHERE jobId=?'
  598. row = self.conn.execute(sql, (jobId,)).fetchone()
  599. if row is None:
  600. logging.warn("jobDone - job %s has been removed already" % jobId)
  601. self.conn.commit()
  602. return
  603. sql = 'DELETE FROM queue WHERE jobId=?'
  604. self.conn.execute(sql, (jobId,))
  605. self.conn.commit()
  606. except:
  607. self.conn.rollback()
  608. raise
  609. # good to have a log file of the old jobs
  610. with open("doneJobs.tsv", "a") as ofh: # if this triggers an error: run 'touch doneJobs.tsv && chmod a+rw doneJobs.tsv' in the crispor dir.
  611. row = [str(x) for x in row]
  612. line = "\t".join(row)
  613. ofh.write(line)
  614. ofh.write("\n")
  615. def waitCount(self):
  616. " return number of waiting jobs "
  617. sql = 'SELECT count(*) FROM queue WHERE isRunning=0'
  618. return self.conn.execute(sql).fetchone()[0]
  619. def popJob(self):
  620. " return (jobType, jobId, params) of first waiting job and set it to running state "
  621. print('pop job<br>')
  622. try:
  623. self.conn.execute('BEGIN IMMEDIATE')
  624. sql = 'SELECT jobType, jobId, paramStr FROM queue WHERE isRunning=0 ORDER BY lastUpdate LIMIT 1'
  625. row = self.conn.execute(sql).fetchone()
  626. if row is None:
  627. logging.debug("popJob: no waiting jobs")
  628. self.conn.commit()
  629. return None, None, None
  630. jobType, jobId, paramStr = row
  631. sql = 'UPDATE queue SET isRunning=1 where jobId=?'
  632. self.conn.execute(sql, (jobId,))
  633. self.conn.commit()
  634. except:
  635. self.conn.rollback()
  636. raise
  637. return jobType, jobId, paramStr
  638. def clearJobs(self):
  639. " clear the job table, removing running jobs, too "
  640. self.conn.execute("DELETE from queue")
  641. def close(self):
  642. " "
  643. self.conn.close()
  644. # ====== FUNCTIONS =====
  645. contentLineDone = False
  646. # the queue workers should be able to never abort
  647. doAbort = True
  648. def getTwoBitFname(db):
  649. " return the name of the twoBit file for a genome "
  650. # at UCSC, try to use local disk, if possible
  651. locPath = join("/scratch", "data", db, db+".2bit")
  652. if isfile(locPath):
  653. return locPath
  654. path = join(genomesDir, db, db+".2bit")
  655. return path
  656. def errAbort(msg, isWarn=False):
  657. " print err msg and exit "
  658. if commandLineMode:
  659. raise Exception(msg)
  660. if not contentLineDone:
  661. print("Content-type: text/html\n")
  662. print('<div style="position: absolute; padding: 10px; left: 100; top: 100; border: 10px solid black; background-color: white; text-align:left; width: 800px; font-size: 18px">')
  663. if isWarn:
  664. print("<strong>Warning:</strong><p> ")
  665. else:
  666. print("<strong>Error:</strong><p> ")
  667. print((msg+"<p>"))
  668. print(("If you think this is a bug or you have any other suggestions, please do not hesitate to contact us %s<p>" % contactEmail))
  669. if isWarn:
  670. print("In the email, please also send us the full URL of the page.")
  671. else:
  672. print("Please also send us the full URL of the page where you see the error. Thanks!")
  673. print('</div>')
  674. if doAbort:
  675. sys.exit(0) # cgi must not exit with 1
  676. # allow only dashes, digits, characters, underscores and colons in the CGI parameters
  677. # and +
  678. notOkChars = re.compile(r'[^+a-zA-Z0-9/:\n\r_. -]')
  679. def checkVal(key, inStr):
  680. """ remove special characters from input string, to protect against injection attacks """
  681. if key!="geneIds":
  682. if len(inStr) > 10000:
  683. errAbort("input parameter %s is too long" % key)
  684. else:
  685. if len(inStr) > 100000:
  686. errAbort("Pasting more than tens of thousands of gene IDs makes little sense. Copy/paste error?")
  687. matchObj =notOkChars.search(inStr)
  688. if matchObj!=None:
  689. errAbort("input parameter %s contains an invalid character %s (ASCII %d)" % (key, repr(matchObj.group()), ord(matchObj.group())))
  690. return inStr
  691. def cgiGetParams():
  692. " get CGI parameters and return as dict "
  693. form = cgi.FieldStorage()
  694. global cgiParams
  695. cgiParams = {}
  696. # parameters are:
  697. #"pamId", "batchId", "pam", "seq", "org", "download", "sortBy", "format", "ajax
  698. for key in list(form.keys()):
  699. val = form.getfirst(key)
  700. if val!=None:
  701. # "seq" is cleaned by cleanSeq later
  702. val = urllib.parse.unquote(val)
  703. if key not in ["seq", "name"]:
  704. checkVal(key, val)
  705. cgiParams[key] = val
  706. if "pam" in cgiParams:
  707. legalChars = set("ACTGNMKRYVBE120345/-")
  708. illegalChars = set(cgiParams["pam"])-legalChars
  709. if len(illegalChars)!=0:
  710. errAbort("Illegal character in PAM-sequence. Only %s are allowed."+"".join(legalChars))
  711. if "batchId" in cgiParams:
  712. batchId = cgiParams["batchId"]
  713. if not batchId.isalnum() or len(batchId) > 30:
  714. errAbort("Invalid batchId")
  715. return cgiParams
  716. def cgiGetStr(params, argName, default=None):
  717. val = params.get(argName, None)
  718. if val==None and default==None:
  719. errAbort("'%s' parameter must be specified" % argName)
  720. if val==None:
  721. return default
  722. return val
  723. def cgiGetNum(params, argName, default):
  724. " get CGI parameter which must be a number "
  725. val = params.get(argName, None)
  726. if val==None:
  727. return default
  728. if not val.isdigit():
  729. errAbort("'%s' parameter must be a number" % argName)
  730. val = int(val)
  731. return val
  732. transTab = str.maketrans("-=/+_", "abcde")
  733. def makeTempBase(seq, org, pam, batchName):
  734. "create the base name of temp files using a hash function and some prettyfication "
  735. hasher = hashlib.sha1(seq.encode("latin1")+org.encode("latin1")+pam.encode("latin1")+batchName.encode("latin1"))
  736. shortHash = hasher.digest()[0:20]
  737. batchId = base64.urlsafe_b64encode(shortHash).decode('latin1').translate(transTab)[:20]
  738. return batchId
  739. def makeTempFile(prefix, suffix):
  740. " return a temporary file that is deleted upon exit, unless DEBUG is set "
  741. if DEBUG:
  742. fname = join("/tmp", prefix+suffix)
  743. fh = open(fname, "wt")
  744. else:
  745. fh = tempfile.NamedTemporaryFile(mode="wt", dir=TEMPDIR, prefix="primer3In", suffix=".txt")
  746. return fh
  747. def pamIsCpf1(pam):
  748. " if you change this, also change bin/filterFaToBed and bin/samToBed!!! "
  749. return (pam in ["TNN", "TTN", "TTTN", "TYCV", "TATV", "TTTV", "TTTR", "ATTN", "TTTA", "TCTA", "TCCA", "CCCA", "YTTV", "TTYN"])
  750. # test : modify cleavage sites for hfCas12max according to Synthego specifications
  751. def pamIsCas12max(pam):
  752. return (pam in ["TNN", "TTN"])
  753. def pamIsCasX(pam):
  754. " if you change this, also change bin/filterFaToBed and bin/samToBed!!! "
  755. return (pam in ["TTCN"])
  756. def pamIsSaCas9(pam):
  757. " only used for notes and efficiency scores, unlike its Cpf1 cousin function "
  758. return (pam.split("-")[0] in ["NNGRRT", "NNNRRT"])
  759. def isSlowPam(pam):
  760. " do not allow input sequences > 500 bp "
  761. if pamIsXCas9(pam) or pam=="TTYN" or pam=="NNG" or pam=="TNN":
  762. return True
  763. else:
  764. return False
  765. def pamIsXCas9(pam):
  766. " "
  767. return (pam in ["NGK", "NGN"])
  768. def pamIsSpCas9(pam):
  769. " only used for notes and efficiency scores, unlike its Cpf1 cousin function "
  770. return (pam in ["NGG", "NGA", "NGCG"])
  771. def saveSeqOrgPamToCookies(seq, org, pam):
  772. " create a cookie with seq, org and pam and print it"
  773. cookies=http.cookies.SimpleCookie()
  774. expires = 365 * 24 * 60 * 60
  775. if len(seq)<3000:
  776. cookies['lastseq'] = seq
  777. else:
  778. cookies['lastseq'] = "(last sequence was too long, could not be saved in Internet Browser cookie)"
  779. cookies['lastseq']['expires'] = expires
  780. cookies['lastorg'] = org
  781. cookies['lastorg']['expires'] = expires
  782. cookies['lastpam'] = pam
  783. cookies['lastpam']['expires'] = expires
  784. print(cookies)
  785. def debug(msg):
  786. if commandLineMode:
  787. logging.debug(msg)
  788. elif DEBUG:
  789. print(msg)
  790. print("<br>")
  791. def gcContent(seq):
  792. " return GC content as a float "
  793. c = 0
  794. for x in seq:
  795. if x in ["G","C"]:
  796. c+= 1
  797. return (float(c)/len(seq))
  798. def findPat(seq, pat):
  799. """ yield positions where pat matches seq, stupid brute force search
  800. """
  801. seq = seq.upper()
  802. pat = pat.upper()
  803. patLen = len(pat)
  804. for i in range(0, len(seq)-patLen+1):
  805. subseq = seq[i:i+patLen]
  806. if patMatch(subseq, pat):
  807. yield i
  808. def rndSeq(seqLen):
  809. " return random seq "
  810. seq = []
  811. alf = "ACTG"
  812. for i in range(0, seqLen):
  813. seq.append(alf[random.randint(0,3)])
  814. return "".join(seq)
  815. def cleanSeq(seq, db):
  816. """ remove fasta header, check seq for illegal chars and return (filtered
  817. seq, user message) special value "random" returns a random sequence.
  818. """
  819. #print repr(seq)
  820. if seq.startswith("random"):
  821. seq = rndSeq(800)
  822. lines = seq.strip().splitlines()
  823. #print "<br>"
  824. #print "before fasta cleaning", "|".join(lines)
  825. if len(lines)>0 and lines[0].startswith(">"):
  826. line1 = lines.pop(0)
  827. #print "<br>"
  828. #print "after fasta cleaning", "|".join(lines)
  829. #print "<br>"
  830. newSeq = []
  831. nCount = 0
  832. for l in lines:
  833. if len(l)==0:
  834. continue
  835. for c in l:
  836. if c not in "actgACTGNn":
  837. nCount +=1
  838. else:
  839. newSeq.append(c)
  840. seq = "".join(newSeq)
  841. msgs = []
  842. tooLongHint = """
  843. Please split your input sequence into shorter sequences or use
  844. the <a href='downloads/'>stand-alone version</a> on your own Linux or Mac server to process longer sequences in batch.<br>
  845. """
  846. if len(seq)>MAXSEQLEN and db!="noGenome":
  847. errMsg = "<strong>Sorry, this tool cannot handle sequences longer than %d bp</strong><br>" % (MAXSEQLEN)
  848. errAbort(errMsg+tooLongHint)
  849. if len(seq)>MAXSEQLEN_NOGENOME and db=="noGenome":
  850. errMsg = "<strong>Sorry, this tool cannot handle sequences longer than %d bp when using the 'No Genome' option.</strong><br>" % (MAXSEQLEN_NOGENOME)
  851. errAbort(errMsg+tooLongHint)
  852. if nCount!=0:
  853. msgs.append("Sequence contained %d non-ACTGN letters. They were removed." % nCount)
  854. return seq, "<br>".join(msgs)
  855. revTbl = {'A' : 'T', 'C' : 'G', 'G' : 'C', 'T' : 'A', 'N' : 'N' , 'M' : 'K', 'K' : 'M',
  856. "R" : "Y" , "Y":"R" , "g":"c", "a":"t", "c":"g","t":"a", "n":"n", "V" : "B", "v":"b",
  857. "B" : "V", "b": "v", "W" : "W", "w" : "w"}
  858. def revComp(seq):
  859. " rev-comp a dna sequence with UIPAC characters "
  860. newSeq = []
  861. for c in reversed(seq):
  862. newSeq.append(revTbl[c])
  863. return "".join(newSeq)
  864. def docTestInit(isCpf1, guideLen):
  865. global pamIsFirst
  866. global GUIDELEN
  867. pamIsFirst=isCpf1
  868. GUIDELEN=guideLen
  869. def findPams (seq, pam, strand, startDict, endSet):
  870. """ return two values: dict with pos -> strand of PAM and set of end positions of PAMs
  871. Makes sure to return only values with at least GUIDELEN bp left (if strand "+") or to the
  872. right of the match (if strand "-")
  873. If the PAM is cpf1, then this is inversed: pos-strand matches must have at least GUIDELEN
  874. basepairs to the right, neg-strand matches must have at least GUIDELEN bp on their left
  875. >>> docTestInit(False, 20)
  876. >>> findPams("GGGGGGGGGGGGGGGGGGGGGGG", "NGG", "+", {}, set())
  877. ({20: '+'}, {23})
  878. >>> findPams("CCAGCCCCCCCCCCCCCCCCCCC", "CCA", "-", {}, set())
  879. ({0: '-'}, {3})
  880. >>> docTestInit(True, 20)
  881. >>> findPams("TTTNCCCCCCCCCCCCCCCCCTTTN", "TTTN", "+", {}, set())
  882. ({0: '+'}, {4})
  883. >>> docTestInit(False, 20)
  884. >>> findPams("CCCCCCCCCCCCCCCCCCCCCAAAA", "NAA", "-", {}, set())
  885. ({}, set())
  886. >>> findPams("AAACCCCCCCCCCCCCCCCCCCCC", "NAA", "-", {}, set())
  887. ({0: '-'}, {3})
  888. >>> findPams("CCCCCCCCCCCCCCCCCCCCCCCCCAA", "NAA", "-", {}, set())
  889. ({}, set())
  890. >>> findPams("GTTGTGTTTTACAATGCAGAGAGTGGAGGATGCTTTTTATACATTGGTGAGAGAGATCCGACAGTACAGATTGAAAAAAATCAGCAAAGAAGAAAAGACTCCTGGCTGTGTGAAAATTAAAAAATGCGTTATAATGTAATCTGGTAAGTTGAGCATATTCATTCTGGTACAAAGCAGATGTCTTCAGAGGTAACA", "TATV", "-", {}, set())
  891. ({37: '-', 129: '-'}, {41, 133})
  892. >>> findPams("GTTGTGTTTTACAATGCAGAGAGTGGAGGATGCTTTTTATACATTGGTGAGAGAGATCCGACAGTACAGATTGAAAAAAATCAGCAAAGAAGAAAAGACTCCTGGCTGTGTGAAAATTAAAAAATGCGTTATAATGTAATCTGGTAAGTTGAGCATATTCATTCTGGTACAAAGCAGATGTCTTCAGAGGTAACA", "TATV", "+", {}, set())
  893. ({37: '+', 129: '+'}, {41, 133})
  894. """
  895. assert(pamIsFirst is not None)
  896. if pamIsFirst:
  897. maxPosPlus = len(seq)-(GUIDELEN+len(pam))
  898. minPosMinus = GUIDELEN
  899. else:
  900. # -------------------
  901. # OKOKOKOKOK
  902. minPosPlus = GUIDELEN
  903. # -------------------
  904. # OKOKOKOKOK
  905. maxPosMinus = len(seq)-(GUIDELEN+len(pam))
  906. #print "new search", seq, pam, "minPosPlus=",minPosPlus, "guideLen=", GUIDELEN, "<br>"
  907. for start in findPat(seq, pam):
  908. if pamIsFirst:
  909. # need enough flanking seq on one side
  910. #return("Cpf1 mode found", start,"<br>")
  911. if strand == "+" and start > maxPosPlus:
  912. continue
  913. if strand == "-" and start < minPosMinus:
  914. continue
  915. else:
  916. # return("non-Cpf1 mode found", start,"<br>")
  917. if strand=="+" and start < minPosPlus:
  918. continue
  919. if strand=="-" and start > maxPosMinus:
  920. continue
  921. #print "match", strand, start, end, "<br>"
  922. startDict[start] = strand
  923. end = start+len(pam)
  924. endSet.add(end)
  925. return startDict, endSet
  926. def rulerString(maxLen):
  927. " return line with positions every 10 chars "
  928. texts = []
  929. for i in range(0, maxLen, 10):
  930. numStr = str(i)
  931. texts.append(numStr)
  932. spacer = "".join([" "]*(10-len(numStr)))
  933. texts.append(spacer)
  934. return "".join(texts)
  935. def varDictToHtml(varDict, seq, varShortLabel):
  936. " make a list of one html string per position in the sequence "
  937. if varDict is None:
  938. return None
  939. varHtmls = []
  940. for i in range(0, len(seq)):
  941. if not i in varDict:
  942. varHtmls.append(".")
  943. else:
  944. varHooverLines = []
  945. showStar = False # show a star if change is non-simple SNP
  946. varInfos = varDict[i]
  947. for chrom, pos, refAll, altAll, infoDict in varInfos:
  948. varHooverLines.append("%s: %s &rarr; %s<br>" % (varShortLabel, refAll, altAll))
  949. if "freq" in infoDict:
  950. varHooverLines.append("&nbsp;<b>Freq:</b> %s<br>" % infoDict["freq"])
  951. #if "dbg" in infoDict:
  952. #varHooverLines.append("%s<br>" % infoDict["dbg"])
  953. if "varId" in infoDict:
  954. varHooverLines.append("&nbsp;<b>ID:</b> %s<br>" % infoDict["varId"])
  955. if "ExAC" in varDict["label"]:
  956. endPos = int(pos)+len(refAll)
  957. varHooverLines.append('&nbsp;<a target=_blank href="http://exac.broadinstitute.org/region/%s-%s-%d">ExAC Browser</a><br>' % (chrom, pos, endPos))
  958. if len(refAll)!=1 or len(altAll)!=1:
  959. showStar = True
  960. if len(varInfos)!=1:
  961. showStar = True
  962. varDesc = "".join(varHooverLines)
  963. if showStar:
  964. dispChar = "*"
  965. else:
  966. dispChar = altAll
  967. varHtmls.append("<u class='tooltipsterInteract' title='%s'>%s</u>" % (varDesc, dispChar))
  968. return varHtmls
  969. def cssClassesFromSeq(guideSeq, suffix=""):
  970. " The CSS class of guide row and links in seq viewer depend on the first nucl of guide "
  971. classNames = ["guideRow"]
  972. if guideSeq[0].upper()!="G":
  973. classNames.append("guideRowNoPrefixG"+suffix)
  974. if not guideSeq.startswith("GG"):
  975. classNames.append("guideRowNoPrefixGG"+suffix)
  976. if guideSeq[0].upper()!="A":
  977. classNames.append("guideRowNoPrefixA"+suffix)
  978. classStr = " ".join(classNames)
  979. return classStr
  980. def buildCodonTable():
  981. " from http://www.petercollingridge.co.uk/tutorials/bioinformatics/codon-table/ "
  982. bases = "TCAG"
  983. codons = [a + b + c for a in bases for b in bases for c in bases]
  984. amino_acids = 'FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG'
  985. codon_table = dict(list(zip(codons, amino_acids)))
  986. return codon_table
  987. def buildOneToThree():
  988. " return one-letter -> three-letter conversion table for amino acids "
  989. oneToThree = \
  990. {'C':'Cys', 'D':'Asp', 'S':'Ser', 'Q':'Gln', 'K':'Lys',
  991. 'I':'Ile', 'P':'Pro', 'T':'Thr', 'F':'Phe', 'N':'Asn',
  992. 'G':'Gly', 'H':'His', 'L':'Leu', 'R':'Arg', 'W':'Trp',
  993. 'A':'Ala', 'V':'Val', 'E':'Glu', 'Y':'Tyr', 'M':'Met',
  994. 'U':'Sec', '*':'Stop',
  995. 'X':'Stop', # is this really used like that?
  996. 'Z':'Glx', # special case: asparagine or aspartic acid
  997. 'B':'Asx' # special case: glutamine or glutamic acid
  998. }
  999. return oneToThree
  1000. def makeExonLines(exonInfo, seq, selTransId):
  1001. """ create text that draws exons, input is transId -> (exonNumber, exStart, exEnd, exFrame).
  1002. returns a list of (transId (=label), symbol (=mouseover), ASCII-line) """
  1003. lines = []
  1004. #maxLabelLen = 0
  1005. codonTable = buildCodonTable()
  1006. #oneToThree = buildOneToThree()
  1007. seqLen = len(seq)
  1008. seq = seq.upper()
  1009. for (transId, symbol), exRows in exonInfo.items():
  1010. if selTransId!="allTrans" and transId!=selTransId:
  1011. continue
  1012. line = [" "]*seqLen
  1013. mouseOvers = {} # position -> mouseOver-text or null, for end of mouse over
  1014. for exIdx, (exNum, exStart, exEnd, exFrame, nextFrame, exStrand) in enumerate(exRows):
  1015. if exFrame==-1:
  1016. for i in range(exStart, exEnd):
  1017. line[i]="="
  1018. exonLabel = "noncoding"
  1019. if (exEnd-exStart)>len(exonLabel)+4:
  1020. # center the exon label on the exon
  1021. mid = exStart+int((exEnd-exStart)*0.5)
  1022. halfLen = int(len(exonLabel)*0.5)
  1023. labStart = mid-halfLen
  1024. for i in range(0, len(exonLabel)):
  1025. line[labStart+i] = exonLabel[i]
  1026. line[labStart-1] = " "
  1027. line[labStart+len(exonLabel)] = " "
  1028. else:
  1029. exonDesc = "gene %s<br>transcript %s<br>exon %d<br>start phase %s" % (symbol, transId, exNum+1, exFrame)
  1030. if nextFrame is not None:
  1031. exonDesc += "<br>end phase %s" % nextFrame
  1032. if (exFrame+nextFrame) % 3 == 0:
  1033. exonDesc += "<br>Removing the exon retains the reading frame"
  1034. else:
  1035. exonDesc += "<br>Removing the exon will destroy the reading frame"
  1036. mouseOvers[exStart] = exonDesc
  1037. mouseOvers[exEnd] = None
  1038. for i in range(exStart, exStart+exFrame):
  1039. line[i] = "-"
  1040. for i in range(exStart+exFrame, exEnd, 3):
  1041. codon = seq[i:i+3]
  1042. if len(codon)==3:
  1043. shortAa = codonTable[codon]
  1044. if exStrand=="+":
  1045. longAa = shortAa+"]]"
  1046. else:
  1047. longAa = "[["+shortAa # highlighting rev. dir. more
  1048. else:
  1049. # codon is split by splice site
  1050. longAa = "-"
  1051. for j in range(0, len(longAa)):
  1052. line[i+j] = longAa[j]
  1053. # now merge the mouse overs as span tags into the ASCII line
  1054. newLine = []
  1055. for pos, char in enumerate(line):
  1056. if pos in mouseOvers:
  1057. overString = mouseOvers[pos]
  1058. if overString is not None:
  1059. newLine.append("<span class='tooltipsterInteract' title='%s'>" % overString)
  1060. if char=="-":
  1061. newLine.append("<span class='tooltipsterInteract' title='This codon goes over a splice site. The nucleotides of the split codon are not translated to amino acids but shown as dashes.'>-</span>")
  1062. else:
  1063. newLine.append(char)
  1064. if overString is None:
  1065. newLine.append("</span>")
  1066. else:
  1067. newLine.append(char)
  1068. # and fix up the < signs
  1069. newLineStr = "".join(newLine)
  1070. newLineStr = newLineStr.replace("[[", "<lc>&lt;&lt;</lc>")
  1071. newLineStr = newLineStr.replace("]]", "<lc>&gt;&gt;</lc>")
  1072. lines.append((symbol, transId, newLineStr))
  1073. #maxLabelLen = max(maxLabelLen, len(transId))
  1074. return lines
  1075. def getGeneModels(org):
  1076. " read possible gene models for org and return as list (name, desc) or None if no gene models "
  1077. mask = join(genomesDir, org, "*.bb")
  1078. fnames = glob.glob(mask)
  1079. descFname = join(genomesDir, org, "genes.tsv")
  1080. if not isfile(descFname):
  1081. return None
  1082. geneDescs = {}
  1083. for line in open(descFname):
  1084. fname, desc = line.split(maxsplit=1)
  1085. geneDescs[fname] = desc
  1086. ret = []
  1087. for fname in fnames:
  1088. baseName = basename(fname)
  1089. name = baseName.split('.')[0]
  1090. desc = geneDescs.get(baseName, name)
  1091. ret.append((name, desc))
  1092. return ret
  1093. def getSelGeneModel(org):
  1094. " return (list of (name, desc) of models, selected gene model name) "
  1095. geneModels = getGeneModels(org)
  1096. selGeneModel = None
  1097. selTransId = None
  1098. if geneModels:
  1099. #selGeneModel = cgiParams.get("geneModel", geneModels[0][0])
  1100. selGeneModel = cgiParams.get("geneModel", "noGenes")
  1101. geneModels.insert(0, ("noGenes", "Do not show"))
  1102. possNames = [x for x,y in geneModels]
  1103. if not selGeneModel in possNames:
  1104. errAbort("The gene model name specified with the argument geneModel is invalid")
  1105. selTransId = cgiParams.get("transId", "allTrans")
  1106. return geneModels, selGeneModel, selTransId
  1107. def printSeqForCopy(seq):
  1108. " print a hidden text area so we can copy the sequence to the clipboard "
  1109. print('<input id="seqAsText" type="text" style="display:none">')
  1110. print(seq)
  1111. print("</input>")
  1112. def calcKomorScore(guideSeq, pos):
  1113. " return base editing score given the guide sequence and the position "
  1114. return pos/7.0 # temporary hack
  1115. def makeEditLines(seq, pamSeqs, winStart, winEnd, guideScores):
  1116. " create the lines that show the possible baseEditor edits "
  1117. editInfos = []
  1118. for i in range(0, len(seq)):
  1119. editInfos.append(defaultdict(list))
  1120. upSeq = seq.upper()
  1121. for pamId, pamStart, guideStart, strand, guideSeq, pamSeq, pamPlusSeq in pamSeqs:
  1122. specScore = guideScores[pamId]
  1123. if strand=="+":
  1124. fromPos = guideStart+winStart
  1125. toPos = guideStart+winEnd
  1126. fromNucl = "C"
  1127. toNucl = "T"
  1128. else:
  1129. guideEnd = guideStart+GUIDELEN
  1130. fromPos = guideEnd-winEnd
  1131. toPos = guideEnd-winStart
  1132. fromNucl = "G"
  1133. toNucl = "A"
  1134. for pos in range(fromPos, toPos):
  1135. # position of mutated nucl on forw strand guide
  1136. if strand=="+":
  1137. mutPos = pos-guideStart
  1138. else:
  1139. mutPos = GUIDELEN - (pos - guideStart) - 1
  1140. if upSeq[pos]==fromNucl:
  1141. beScore = calcKomorScore(guideSeq, mutPos)
  1142. editInfos[pos][toNucl].append((pamId, guideSeq, pamSeq, mutPos, beScore, specScore))
  1143. altNucls = ["A", "T"]
  1144. editLabels = []
  1145. for an in altNucls:
  1146. editLabels.append("Edits to "+an)
  1147. editLines = []
  1148. for i in range(0, len(altNucls)):
  1149. editLines.append([" "]*len(seq))
  1150. # rearrange into lines of text + JSON
  1151. jsonData = defaultdict(list)
  1152. for pos, eiDict in enumerate(editInfos):
  1153. if not eiDict:
  1154. continue
  1155. jsonData[pos] = eiDict
  1156. for nucl, guideData in eiDict.items():
  1157. yPos = altNucls.index(nucl)
  1158. editLines[yPos][pos] = "<d pos=%d>%s</d>" % (pos, nucl)
  1159. ret = []
  1160. for label, lineChars in zip(editLabels, editLines):
  1161. ret.append( (label, None, "".join(lineChars)) )
  1162. return ret, jsonData
  1163. def makePamLines(lines, maxY, pamIdToSeq, guideScores):
  1164. for y in range(0, maxY+1):
  1165. texts = []
  1166. lastEnd = 0
  1167. for start, end, name, strand, pamId in lines[y]:
  1168. guideSeq = pamIdToSeq.get(pamId)
  1169. if guideSeq==None:
  1170. # when there is an N in the guide, the PAM is valid, but the guide is not
  1171. continue
  1172. classStr = cssClassesFromSeq(guideSeq, suffix="Seq")
  1173. spacer = "".join([" "]*((start-lastEnd)))
  1174. lastEnd = end
  1175. texts.append(spacer)
  1176. score = guideScores[pamId]
  1177. # XX How can this happen for non-Cpf1 enzymes? Can this ever happen?
  1178. if score is None and not pamIsFirst:
  1179. continue
  1180. color = scoreToColor(score)[0]
  1181. texts.append('''<a class='%s' style="text-shadow: 1px 1px 1px #bbb; color: %s" id="list%s" href="#%s">''' % (classStr, color, pamId,pamId))
  1182. texts.append(name)
  1183. texts.append("</a>")
  1184. yield ("", None, ''.join(texts))
  1185. def getBeWin(winVal):
  1186. " return (start, end) of base editor window given CGI variable "
  1187. fs = winVal.split("-")
  1188. if len(fs)!=2:
  1189. errAbort("parameter beWin must contain only one dash")
  1190. start = fs[0].strip()
  1191. end = fs[1].strip()
  1192. if not start.isdigit() or not end.isdigit():
  1193. errAbort("parameter beWin must be two dash-separated numbers")
  1194. start = int(start)
  1195. end = int(end)
  1196. return start, end
  1197. def printLines(lines, labelLen):
  1198. " print list of (label, string) such that label is at least labelLen characters long "
  1199. for label, mouseOver, line in lines:
  1200. if mouseOver is not None:
  1201. labelStr =('<span class="tooltipsterInteract" title="{:s}">{:'+str(labelLen)+'s} </span>').format(label, mouseOver)
  1202. else:
  1203. labelStr = ('{:'+str(labelLen)+'s} ').format(label)
  1204. print((labelStr), end=' ')
  1205. print(line)
  1206. def getMaxLen(lines):
  1207. " given a list of tuples where first element is the label, return the longest label len "
  1208. maxLen = 0
  1209. for l in lines:
  1210. label = l[0]
  1211. maxLen = max(maxLen, len(label))
  1212. return maxLen
  1213. def printJson(name, obj):
  1214. print("<script>")
  1215. print((name), end=' ')
  1216. print(("="), end=' ')
  1217. print((json.dumps(obj)))
  1218. print("</script>")
  1219. def showSeqAndPams(org, seq, startDict, pam, guideScores, varHtmls, varDbs, varDb, minFreq, position, pamIdToSeq):
  1220. " show the sequence and the PAM sites underneath in a sequence viewer "
  1221. pamSeqs = list(flankSeqIter(seq, startDict, len(pam), True))
  1222. lines, maxY = distrOnLines(seq.upper(), startDict, len(pam), pam)
  1223. posLabel = "Position"
  1224. varLabel = "Variants"
  1225. seqLabel = "Sequence"
  1226. exonLabelLen = 0
  1227. editLines = []
  1228. exonLines = []
  1229. geneModels, selGeneModel, selTransId = getSelGeneModel(org)
  1230. #selGeneModel = None
  1231. #geneModels = None
  1232. if baseEditor:
  1233. beWinStart, beWinEnd = getBeWin(cgiParams.get("beWin", DEFAULTBEWIN))
  1234. editLines, jsonData = makeEditLines(seq, pamSeqs, beWinStart, beWinEnd, guideScores)
  1235. printJson("editData", jsonData)
  1236. pamLines = list(makePamLines(lines, maxY, pamIdToSeq, guideScores))
  1237. labelLen = max(len(varLabel), len(seqLabel), len(posLabel), getMaxLen(pamLines))
  1238. if selGeneModel!=None:
  1239. exonInfo, maxTransIdLen = getExonInfo(org, selGeneModel, position)
  1240. labelLen = max(labelLen, maxTransIdLen)
  1241. if baseEditor:
  1242. labelLen = max(labelLen, getMaxLen(editLines))
  1243. if selGeneModel:
  1244. labelLen = max(labelLen, exonLabelLen)
  1245. print("<div class='substep'>")
  1246. print('<a id="seqStart"></a>')
  1247. print("Your input sequence is %d bp long. It contains %d possible guide sequences.<br>" % (len(seq), len(guideScores)))
  1248. if not pamIsFirst:
  1249. print("Shown below are their PAM sites and the expected cleavage position located -3bp 5' of the PAM site.<br>")
  1250. print("Click on a match for the PAM %s below to show its %d bp-long guide sequence. " % (pam, GUIDELEN))
  1251. print("(Need help? Look at the <a target=_blank href='manual/#annotseq'>CRISPOR manual</a>)<br>")
  1252. print('''Colors <span style="color:#32cd32; text-shadow: 1px 1px 1px #bbb">green</span>, <span style="color:#ffff00; text-shadow: 1px 1px 1px #888">yellow</span> and <span style="text-shadow: 1px 1px 1px #f01; color:#aa0014">red</span> indicate high, medium and low specificity of the PAM's guide sequence in the genome.<p>''')
  1253. else:
  1254. print("Click on a match for the PAM %s below to show its %d bp-long guide sequence.<br>" % (pam, GUIDELEN))
  1255. if baseEditor or varDb or selGeneModel:
  1256. print(("""<form style="display:inline" id="paramForm" action="%s" method="GET">""" % basename(__file__)))
  1257. if geneModels:
  1258. print ("Gene Models:")
  1259. printDropDown("geneModel", geneModels, selGeneModel, style="width:20em")
  1260. if selGeneModel!="noGenes":
  1261. print ("Transcript:")
  1262. transIdInfo = [("allTrans", "All Transcripts")]
  1263. for transId, sym in list(exonInfo.keys()):
  1264. transIdInfo.append( (transId, sym+" / "+transId) )
  1265. printDropDown("transId", transIdInfo, selTransId, style="width:20em")
  1266. # XX XXXXXX
  1267. exonLines = makeExonLines(exonInfo, seq, selTransId)
  1268. #exonLines = []
  1269. print("""<input style="height:18px;margin:0px;font-size:10px;line-height:normal" type="submit" name="submit" value="Update">""")
  1270. print("""<br>""")
  1271. if baseEditor:
  1272. print ("Base Editor modification window:")
  1273. print(("""<input type="text" name="beWin" size="10" value="%s">""" % DEFAULTBEWIN))
  1274. print("""<input style="height:18px;margin:0px;font-size:10px;line-height:normal" type="submit" name="submit" value="Update">""")
  1275. print("<br>")
  1276. if varDb is not None:
  1277. print ("Variant database:")
  1278. varDbList = [(b,c) for a,b,c,d in varDbs] # only keep fname+label
  1279. printDropDown("varDb", varDbList, varDb)
  1280. if minFreq==0.0:
  1281. minFreq="0.0"
  1282. else:
  1283. minFreq = str(minFreq)
  1284. # pull out the hasAF field for this varDb
  1285. varDbHasAF = False
  1286. for shortLabel, fname, desc, hasAF in varDbs:
  1287. if fname==varDb:
  1288. varDbHasAF = hasAF
  1289. break
  1290. if varDbHasAF:
  1291. print("""&nbsp; Min. frequency: """)
  1292. print(("""<input type="text" name="minFreq" size="8" value="%s">""" % minFreq))
  1293. print("""<input style="height:18px;margin:0px;font-size:10px;line-height:normal" type="submit" name="submit" value="Update">""")
  1294. print(("<small style='margin-left:30px'><a href='mailto:%s'>Missing a variant database? We can add it.</a></small>" % contactEmail))
  1295. if position=="?":
  1296. print("<small style=''>Input sequence not in genome, cannot show genome variants.</small>")
  1297. elif varDb is None:
  1298. print(("<small style=''><a href='mailto:%s'>Suggest a genome variants database to show on this page</a></small>" % contactEmail))
  1299. print("</div>")
  1300. print('''<div class="blueHighlight" style="text-align: left; overflow-x:scroll; width:98vw; background:#DDDDDD; border-style: solid; border-width: 1px">''')
  1301. print('''<pre style="font-family: Source Code Pro; font-size: 80%; display:inline; line-height: 0.95em; text-align:left">''')
  1302. print(('{:'+str(labelLen)+'s} ').format(posLabel), end=' ')
  1303. print(rulerString(len(seq)))
  1304. if varHtmls is not None:
  1305. print(('{:'+str(labelLen)+'s} ').format(varLabel), end=' ')
  1306. print("".join(varHtmls))
  1307. print(('{:'+str(labelLen)+'s} ').format(seqLabel), end=' ')
  1308. print (seq)
  1309. printLines(exonLines, labelLen)
  1310. if baseEditor:
  1311. printLines(editLines, labelLen)
  1312. printLines(pamLines, labelLen)
  1313. print("</pre><br>")
  1314. print('''</div>''')
  1315. #printSeqForCopy(seq)
  1316. if pamIsCas12max(pam):
  1317. print('<div style="line-height: 1.0; padding-top: 5px; font-size: 15px">Cpf1 has a staggered site: cleavage occurs between the 14th and 16th base on the non-targeted strand (indicated by "\\" in the schema above). Cleavage mostly occurs after the 24rd base on the targeted strand (indicated by "/" in the schema above). See on <a target=_blank href="https://www.synthego.com/products/nuclease/hfcas12max-hifi">Synthego</a></div>')
  1318. elif pamIsCpf1(pam):
  1319. print('<div style="line-height: 1.0; padding-top: 5px; font-size: 15px">Cpf1 has a staggered site: cleavage occurs usually - but not always - after the 18th base on the non-targeted strand which has the TTTV PAM motif (indicated by "\\" in the schema above). Cleavage mostly occurs after the 23rd base on the targeted strand which has the AAAN motif (indicated by "/" in the schema above). See <a target=_blank href="http://www.sciencedirect.com/science/article/pii/S0092867415012003">Zetsche et al 2015</a>, in particular <a target=_blank href="http://www.sciencedirect.com/science?_ob=MiamiCaptionURL&_method=retrieve&_eid=1-s2.0-S0092867415012003&_image=1-s2.0-S0092867415012003-gr3.jpg&_cid=272196&_explode=defaultEXP_LIST&_idxType=defaultREF_WORK_INDEX_TYPE&_alpha=defaultALPHA&_ba=&_rdoc=1&_fmt=FULL&_issn=00928674&_pii=S0092867415012003&md5=11771263f3e390e444320cacbcfae323">Fig 3</a>.</div>')
  1320. elif pamIsCasX(pam):
  1321. print('<div style="line-height: 1.0; padding-top: 5px; font-size: 15px">We have no description yet on how exactly the CasX cleavage looks like. Please contact [email hidden] if you have an idea how to describe the cleavage site.</div>')
  1322. def iterOneDelSeqs(seq):
  1323. """ given a seq, create versions with each bp removed. Avoid duplicates
  1324. yields (delPos, seq)
  1325. >>> list(iterOneDelSeqs("AATGG"))
  1326. [(0, 'ATGG'), (2, 'AAGG'), (3, 'AATG')]
  1327. """
  1328. doneSeqs = set()
  1329. for i in range(0, len(seq)):
  1330. delSeq = seq[:i]+seq[i+1:]
  1331. if delSeq not in doneSeqs:
  1332. yield i, delSeq
  1333. doneSeqs.add(delSeq)
  1334. def flankSeqIter(seq, startDict, pamLen, doFilterNs):
  1335. """ given a seq and dictionary of pamPos -> strand and the length of the pamSite
  1336. yield tuples of (name, pamStart, guideStart, strand, flankSeq, pamSeq)
  1337. flankSeq is the guide sequence (=flanking the PAM).
  1338. if doFilterNs is set, will not return any sequences that contain an N character
  1339. pamPlusSeq are the 5bp after the PAM. If not enough space, pamPlusSeq is None
  1340. """
  1341. startList = sorted(startDict.keys())
  1342. for pamStart in startList:
  1343. strand = startDict[pamStart]
  1344. pamPlusSeq = None
  1345. if pamIsFirst: # Cpf1: get the sequence to the right of the PAM
  1346. if strand=="+":
  1347. guideStart = pamStart+pamLen
  1348. flankSeq = seq[guideStart:guideStart+GUIDELEN]
  1349. pamSeq = seq[pamStart:pamStart+pamLen]
  1350. if pamStart-pamPlusLen >= 0:
  1351. pamPlusSeq = seq[pamStart-pamPlusLen:pamStart]
  1352. else: # strand is minus
  1353. guideStart = pamStart-GUIDELEN
  1354. flankSeq = revComp(seq[guideStart:pamStart])
  1355. pamSeq = revComp(seq[pamStart:pamStart+pamLen])
  1356. if pamStart+pamLen+pamPlusLen < len(seq):
  1357. pamPlusSeq = revComp(seq[pamStart+pamLen:pamStart+pamLen+pamPlusLen])
  1358. else: # common case: get the sequence on the left side of the PAM
  1359. if strand=="+":
  1360. guideStart = pamStart-GUIDELEN
  1361. flankSeq = seq[guideStart:pamStart]
  1362. pamSeq = seq[pamStart:pamStart+pamLen]
  1363. if pamStart+pamLen+pamPlusLen < len(seq):
  1364. pamPlusSeq = seq[pamStart+pamLen:pamStart+pamLen+pamPlusLen]
  1365. else: # strand is minus
  1366. guideStart = pamStart+pamLen
  1367. flankSeq = revComp(seq[guideStart:guideStart+GUIDELEN])
  1368. pamSeq = revComp(seq[pamStart:pamStart+pamLen])
  1369. if pamStart-pamPlusLen >= 0:
  1370. pamPlusSeq = revComp(seq[pamStart-pamPlusLen:pamStart])
  1371. if "N" in flankSeq and doFilterNs:
  1372. continue
  1373. yield "s%d%s" % (pamStart, strand), pamStart, guideStart, strand, flankSeq, pamSeq, pamPlusSeq
  1374. def makeBrowserLink(dbInfo, pos, text, title, cssClasses, ctUrl=None):
  1375. " return link to genome browser (ucsc or ensembl) at pos, with given text "
  1376. if dbInfo is None:
  1377. errAbort("Your batchID relates to a genome that is not present anymore. You will have to change the version of the site. Or contact us and send us the full URL of this page.")
  1378. if dbInfo.server.startswith("Ensembl"):
  1379. baseUrl = "www.ensembl.org"
  1380. urlLabel = "Ensembl"
  1381. # link back to archive, if possible
  1382. if dbInfo.description.startswith("Ensembl "):
  1383. ensVersion = dbInfo.description.split()[1]
  1384. if ensVersion.isdigit():
  1385. baseUrl = "e%s.ensembl.org" % ensVersion
  1386. elif dbInfo.server=="EnsemblPlants":
  1387. baseUrl = "plants.ensembl.org"
  1388. elif dbInfo.server=="EnsemblMetazoa":
  1389. baseUrl = "metazoa.ensembl.org"
  1390. elif dbInfo.server=="EnsemblProtists":
  1391. baseUrl = "protists.ensembl.org"
  1392. org = dbInfo.scientificName.replace(" ", "_")
  1393. pos = pos.replace(":+","").replace(":-","") # remove the strand
  1394. url = "http://%s/%s/Location/View?r=%s" % (baseUrl, org, pos)
  1395. elif dbInfo.server=="ucsc" or dbInfo.name.startswith("GCA_") or dbInfo.name.startswith("GCF_"):
  1396. urlLabel = "UCSC"
  1397. if len(pos)>0 and pos[0].isdigit():
  1398. pos = "chr"+pos
  1399. # remove the strand
  1400. pos = pos.replace(":+","").replace(":-","")
  1401. url = "http://genome.ucsc.edu/cgi-bin/hgTracks?db=%s&position=%s" % (dbInfo.name, pos)
  1402. if ctUrl is not None:
  1403. url+= "&hgt.customText=%s" % ctUrl
  1404. # some limited support for gbrowse
  1405. elif dbInfo.server.startswith("http://"):
  1406. urlLabel = "GBrowse"
  1407. chrom, start, end, strand = parsePos(pos)
  1408. start = start+1
  1409. url = "%s/?name=%s:%d..%d" % (dbInfo.server, chrom, start, end)
  1410. else:
  1411. chrom, start, end, strand = parsePos(pos)
  1412. if chrom is not None and chrom.startswith("NC_"):
  1413. start = start+1
  1414. url = "https://www.ncbi.nlm.nih.gov/nuccore/%s?report=graph&log$=seqview&v=%d-%d" % \
  1415. (chrom, start, end)
  1416. urlLabel = "NCBI "
  1417. else:
  1418. #return "unknown genome browser server %s, please email [email hidden]" % dbInfo.server
  1419. urlLabel = None
  1420. url = "javascript:void(0)"
  1421. classStr = ""
  1422. if len(cssClasses)!=0:
  1423. classStr = ' class="%s"' % (" ".join(cssClasses))
  1424. if title is None:
  1425. if urlLabel != None:
  1426. title = "Link to %s Genome Browser" % urlLabel
  1427. else:
  1428. title = "No Genome Browser link available yet for this organism"
  1429. return '''<a title="%s"%s target="_blank" href="%s">%s</a>''' % (title, classStr, url, text)
  1430. def highlightMismatches(guide, offTarget, pamLen):
  1431. " return a string that marks mismatches between guide and offtarget with * "
  1432. if pamLen!=0:
  1433. if pamIsFirst:
  1434. offTarget = offTarget[pamLen:]
  1435. else:
  1436. offTarget = offTarget[:-pamLen]
  1437. assert(len(guide)==len(offTarget))
  1438. s = []
  1439. for x, y in zip(guide, offTarget):
  1440. if x==y:
  1441. s.append(".")
  1442. else:
  1443. s.append("*")
  1444. return "".join(s)
  1445. def parseNewAlias(ifh):
  1446. " part of parseAlias(): IGV-compatible format: first is UCSC, all other columns are aliases "
  1447. toUcsc = {}
  1448. for line in ifh:
  1449. if line.startswith("#"):
  1450. continue
  1451. row = line.rstrip("\n").split("\t")
  1452. for i in range(1, len(row)):
  1453. toUcsc[row[i]] = row[0]
  1454. return toUcsc
  1455. def parseAlias(fname):
  1456. " parse tsv file with at least two columns, orig chrom name and new chrom name. copied from chromToUcsc script from the UCSC tools. "
  1457. logging.debug("alias file is in IGV-format")
  1458. toUcsc = {}
  1459. if fname.startswith("http://") or fname.startswith("https://"):
  1460. ifh = urlopen(fname)
  1461. if fname.endswith(".gz"):
  1462. data = gzip.GzipFile(fileobj=ifh).read().decode()
  1463. ifh = data.splitlines()
  1464. elif fname.endswith(".gz"):
  1465. ifh = gzip.open(fname, "rt")
  1466. else:
  1467. ifh = open(fname)
  1468. firstLine = True
  1469. for line in ifh:
  1470. if line.startswith("#") and firstLine:
  1471. return parseNewAlias(ifh)
  1472. if line.startswith("alias"):
  1473. continue
  1474. row = line.rstrip("\n").split("\t")
  1475. toUcsc[row[0]] = row[1]
  1476. firstLine = False
  1477. return toUcsc
  1478. chromAlias = None
  1479. def applyChromAlias(db, chrom):
  1480. " if chrom is in chromAlias, return the human-readable name "
  1481. global chromAlias
  1482. if chromAlias==-1: # == chromAlias file not present
  1483. return chrom
  1484. elif chromAlias is None:
  1485. chromAliasFname = join("genomes", db, db+".chromAlias.txt")
  1486. if not isfile(chromAliasFname):
  1487. chromAlias = -1
  1488. return chrom
  1489. else:
  1490. chromAlias = parseAlias(chromAliasFname)
  1491. return chromAlias.get(chrom, chrom)
  1492. def makeAlnStr(org, seq1, seq2, pam, mitScore, cfdScore, posStr, chromDist):
  1493. """ given two strings of equal length, return a html-formatted string of several lines
  1494. that show the two sequences and a line that highlights where they differ
  1495. """
  1496. lines = [ [], [], [] ]
  1497. last12MmCount = 0
  1498. inLinkage = False
  1499. hlSeed = False
  1500. if pamIsSpCas9(pam):
  1501. hlSeed = True
  1502. if pamIsFirst:
  1503. lines[0].append("<i>"+seq1[:len(pam)]+"</i> ")
  1504. lines[1].append("<i>"+seq2[:len(pam)]+"</i> ")
  1505. lines[2].append("".join([" "]*(len(pam)+1)))
  1506. if pamIsFirst:
  1507. guideStart = len(pam)
  1508. guideEnd = len(seq1)
  1509. else:
  1510. guideStart = 0
  1511. guideEnd = len(seq1)-len(pam)
  1512. for i in range(guideStart, guideEnd):
  1513. if hlSeed and i==10:
  1514. lines[1].append("<u>")
  1515. if seq1[i]==seq2[i]:
  1516. lines[0].append(seq1[i])
  1517. lines[1].append(seq2[i])
  1518. lines[2].append(" ")
  1519. else:
  1520. lines[0].append("<b>%s</b>" % seq1[i])
  1521. lines[1].append("<b>%s</b>" % seq2[i])
  1522. lines[2].append("*")
  1523. if i>7:
  1524. last12MmCount += 1
  1525. if hlSeed and i==guideEnd-1:
  1526. lines[1].append("</u>")
  1527. if not pamIsFirst:
  1528. lines[0].append(" <i>"+seq1[-len(pam):]+"</i>")
  1529. lines[1].append(" <i>"+seq2[-len(pam):]+"</i>")
  1530. lines = ["".join(l) for l in lines]
  1531. chrom, chromPos, strand = posStr.split(":")
  1532. chrom = applyChromAlias(org, chrom)
  1533. posStr = ":".join((chrom, chromPos, strand))
  1534. if len(posStr)>1 and posStr[0].isdigit():
  1535. posStr = "chr"+posStr
  1536. htmlText1 = "<small><pre>guide: %s<br>off-target: %s<br> %s</pre>" \
  1537. % (lines[0], lines[1], lines[2])
  1538. if pamIsCpf1(pam) or pamIsCasX(pam):
  1539. htmlText2 = "Cpf1/CasX: No off-target scores available</small>"
  1540. elif saCas9Mode:
  1541. htmlText2 = "SaCas9 Tycko Score: %s" % mitScore
  1542. else:
  1543. if cfdScore==None:
  1544. cfdStr = "Cannot calculate CFD score on non-ACTG characters"
  1545. else:
  1546. cfdStr = "%f" % cfdScore
  1547. htmlText2 = "CFD Off-target score: %s<br>MIT Off-target score: %.2f<br>Position: %s</small>" % (cfdStr, mitScore, posStr)
  1548. if chromDist!=None and org!=None:
  1549. htmlText2 += "<br><small>Distance from target: %.3f Mbp</small>" % (float(chromDist)/1000000.0)
  1550. if org.startswith("mm") or org.startswith("hg") or org.startswith("rn"):
  1551. if chromDist > 20000000:
  1552. htmlText2 += "<br><small>&gt;20Mbp = unlikely to be in linkage with target</small>"
  1553. else:
  1554. htmlText2 += "<br><small>&lt;20Mbp= likely to be in linkage with "
  1555. "target! Even if no linkage: beware of chromosomal rearrangements "
  1556. "when using this guide!</small>"
  1557. inLinkage = True
  1558. hasLast12Mm = last12MmCount>0
  1559. return htmlText1+htmlText2, hasLast12Mm, inLinkage
  1560. def parsePos(text):
  1561. """ parse a string of format chr:start-end:strand and return a 4-tuple
  1562. Strand defaults to + and end defaults to start+23
  1563. """
  1564. if text!=None and len(text)!=0 and text!="?":
  1565. fields = text.split(":")
  1566. if len(fields)==2:
  1567. chrom, posRange = fields
  1568. strand = "+"
  1569. else:
  1570. chrom, posRange, strand = fields
  1571. posRange = posRange.replace(",","")
  1572. if "-" in posRange:
  1573. start, end = posRange.split("-")
  1574. start, end = int(start), int(end)
  1575. if start > end:
  1576. start, end = end, start
  1577. strand = "-"
  1578. else:
  1579. # if the end position is not specified (as by default done by UCSC outlinks), use start+23
  1580. start = int(posRange)
  1581. end = start+23
  1582. else:
  1583. chrom, start, end, strand = "", 0, 0, "+"
  1584. return chrom, start, end, strand
  1585. def annotateOfftargets(org, countDict, guideSeq, pam, inputPos):
  1586. """ for a given guide sequence, return a list of tuples that
  1587. describes the offtargets sorted by score and a string to describe the offtargets in the
  1588. format x/y/z/w of mismatch counts
  1589. inputPos has format "chrom:start-end:strand". All 0MM matches in this range
  1590. are ignored from scoring ("ontargets")
  1591. Also return the same description for just the last 12 bp and the score
  1592. of the guide sequence (calculated using all offtargets).
  1593. """
  1594. inChrom, inStart, inEnd, inStrand = parsePos(inputPos)
  1595. count = 0
  1596. otCounts = []
  1597. posList = []
  1598. mitOtScores = []
  1599. cfdScores = []
  1600. last12MmCounts = []
  1601. ontargetDesc = ""
  1602. repCount = 0 # if repCount for a guide is !=0, then the guide should not be used. repCount is then the number
  1603. # of matches for the guide in the genome (not looking at the PAM)
  1604. # for each edit distance, get the off targets and iterate over them
  1605. foundOneOntarget = False
  1606. isSaCas9 = pamIsSaCas9(pam)
  1607. isCpf1 = pamIsCpf1(pam)
  1608. for editDist in range(0, maxMMs+1):
  1609. #print countDict,"<p>"
  1610. matches = countDict.get(editDist, [])
  1611. #print otCounts,"<p>"
  1612. last12MmOtCount = 0
  1613. # create html and score for every offtarget
  1614. otCount = 0
  1615. for chrom, start, end, otSeq, strand, segType, geneNameStr, totalAlnCount, isRep in matches:
  1616. # if repCount is > 0, then this means that the guide should not be used and we cannot
  1617. # even get any off-targets
  1618. if (totalAlnCount > MAXOCC) or (totalAlnCount > 1 and isRep):
  1619. repCount = totalAlnCount # any off-target with this condition will trigger the whole guide to be suppressed
  1620. # skip on-targets
  1621. if segType!="":
  1622. segTypeDesc = segTypeConv[segType]
  1623. geneDesc = segTypeDesc+":"+geneNameStr
  1624. geneDesc = geneDesc.replace("|", "-")
  1625. else:
  1626. geneDesc = geneNameStr
  1627. # is this not an off-target but the on-target?
  1628. # if we got a genome position, use it. Otherwise use a random off-target with 0MMs
  1629. # as the on-target ("auto-ontarget" mode)
  1630. if editDist==0 and \
  1631. repCount==0 and \
  1632. ((chrom==inChrom and start >= inStart and end <= inEnd) \
  1633. or (inChrom=='' and foundOneOntarget==False)):
  1634. foundOneOntarget = True
  1635. ontargetDesc = geneDesc
  1636. continue
  1637. otCount += 1
  1638. guideNoPam = guideSeq[:len(guideSeq)-len(pam)]
  1639. otSeqNoPam = otSeq[:len(otSeq)-len(pam)]
  1640. if len(otSeqNoPam)==19:
  1641. otSeqNoPam = "A"+otSeqNoPam # should not change the score a lot, weight0 is very low
  1642. guideNoPam = "A"+guideNoPam
  1643. if isCpf1:
  1644. # Cpf1 has no off-target scores yet
  1645. mitScore=0.0
  1646. cfdScore=0.0
  1647. elif isSaCas9:
  1648. mitScore = calcSaHitScore(guideNoPam, otSeqNoPam)
  1649. cfdScore = -1
  1650. else:
  1651. # MIT score must not include the PAM
  1652. mitScore = calcHitScore(guideNoPam, otSeqNoPam)
  1653. # this is a heuristic based on the guideSeq data where alternative
  1654. # PAMs represent only ~10% of all cleaveage events.
  1655. # We divide the MIT score by 5 to make sure that these off-targets
  1656. # are not ranked among the top but still appear in the list somewhat
  1657. if pam=="NGG" and otSeq[-2:]!="GG":
  1658. mitScore = mitScore * 0.2
  1659. # CFD score must include the PAM
  1660. cfdScore = calcCfdScore(guideSeq, otSeq)
  1661. mitOtScores.append(mitScore)
  1662. if cfdScore != -1:
  1663. cfdScores.append(cfdScore)
  1664. posStr = "%s:%d-%s:%s" % (chrom, start+1,end, strand)
  1665. if (chrom==inChrom):
  1666. dist = abs(start-inStart)
  1667. else:
  1668. dist = None
  1669. parNum = isInPar(org, chrom, start, end)
  1670. if parNum is not None:
  1671. posStr += " PAR%s" % parNum
  1672. alnHtml, hasLast12Mm, inLinkage = makeAlnStr(org, guideSeq, otSeq,
  1673. pam, mitScore, cfdScore, posStr, dist)
  1674. if not hasLast12Mm:
  1675. last12MmOtCount+=1
  1676. posList.append( (otSeq, mitScore, cfdScore, editDist, posStr, geneDesc,
  1677. alnHtml, inLinkage) )
  1678. last12MmCounts.append(str(last12MmOtCount))
  1679. # create a list of number of offtargets for this edit dist
  1680. otCounts.append( str(otCount) )
  1681. # calculate the guide scores
  1682. if pamIsCpf1(pam):
  1683. guideScore = -1
  1684. guideCfdScore = -1
  1685. else:
  1686. if repCount>0:
  1687. guideScore = 0
  1688. guideCfdScore = 0
  1689. else:
  1690. guideScore = calcMitGuideScore(sum(mitOtScores))
  1691. if doCfdFix:
  1692. guideCfdScore = calcCfdGuideScore(sum(cfdScores))
  1693. else:
  1694. guideCfdScore = calcMitGuideScore(sum(cfdScores))
  1695. # obtain the off-target info: coordinates, descriptions and off-target counts
  1696. if repCount>0:
  1697. posList = []
  1698. ontargetDesc = ""
  1699. last12DescStr = ""
  1700. otDescStr = ""
  1701. else:
  1702. otDescStr = "&thinsp;-&thinsp;".join(otCounts)
  1703. last12DescStr = "&thinsp;-&thinsp;".join(last12MmCounts)
  1704. if pamIsCpf1(pam):
  1705. # sort by edit dist if using Cfp1
  1706. posList.sort(key=operator.itemgetter(3))
  1707. else:
  1708. # sort by CFD score if we have it
  1709. posList.sort(reverse=True, key=operator.itemgetter(2))
  1710. return posList, otDescStr, guideScore, guideCfdScore, last12DescStr, \
  1711. ontargetDesc, repCount
  1712. # --- START OF SCORING ROUTINES
  1713. saGuide = None
  1714. saScorer = None
  1715. def calcSaHitScore(guideSeq, otSeq):
  1716. """
  1717. saCas9 offtarget scoring from Tycko et al, https://www.nature.com/articles/s41467-018-05391-2
  1718. see bin/src/pairwise-library-screen/
  1719. """
  1720. global saScorer
  1721. global saGuide
  1722. if guideSeq!=saGuide:
  1723. sys.path.append("bin/src/pairwise-library-screen")
  1724. import predictSingle
  1725. saGuide = guideSeq
  1726. saScorer = predictSingle.SaCas9Scorer(len(guideSeq))
  1727. # to be compatible with the MIT score, has to be in the range 0-100
  1728. # for the MIT aggregate guide specificity score
  1729. return 100.0*saScorer.calcScore(guideSeq, otSeq)
  1730. # MIT offtarget scoring, "Hsu score"
  1731. # aka Matrix "M"
  1732. hitScoreM = [0,0,0.014,0,0,0.395,0.317,0,0.389,0.079,0.445,0.508,0.613,0.851,0.732,0.828,0.615,0.804,0.685,0.583]
  1733. def calcHitScore(string1,string2):
  1734. " see 'Scores of single hits' on http://crispr.mit.edu/about "
  1735. # The Patrick Hsu weighting scheme
  1736. # S. aureus requires 21bp long guides. We fudge by using only last 20bp
  1737. matrixStart = 0
  1738. maxDist = 19
  1739. assert(string1[0].isupper())
  1740. assert(len(string1)==len(string2))
  1741. #for nmCas9 and a few others with longer guides, we limit ourselves to 20bp
  1742. if len(string1)>20:
  1743. string1 = string1[-20:]
  1744. string2 = string2[-20:]
  1745. # for 19bp guides, we fudge a little, but first pos has no weight anyways
  1746. elif len(string1)==19:
  1747. string1 = "A"+string1
  1748. string2 = "A"+string2
  1749. # for shorter guides, I'm not sure if this score makes sense anymore, we force things
  1750. elif len(string1)<19:
  1751. matrixStart = 20-len(string1)
  1752. maxDist = len(string1)-1
  1753. assert(len(string1)==len(string2))
  1754. dists = [] # distances between mismatches, for part 2
  1755. mmCount = 0 # number of mismatches, for part 3
  1756. lastMmPos = None # position of last mismatch, used to calculate distance
  1757. score1 = 1.0
  1758. for pos in range(matrixStart, len(string1)):
  1759. if string1[pos]!=string2[pos]:
  1760. mmCount+=1
  1761. if lastMmPos!=None:
  1762. dists.append(pos-lastMmPos)
  1763. score1 *= 1-hitScoreM[pos]
  1764. lastMmPos = pos
  1765. # 2nd part of the score
  1766. if mmCount<2: # special case, not shown in the paper
  1767. score2 = 1.0
  1768. else:
  1769. avgDist = sum(dists)/len(dists)
  1770. score2 = 1.0 / (((maxDist-avgDist)/float(maxDist)) * 4 + 1)
  1771. # 3rd part of the score
  1772. if mmCount==0: # special case, not shown in the paper
  1773. score3 = 1.0
  1774. else:
  1775. score3 = 1.0 / (mmCount**2)
  1776. score = score1 * score2 * score3 * 100
  1777. return score
  1778. def calcMitGuideScore(hitSum):
  1779. """ Sguide defined on http://crispr.mit.edu/about
  1780. Input is the sum of all off-target hit scores. Returns the specificity of the guide.
  1781. """
  1782. score = 100 / (100+hitSum)
  1783. score = int(round(score*100))
  1784. return score
  1785. def calcCfdGuideScore(hitSum):
  1786. " suggested by Nicholas Parkinson "
  1787. norm_score = 100.* 100.0 / (hitSum)
  1788. return norm_score
  1789. # === SOURCE CODE cfd-score-calculator.py provided by John Doench =====
  1790. # The CFD score is an improved specificity score
  1791. def get_mm_pam_scores():
  1792. """
  1793. """
  1794. import pickle
  1795. dataDir = join(dirname(__file__), 'CFD_Scoring')
  1796. mm_scores = pickle.load(open(join(dataDir, 'mismatch_score.pkl'),'rb'))
  1797. pam_scores = pickle.load(open(join(dataDir, 'pam_scores.pkl'),'rb'))
  1798. return (mm_scores,pam_scores)
  1799. #Reverse complements a given string
  1800. def revcom(s):
  1801. basecomp = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A','U':'A'}
  1802. letters = list(s[::-1])
  1803. letters = [basecomp[base] for base in letters]
  1804. return ''.join(letters)
  1805. #Calculates CFD score
  1806. def calc_cfd(wt,sg,pam):
  1807. #mm_scores,pam_scores = get_mm_pam_scores()
  1808. score = 1
  1809. sg = sg.replace('T','U')
  1810. wt = wt.replace('T','U')
  1811. s_list = list(sg)
  1812. wt_list = list(wt)
  1813. for i,sl in enumerate(s_list):
  1814. if wt_list[i] == sl:
  1815. score*=1
  1816. else:
  1817. key = 'r'+wt_list[i]+':d'+revcom(sl)+','+str(i+1)
  1818. score*= mm_scores[key]
  1819. score*=pam_scores[pam]
  1820. return (score)
  1821. mm_scores, pam_scores = None, None
  1822. def calcCfdScore(guideSeq, otSeq):
  1823. """ based on source code provided by John Doench
  1824. >>> calcCfdScore("GGGGGGGGGGGGGGGGGGGGGGG", "GGGGGGGGGGGGGGGGGAAAGGG")
  1825. 0.4635989007074176
  1826. >>> calcCfdScore("GGGGGGGGGGGGGGGGGGGGGGG", "GGGGGGGGGGGGGGGGGGGGGGG")
  1827. 1.0
  1828. >>> calcCfdScore("GGGGGGGGGGGGGGGGGGGGGGG", "aaaaGaGaGGGGGGGGGGGGGGG")
  1829. 0.5140384614450001
  1830. # mismatches: * !!
  1831. >>> calcCfdScore("ATGGTCGGACTCCCTGCCAGAGG", "ATGGTGGGACTCCCTGCCAGAGG")
  1832. 0.5
  1833. # mismatches: * ** *
  1834. >>> calcCfdScore("ATGGTCGGACTCCCTGCCAGAGG", "ATGATCCAAATCCCTGCCAGAGG")
  1835. 0.53625000020625
  1836. >>> calcCfdScore("ATGTGGAGATTGCCACCTACCGG", "ATCTGGAGATTGCCACCTACAGG")
  1837. 0.384615385
  1838. """
  1839. global mm_scores, pam_scores
  1840. if mm_scores is None:
  1841. mm_scores,pam_scores = get_mm_pam_scores()
  1842. wt = guideSeq.upper()
  1843. off = otSeq.upper()
  1844. m_wt = re.search('[^ATCG]',wt)
  1845. m_off = re.search('[^ATCG]',off)
  1846. if (m_wt is None) and (m_off is None):
  1847. pam = off[-2:]
  1848. sg = off[:20]
  1849. cfd_score = calc_cfd(wt,sg,pam)
  1850. if doCfdFix:
  1851. cfd_score = cfd_score*100.0
  1852. return cfd_score
  1853. return -1
  1854. # ==== END CFD score source provided by John Doench
  1855. # --- END OF SCORING ROUTINES
  1856. def getSizeFname(genome):
  1857. " return name of chrom.sizes file "
  1858. genomeDir = genomesDir # make local
  1859. sizeFname = "%(genomeDir)s/%(genome)s/%(genome)s.sizes" % locals()
  1860. return sizeFname
  1861. def parseChromSizes(genome):
  1862. " return chrom sizes as dict chrom -> size "
  1863. sizeFname = getSizeFname(genome)
  1864. ret = {}
  1865. for line in open(sizeFname).read().splitlines():
  1866. fields = line.split()
  1867. chrom, size = fields[:2]
  1868. ret[chrom] = int(size)
  1869. return ret
  1870. def extendAndGetSeq(db, chrom, start, end, strand, oldSeq, flank=FLANKLEN):
  1871. """ extend (start, end) by flank and get sequence for it using twoBitTwoFa.
  1872. Return None if not possible to extend.
  1873. #>>> extendAndGetSeq("hg19", "chr21", 10000000, 10000005, "+", flank=3)
  1874. #'AAGGAATGTAG'
  1875. """
  1876. assert("|" not in chrom) # we are using | to split info in BED files. | is not allowed in the fasta
  1877. chromSizes = parseChromSizes(db)
  1878. maxEnd = chromSizes[chrom]+1
  1879. start -= flank
  1880. end += flank
  1881. if start < 0 or end > maxEnd:
  1882. return None
  1883. genomeDir = genomesDir
  1884. twoBitFname = "%(genomeDir)s/%(db)s/%(db)s.2bit" % locals()
  1885. progDir = binDir
  1886. genome = db
  1887. cmd = "%(progDir)s/twoBitToFa %(genomeDir)s/%(genome)s/%(genome)s.2bit stdout -seq='%(chrom)s' -start=%(start)s -end=%(end)s" % locals()
  1888. proc = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, encoding="utf8")
  1889. seqStr = proc.stdout.read()
  1890. proc.wait()
  1891. if proc.returncode!=0:
  1892. errAbort("Could not run '%s'. Return code %s" % (cmd, str(proc.returncode)))
  1893. faFile = StringIO(seqStr)
  1894. seqs = parseFasta(faFile)
  1895. assert(len(seqs)==1)
  1896. seq = list(seqs.values())[0].upper()
  1897. if strand=="-":
  1898. seq = revComp(seq)
  1899. genomeSeq = seq[FLANKLEN:(FLANKLEN+len(oldSeq))].upper()
  1900. if oldSeq.upper() not in genomeSeq:
  1901. logging.warn("Input sequence has SNPs compared to genome, not returning extended seq:")
  1902. logging.warn("- Input sequence: %s" % oldSeq)
  1903. logging.warn("- Genome sequence: %s" % genomeSeq)
  1904. logging.warn("- Diff String : %s" % highlightMismatches(oldSeq, genomeSeq, 0))
  1905. return None
  1906. # ? make sure that user annotations, like added Ns, are retained in the long sequence
  1907. #fixedSeq = seq[:100]+oldSeq+seq[-100:]
  1908. #assert(len(fixedSeq)==len(seq))
  1909. return seq
  1910. def getExtSeq(seq, start, end, strand, extUpstream, extDownstream, extSeq=None, extFlank=FLANKLEN):
  1911. """ extend (start,end) by extUpstream and extDownstream and return the subsequence
  1912. at this position in seq.
  1913. Return None if there is not enough space to extend (start, end).
  1914. extSeq is a sequence with extFlank additional flanking bases on each side. It can be provided
  1915. optionally and is used if needed to return a subseq.
  1916. Careful: returned sequence might contain lowercase letters.
  1917. >>> getExtSeq("AACCTTGG", 2, 4, "+", 2, 4)
  1918. 'AACCTTGG'
  1919. >>> getExtSeq("CCAACCTTGGCC", 4, 6, "-", 2, 3)
  1920. 'AAGGTTG'
  1921. >>> getExtSeq("AA", 0, 2, "+", 2, 3)
  1922. >>> getExtSeq("AA", 0, 2, "+", 2, 3, extSeq="CAGAATGA", extFlank=3)
  1923. 'AGAATGA'
  1924. >>> getExtSeq("AA", 0, 2, "-", 2, 3, extSeq="CAGAATGA", extFlank=3)
  1925. 'CATTCTG'
  1926. """
  1927. assert(start>=0)
  1928. assert(end<=len(seq))
  1929. # extend
  1930. if strand=="+":
  1931. extStart, extEnd = start-extUpstream, end+extDownstream
  1932. else:
  1933. extStart, extEnd = start-extDownstream, end+extUpstream
  1934. # check for out of bounds and get seq
  1935. if extStart >= 0 and extEnd <= len(seq):
  1936. logging.debug("using input seq, pos %d-%d" % (extStart, extEnd))
  1937. subSeq = seq[extStart:extEnd]
  1938. else:
  1939. if extSeq==None:
  1940. return None
  1941. # lift to extSeq coords and get seq
  1942. extStart += extFlank
  1943. extEnd += extFlank
  1944. assert(extStart >= 0)
  1945. assert(extEnd <= len(extSeq))
  1946. subSeq = extSeq[extStart:extEnd]
  1947. logging.debug("using extended seq, pos %d-%d" % (extStart, extEnd))
  1948. if strand=="-":
  1949. logging.debug("revcomp'ing result")
  1950. subSeq = revComp(subSeq)
  1951. # check that the extended sequence really contains the whole input seq
  1952. # e.g. when user has added nucleotides to a otherwise matching sequence
  1953. #if seq.upper() not in subSeq.upper():
  1954. #debug("seq is not in extSeq")
  1955. #subSeq = None
  1956. logging.debug("Got -%d/+%d-extended seq for (%d, %d, %s) = %s. Result: %s." %
  1957. (extUpstream, extDownstream, start, end, strand, seq[start:end], subSeq))
  1958. return subSeq
  1959. def pamStartToGuideRange(startPos, strand, pamLen):
  1960. """ given a PAM start position and its strand, return the (start,end) of the guide.
  1961. Coords can be negative or exceed the length of the input sequence.
  1962. """
  1963. if not pamIsFirst:
  1964. if strand=="+":
  1965. return (startPos-GUIDELEN, startPos)
  1966. else: # strand is minus
  1967. return (startPos+pamLen, startPos+pamLen+GUIDELEN)
  1968. else:
  1969. if strand=="+":
  1970. return (startPos+pamLen, startPos+pamLen+GUIDELEN)
  1971. else: # strand is minus
  1972. return (startPos-GUIDELEN, startPos)
  1973. def htmlHelp(text):
  1974. " show help text with tooltip or modal dialog "
  1975. className = "tooltipster"
  1976. if "href" in text:
  1977. className = "tooltipsterInteract"
  1978. print('''<img style="padding-bottom: 3px; height:1.1em; width:1.0em" src="%simage/info-small.png" class="help %s" title="%s" />''' % (HTMLPREFIX, className, text))
  1979. def htmlWarn(text):
  1980. " show help text with tooltip "
  1981. print('''<img style="height:0.9em; width:0.8em; padding-bottom: 2px" src="%simage/warning-32.png" class="help tooltipster" title="%s" />''' % (HTMLPREFIX, text))
  1982. def readRestrEnzymes():
  1983. """ parse restrSites.txt and
  1984. return as dict length -> list of (name, suppliers, seq) """
  1985. fname = "restrSites.txt"
  1986. enzList = {}
  1987. for line in open(join(baseDir, fname)):
  1988. if line.startswith("#"):
  1989. continue
  1990. seq, name, suppliers = line.rstrip("\n").rstrip("\r").split("\t")
  1991. suppliers = tuple(suppliers.split(","))
  1992. enzList.setdefault(len(seq), []).append( (name, suppliers, seq) )
  1993. return enzList
  1994. def patMatch(seq, pat, notDegPos=None):
  1995. """ return true if pat matches seq, both have to be same length
  1996. do not match degenerate codes at position notDegPos (0-based)
  1997. """
  1998. assert(len(seq)==len(pat))
  1999. for x in range(0, len(pat)):
  2000. patChar = pat[x]
  2001. nuc = seq[x]
  2002. assert(patChar in "MKYRACTGNWSDVBH")
  2003. assert(nuc in "MKYRACTGNWSDX")
  2004. if notDegPos!=None and x==notDegPos and patChar!=nuc:
  2005. return False
  2006. if nuc=="X":
  2007. return False
  2008. if patChar=="N":
  2009. continue
  2010. if patChar=="H" and nuc in "ACT":
  2011. continue
  2012. if patChar=="D" and nuc in "AGT":
  2013. continue
  2014. if patChar=="B" and nuc in "CGT":
  2015. continue
  2016. if patChar=="V" and nuc in "ACG":
  2017. continue
  2018. if patChar=="W" and nuc in "AT":
  2019. continue
  2020. if patChar=="S" and nuc in "GC":
  2021. continue
  2022. if patChar=="M" and nuc in "AC":
  2023. continue
  2024. if patChar=="K" and nuc in "TG":
  2025. continue
  2026. if patChar=="R" and nuc in "AG":
  2027. continue
  2028. if patChar=="Y" and nuc in "CT":
  2029. continue
  2030. if patChar!=nuc:
  2031. return False
  2032. return True
  2033. def findSite(seq, restrSite):
  2034. """ return the positions where restrSite matches seq
  2035. seq can be longer than restrSite
  2036. Do not allow degenerate characters to match at position len(restrSite) in seq
  2037. """
  2038. posList = []
  2039. for i in range(0, len(seq)-len(restrSite)+1):
  2040. subseq = seq[i:i+len(restrSite)]
  2041. # JP does not want any potential site to be suppressed
  2042. #if i<len(restrSite):
  2043. #isMatch = patMatch(subseq, restrSite, len(restrSite)-i-1)
  2044. #else:
  2045. #isMatch = patMatch(subseq, restrSite)
  2046. isMatch = patMatch(subseq, restrSite)
  2047. if isMatch:
  2048. posList.append( (i, i+len(restrSite)) )
  2049. return posList
  2050. def matchRestrEnz(allEnzymes, guideSeq, pamSeq, pamPlusSeq, pamPat):
  2051. """ return list of enzymes that overlap the -3 position in guideSeq
  2052. returns dict (name, pattern, suppliers) -> list of matching positions
  2053. """
  2054. matches = defaultdict(set)
  2055. if pamPlusSeq is None:
  2056. pamPlusSeq = "XXXXX" # make sure that we never match a restriction site outside the seq boundaries
  2057. fullSeq = concatGuideAndPam(guideSeq, pamSeq, pamPlusSeq)
  2058. #print guideSeq, pamSeq, pamPlusSeq, fullSeq, "<br>"
  2059. for siteLen, sites in allEnzymes.items():
  2060. if pamIsCpf1(pamPat):
  2061. # most modified position: 4nt from the end
  2062. # see http://www.nature.com/nbt/journal/v34/n8/full/nbt.3620.html
  2063. # Figure 1
  2064. startSeq = len(fullSeq)-4-pamPlusLen-(siteLen)+1
  2065. else:
  2066. # most modified position for Cas9: 3bp from the end
  2067. startSeq = len(fullSeq)-len(pamSeq)-3-pamPlusLen-(siteLen)+1
  2068. seq = fullSeq[startSeq:].upper()
  2069. for name, suppliers, restrSite in sites:
  2070. posList = findSite(seq, restrSite)
  2071. if len(posList)!=0:
  2072. liftOffset = startSeq
  2073. posList = [(liftOffset+x, liftOffset+y) for x,y in posList]
  2074. matches.setdefault((name, restrSite, suppliers), set()).update(posList)
  2075. return matches
  2076. def mergeGuideInfo(seq, startDict, pamPat, otMatches, inputPos, effScores, sortBy=None, org=None):
  2077. """
  2078. merges guide information from the sequence, the efficiency scores and the off-targets.
  2079. creates rows with too many fields. needs refactoring.
  2080. for each pam in startDict, retrieve the guide sequence next to it and score it
  2081. sortBy can be "effScore", "mhScore", "oofScore" or "pos"
  2082. """
  2083. allEnzymes = readRestrEnzymes()
  2084. guideData = []
  2085. guideScores = {}
  2086. hasNotFound = False
  2087. pamIdToSeq = {}
  2088. pamSeqs = list(flankSeqIter(seq.upper(), startDict, len(pamPat), True))
  2089. for pamId, pamStart, guideStart, strand, guideSeq, pamSeq, pamPlusSeq in pamSeqs:
  2090. # matches in genome
  2091. # one desc in last column per OT seq
  2092. if pamId in otMatches:
  2093. pamMatches = otMatches[pamId]
  2094. guideSeqFull = concatGuideAndPam(guideSeq, pamSeq)
  2095. mutEnzymes = matchRestrEnz(allEnzymes, guideSeq, pamSeq, pamPlusSeq, pamPat)
  2096. posList, otDesc, guideScore, guideCfdScore, last12Desc, ontargetDesc, \
  2097. repCount = \
  2098. annotateOfftargets(org, pamMatches, guideSeqFull, pamPat, inputPos)
  2099. if repCount!=0:
  2100. guideScore = 0
  2101. guideCfdScore = 0
  2102. # no off-targets found?
  2103. else:
  2104. posList, otDesc, guideScore = None, "Not found", -1
  2105. guideCfdScore = -1
  2106. last12Desc = ""
  2107. hasNotFound = True
  2108. mutEnzymes = []
  2109. ontargetDesc = ""
  2110. repCount = 0
  2111. seq34Mer = None
  2112. guideRow = [guideScore, guideCfdScore, effScores.get(pamId, {}), pamStart, guideStart, strand, pamId, guideSeq, pamSeq, posList, otDesc, last12Desc, mutEnzymes, ontargetDesc, repCount]
  2113. guideData.append( guideRow )
  2114. guideScores[pamId] = guideScore
  2115. pamIdToSeq[pamId] = guideSeq
  2116. if sortBy == "pos":
  2117. sortFunc = (lambda row: row[3])
  2118. reverse = False
  2119. elif sortBy == "offCount":
  2120. sortFunc = (lambda row: len(row[9]))
  2121. reverse = False
  2122. elif sortBy == "cfdSpec":
  2123. sortFunc = operator.itemgetter(1)
  2124. reverse = True
  2125. elif sortBy == "spec" or sortBy is None:
  2126. sortFunc = (lambda row: row[0])
  2127. reverse = True
  2128. elif sortBy is not None and not sortBy.endswith("pec"):
  2129. sortFunc = (lambda row: row[2].get(sortBy, 0))
  2130. reverse = True
  2131. else:
  2132. errAbort("Unknown sortBy value. This is a bug. Please contact us.")
  2133. guideData.sort(reverse=reverse, key=sortFunc)
  2134. return guideData, guideScores, hasNotFound, pamIdToSeq
  2135. def printDownloadTableLinks(batchId, addTsv=False):
  2136. print('<div id="downloads" style="text-align:left">')
  2137. print("Download as Excel tables: ", end=' ')
  2138. print('<a href="crispor.py?batchId=%s&download=guides&format=xls">Guides</a>&nbsp;/&nbsp;' % batchId, end=' ')
  2139. if not pamIsFirst and not saCas9Mode:
  2140. print('<a href="crispor.py?batchId=%s&showAllScores=1&download=guides&format=xls">Guides, all scores</a>&nbsp;/&nbsp;' % batchId, end=' ')
  2141. print('<a href="crispor.py?batchId=%s&download=offtargets&format=xls">Off-targets</a>&nbsp;/&nbsp;' % batchId, end=' ')
  2142. print(('<a href="crispor.py?batchId=%s&satMut=1">Saturating mutagenesis assistant</a><br>' % batchId))
  2143. #print "<small>Plasmid Editor: ",
  2144. #print '<a href="crispor.py?batchId=%s&download=genbank">Guides</a></small>' % batchId,
  2145. if addTsv:
  2146. print("<small>Tab-sep format: ", end=' ')
  2147. print('<a href="crispor.py?batchId=%s&download=guides&format=tsv">Guides</a>&nbsp;/&nbsp;' % batchId, end=' ')
  2148. print('<a href="crispor.py?batchId=%s&download=offtargets&format=tsv">Off-targets</a></small>' % batchId, end=' ')
  2149. print('</div>')
  2150. def hasGeneModels(org):
  2151. " return true if this organism has gene model information "
  2152. geneFname = join(genomesDir, org, org+".segments.bed")
  2153. return isfile(geneFname)
  2154. def printTableHead(pam, batchId, chrom, org, varHtmls, showColumns):
  2155. " print guide score table description and columns "
  2156. # one row per guide sequence
  2157. if not pamIsCpf1(pam):
  2158. print('''<div class='substep'>Ranked by default from highest to lowest specificity score (<a target='_blank' href='http://dx.doi.org/10.1038/nbt.2647'>Hsu et al., Nat Biot 2013</a>). Click on a column title to rank by a score.<br>''')
  2159. #print("""<b>Our recommendation:</b> Use Fusi for in-vivo (U6) transcribed guides, Moreno-Mateos for in-vitro (T7) guides injected into Zebrafish/Mouse oocytes.<br>""")
  2160. print('''If you use this website, please cite our <a href="https://academic.oup.com/nar/article/46/W1/W242/4995687">paper in NAR 2018</a>.''')
  2161. print("Too much information? Look at the <a target=_blank href='manual/'>CRISPOR manual</a>.<p>")
  2162. print('</div>')
  2163. printDownloadTableLinks(batchId)
  2164. print("""
  2165. <script type="text/javascript">
  2166. function allRows() {
  2167. $("guideRow").show();
  2168. }
  2169. //function copySeq() {
  2170. //var c = new ClipboardJS('#seqAsText');
  2171. //var copyText = document.getElementById("seqAsText");
  2172. //var selRes = copyText.select();
  2173. //var val = copyText.value;
  2174. //var res = document.execCommand("copy");
  2175. //alert("The input sequence is now in your clipboard. You can paste it into other programs.");
  2176. //}
  2177. $(document).ready( function() {
  2178. //#$('#copyLink').click( copySeq );
  2179. var clipboard = new ClipboardJS('#copyLink');
  2180. clipboard.on('success', function(e) {
  2181. alert("The input sequence is now in your clipboard. You can paste it into other programs.");
  2182. console.info('Action:', e.action);
  2183. console.info('Text:', e.text);
  2184. console.info('Trigger:', e.trigger);
  2185. e.clearSelection();
  2186. });
  2187. $('d').mouseenter( onEditHover );
  2188. $('d').mouseleave ( onEditOut );
  2189. });
  2190. function onEditOut() {
  2191. $('#editHover').hide();
  2192. }
  2193. function colorChar(str, pos) {
  2194. /* put a span-color tag around the char at pos in str and return result */
  2195. var prefix = str.substring(0, pos);
  2196. var hlChar = str[pos];
  2197. var suffix = str.substring(pos+1);
  2198. return prefix+"<mut>"+hlChar+"</mut>"+suffix;
  2199. }
  2200. function onEditHover(ev) {
  2201. /* user hovers over an edit letter */
  2202. ev.preventDefault();
  2203. var oldEl = document.getElementById("editHover");
  2204. if (oldEl)
  2205. oldEl.remove();
  2206. console.log(ev.target);
  2207. const boundBox = ev.target.getBoundingClientRect();
  2208. var x = boundBox.left;
  2209. var y = boundBox.top;
  2210. y += 14;
  2211. var div = document.createElement('div');
  2212. div.id = "editHover";
  2213. div.style.width="400px";
  2214. div.style.height="200px";
  2215. div.style.border="1px solid black";
  2216. div.style.padding="10px";
  2217. div.style.position="fixed";
  2218. div.style.backgroundColor="white";
  2219. div.style.left=x+"px";
  2220. div.style.top=y+"px";
  2221. var pos = parseInt(this.getAttribute("pos"));
  2222. var nucl = this.textContent;
  2223. if (nucl.toUpperCase()==="T")
  2224. origNucl = "C";
  2225. else
  2226. origNucl = "G";
  2227. var htmls=[];
  2228. htmls.push("The following guides can mutate "+origNucl+" to "+nucl+" at position "+pos+":<br>");
  2229. htmls.push("<table class='editTable'>");
  2230. htmls.push("<tr><th>Guide ID</th><th>Guide Sequence</th><th>Komor score</th><th>Spec. Score</th></tr>");
  2231. var guides = editData[pos][nucl];
  2232. guides.sort( function (a, b) { a[4] - b[4] } ); // sort by komor score
  2233. for (var i=0; i<guides.length; i++) {
  2234. guide = guides[i];
  2235. pamId = guide[0];
  2236. guideSeq = guide[1];
  2237. pam = guide[2];
  2238. mutPos = guide[3];
  2239. beScore = guide[4];
  2240. specScore = guide[5];
  2241. htmls.push("<tr>");
  2242. htmls.push("<td>"+pamId+"</td>");
  2243. htmls.push("<td><tt>"+colorChar(guideSeq, mutPos)+" "+pam+"</tt></td>");
  2244. htmls.push("<td>"+beScore.toFixed(2)+"</td>");
  2245. htmls.push("<td>"+specScore+"</td>");
  2246. htmls.push("</tr>");
  2247. }
  2248. htmls.push("</table>");
  2249. $(div).append(htmls.join(""));
  2250. document.body.appendChild(div);
  2251. }
  2252. function onlyWith(doPrefix) {
  2253. /* show only guide rows and guide sequence viewer features that start with a prefix */
  2254. if ($("#onlyWith"+doPrefix+"Box").prop("checked"))
  2255. {
  2256. $(".prefixBox").prop("checked", false);
  2257. $("#onlyWith"+doPrefix+"Box").prop("checked", true);
  2258. //$(".guideRow").show();
  2259. $(".guideRow").css("visibility", "visible");
  2260. $(".guideRowNoPrefix"+doPrefix).hide();
  2261. // special handling for sequence viewer: hide() would destroy the layout there
  2262. $(".guideRowNoPrefix"+doPrefix+"Seq").css("visibility", "hidden");
  2263. }
  2264. else
  2265. {
  2266. $(".prefixBox").prop("checked", false);
  2267. $(".guideRow").show();
  2268. $(".guideRowNoPrefix"+doPrefix+"Seq").css("visibility", "visible");
  2269. }
  2270. }
  2271. function displayClass(className, dispVal) {
  2272. /* hide in a loop, works around Safari stack size limits that crash jquery functions */
  2273. var els = document.getElementsByClassName(className);
  2274. for (var el of els) {
  2275. el.style.display = dispVal;
  2276. }
  2277. }
  2278. function onlyExons() {
  2279. /* show only off-targets in exons */
  2280. if ($("#onlyExonBox").prop("checked")) {
  2281. $(".otMore").show();
  2282. $(".otMoreLink").hide();
  2283. $(".otLessLink").hide();
  2284. displayClass("notExon", "none");
  2285. }
  2286. else {
  2287. if ($("#onlySameChromBox").prop("checked")) {
  2288. $(".notExon:not(.diffChrom)").show();
  2289. }
  2290. else {
  2291. displayClass("notExon", "block");
  2292. $(".otMoreLink").show();
  2293. $(".otMore").hide();
  2294. }
  2295. }
  2296. }
  2297. function onlySameChrom() {
  2298. if ($("#onlySameChromBox").prop("checked"))
  2299. {
  2300. $(".otMore").show();
  2301. $(".otMoreLink").hide();
  2302. $(".otLessLink").hide();
  2303. $(".diffChrom").hide();
  2304. }
  2305. else {
  2306. if ($("#onlyExonBox").prop("checked")) {
  2307. $(".diffChrom:not(.notExon)").show();
  2308. }
  2309. else {
  2310. $(".diffChrom").show();
  2311. $(".otMoreLink").show();
  2312. $(".otMore").hide();
  2313. }
  2314. }
  2315. }
  2316. function showAllOts(classId) {
  2317. $("#"+classId).show();
  2318. $("#"+classId+"MoreLink").hide();
  2319. $("#"+classId+"LessLink").show();
  2320. }
  2321. function showLessOts(classId) {
  2322. $("#"+classId).hide();
  2323. $("#"+classId+"MoreLink").show();
  2324. $("#"+classId+"LessLink").hide();
  2325. }
  2326. </script>
  2327. """)
  2328. print('<table id="otTable" style="background:white;table-layout:fixed; overflow:scroll; width:100%">')
  2329. print('<thead>')
  2330. print('<tr style="border-bottom:none; border-left:5px solid black; background-color:#F0F0F0">')
  2331. print('<th style="width:80px; border-bottom:none"><a href="crispor.py?batchId=%s&sortBy=pos" class="tooltipster" title="Click to sort the table by the position of the PAM site">Position/<br>Strand</a>' % batchId)
  2332. htmlHelp("You can click on the links in this column to highlight the <br>PAM site in the sequence viewer at the top of the page.")
  2333. print('</th>')
  2334. print('<th style="width:235px; border-bottom:none">Guide Sequence + <i>PAM</i><br>')
  2335. print ('+ Restriction Enzymes')
  2336. htmlHelp("Restriction enzymes can be very useful for screening mutations induced by the guide RNA using PCR and Restrictrion frament length polymorphism (RFLP).<br>Enzyme sites shown here overlap the main cleavage site 3bp 5' to the PAM.<br>Digestion of the PCR product with these enzymes will not cut the product if the genome was mutated by Cas9. This is a lot easier than screening with the T7 assay, Surveyor or sequencing.")
  2337. print('<br>')
  2338. if varHtmls is not None:
  2339. print(' + Variants')
  2340. htmlHelp("Variants that overlap the guide sequence are shown. You can change the variant database with the drop-down box above the sequence viewer at the top of the page.")
  2341. print('<br>')
  2342. print('''<small>''')
  2343. print('''<input type="checkbox" class="prefixBox" id="onlyWithGBox" onchange="onlyWith('G')">Only G-''')
  2344. print('''<input type="checkbox" class="prefixBox" id="onlyWithGGBox" onchange="onlyWith('GG')">Only GG-''')
  2345. print('''<input type="checkbox" class="prefixBox" id="onlyWithABox" onchange="onlyWith('A')">Only A-''')
  2346. htmlHelp("The three checkboxes allow you to show only guides that start with GG-, G- or A-. While we recommend prefixing a 20bp guide with G for U6 expression with spCas9, some protocols recommend using only guides with a G- prefix for U6 and A- for U3.")
  2347. print('''</small>''')
  2348. if not pamIsCpf1(pam):
  2349. print('<th style="width:80px; border-bottom:none"><a href="crispor.py?batchId=%s&sortBy=spec" class="tooltipster" title="Click to sort the table by specificity score. Hover over the (i) bubble on the right to get more information about the specificity score.">MIT Specificity Score</a>' % batchId)
  2350. if pamIsSaCas9(pam):
  2351. htmlHelp("The higher the specificity score, the lower are off-target effects in the genome.<br>This specificity score has been adapted for SaCas9 and based on the off-target scores shown on mouse-over. The algorithm was provided by Josh Tycko. Like the MIT score for spCas9, it is aggregated from all off-target scores and ranges 0-100. See <a href='https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6063963/'>Tycko et al. Nat Comm 2018</a> for details.")
  2352. else:
  2353. htmlHelp("The higher the specificity score, the lower are off-target effects in the genome.<br>The specificity score ranges from 0-100 and measures the uniqueness of a guide in the genome. See <a href='http://dx.doi.org/10.1038/nbt.2647'>Hsu et al. Nat Biotech 2013</a>. We recommend values &gt;50, where possible. See <a target=_blank href='manual/#offs'>the CRISPOR manual</a>")
  2354. print("</th>")
  2355. if "cfdGuideScore" in showColumns:
  2356. print('<th style="width:60px; border-bottom:none"><a href="crispor.py?batchId=%s&sortBy=cfdSpec" class="tooltipster" title="Click to sort the table by CFD specificity score">CFD Spec. score</a>' % batchId)
  2357. htmlHelp("The CFD specificity score, inspired like guidescan.com, behaves like the MIT specificity score, but it is based on the more accurate CFD off-target model, from <a href='http://www.nature.com/nbt/journal/v34/n2/full/nbt.3437.html'>Doench 2016</a>, which is also used by Crispor to rank the off-targets. The CFD specificity score correlates better than the MIT score with the total off-target cleavage fraction of a guide, see <a target=_blank href='https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6731277/'>Tycko et al, Nat Comm 2019</a> and also the <a target=_blank href='/manual/#faq'>CRISPOR manual</a>.")
  2358. print("</th>")
  2359. if len(scoreNames)==2 or pamIsCpf1(pam) or pamIsSaCas9(pam):
  2360. print('<th style="width:150px; height:100px; border-bottom:none" colspan="%d">Predicted Efficiency' % (len(scoreNames)))
  2361. else:
  2362. print('<th style="width:270px; border-bottom:none" colspan="%d">Predicted Efficiency' % (len(scoreNames))) # -1 because proxGc is in scoreNames but has no column
  2363. htmlHelp("The higher the efficiency score, the more likely is cleavage at this position. For details on the scores, mouseover their titles below.<br>Note that these predictions are not very accurate, they merely enrich for more efficient guides by a factor of 2-3 so you have to test a few guides to see the effect. <a target=_blank href='manual/#onEff'>Read the CRISPOR manual</a>")
  2364. if not pamIsCpf1(pam) and not pamIsSaCas9(pam):
  2365. if cgiParams.get("showAllScores", "0")=="0":
  2366. print(("""<br><a style="font-size:12px" href="%s" class="tooltipsterInteract" title="By default, only the two most relevant scores are shown, based on our study <a href='http://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-1012-2'>Haeussler et al. 2016</a>. Click this link to show all efficiency scores.">Show all scores</a>""" % cgiGetSelfUrl({"showAllScores":"1"}, anchor="otTable")))
  2367. scoreDescs["crisprScan"][0] = "Mor.-Mateos"
  2368. else:
  2369. print(("""<br><a style="font-size:12px" href="%s" class="tooltipsterInteract" title="Show only the two main scores">Show main scores</a>""" % cgiGetSelfUrl({"showAllScores":None}, anchor="otTable")))
  2370. print('</th>')
  2371. mhColName="Outcome"
  2372. if not baseEditor:
  2373. if len(mutScoreNames)<=1:
  2374. mhColName = ""
  2375. oofWidth=45
  2376. #oofDesc = "Click on score to show micro-homology"
  2377. #oofDesc = ""
  2378. else:
  2379. oofWidth=67
  2380. #oofDesc = "Click score for details"
  2381. colSpan = len(mutScoreNames)
  2382. print('<th colspan=%d style="width:%dpx; border-bottom:none"><a href="crispor.py?batchId=%s&sortBy=oof" class="tooltipster" title="Prediction of the DNA sequence after strand break repair. Click to sort the table by frameshift/out-of-frame scores. Hover over the score names to show information about a particular score. Click a score number to see the predicted indel pattern around the guide.">%s</a>' % (colSpan, oofWidth, batchId, mhColName))
  2383. #htmlHelp(scoreDescs["oof"][1])
  2384. #print "<small>%s</small>" % oofDesc
  2385. print('</th>')
  2386. print('<th style="width:117px; border-bottom:none"><a href="crispor.py?batchId=%s&sortBy=offCount" class="tooltipster" title="Click to sort the table by number of off-targets">Off-targets for <br>0-1-2-3-4 mismatches<br></a><span style="color:grey">+ next to PAM </span>' % (batchId))
  2387. altPamsHelp = [pam]
  2388. if pam in offtargetPams:
  2389. altPamsHelp.extend(offtargetPams[pam])
  2390. htmlHelp("For each number of mismatches, the number of off-targets is indicated.<br>Example: 1-3-20-50-60 means 1 off-target with 0 mismatches, 3 off-targets with 1 mismatch, <br>20 off-targets with 2 mismatches, etc.<br>The CRISPOR website only searches up to four mismatches (use the command line version for 5 or 6). Off-targets are considered if they are flanked by one of these motifs: %s .<br>Shown in grey are the off-targets that have no mismatches in the 12 bp adjacent to the PAM. These are the most likely off-targets." % (", ".join(altPamsHelp)))
  2391. print("</th>")
  2392. print('<th style="width:*; border-bottom:none">Genome Browser links to matches sorted by CFD off-target score')
  2393. htmlHelp("For each off-target the number of mismatches is indicated and linked to a genome browser. <br>Matches are ranked by CFD off-target score (see Doench 2016 et al) from most to least likely.<br>Matches can be filtered to show only off-targets in exons or on the same chromosome as the input sequence.<br>On most organisms, you can click the links below to open a window with a genome browser at this position.")
  2394. print('<br><small>')
  2395. print('<input type="hidden" name="batchId" value="%s">' % batchId)
  2396. if hasGeneModels(org):
  2397. print('''<input type="checkbox" id="onlyExonBox" onchange="onlyExons()">exons only''')
  2398. else:
  2399. print('<small title="When this genome was loaded into CRISPOR, gene models were not available. Contact us if you want to filter for off-targets in exons and think that a gene models are now available for this genome." style="color:grey">No exons.</small>')
  2400. if chrom!="":
  2401. if chrom[0].isdigit():
  2402. chrom = "chrom "+chrom
  2403. print('''<input type="checkbox" id="onlySameChromBox" onchange="onlySameChrom()">%s only''' % chrom)
  2404. else:
  2405. print('<small style="color:grey">&nbsp;No match, no chrom filter</small>')
  2406. print("</small>")
  2407. print("</th>")
  2408. print("</tr>")
  2409. # subheaders
  2410. print('<tr style="border-top:none; border-left: solid black 5px; background-color:#F0F0F0">')
  2411. print('<th style="border-top:none"></th>')
  2412. print('<th style="border-top:none"></th>')
  2413. if "cfdGuideScore" in showColumns:
  2414. print('<th style="border-top:none"></th>')
  2415. if not pamIsCpf1(pam):
  2416. print('<th style="border-top:none"></th>')
  2417. for scoreName in scoreNames:
  2418. if scoreName in ["oof", "proxGc"] or "oof" in scoreName:
  2419. continue
  2420. scoreLabel, scoreDesc = scoreDescs[scoreName]
  2421. print('<th style="width: 10px; border: none; border-top:none; border-right: none" class="rotate"><div><span><a title="%s" class="tooltipsterInteract" href="crispor.py?batchId=%s&sortBy=%s">%s</a></span></div></th>' % (scoreDesc, batchId, scoreName, scoreLabel))
  2422. if "proxGc" in scoreNames:
  2423. # the ProxGC score comes next
  2424. print('''<th style="border: none; border-top:none; border-right: none; border-left:none" class="rotate">''')
  2425. print('''<div><span style="border-bottom:none">''')
  2426. print('''<a title="This column shows two heuristics based on observations rather than computational models: <a href='http://www.cell.com/cell-reports/abstract/S2211-1247%2814%2900827-4'>Ren et al</a> 2014 obtained the highest cleavage in Drosophila when the final 6bp contained &gt;= 4 GCs, based on data from 39 guides. <a href='http://www.genetics.org/content/early/2015/02/18/genetics.115.175166.abstract'>Farboud et al.</a> obtained the highest cleavage in C. elegans for the 10 guides that ended with -GG, out of the 50 guides they tested.<br>The column contains + if the final GC count is &gt;= 4 and GG if the guide ends with GG." href="crispor.py?batchId=%s&sortBy=finalGc6" class="tooltipsterInteract">Prox GC</span></div></th>''' % (batchId))
  2427. # these are empty cells to fill up the row and avoid white space
  2428. for scoreName in mutScoreNames:
  2429. scoreLabel, scoreDesc = scoreDescs[scoreName]
  2430. print('<th style="width: 10px; border-top:none; border-right: none" class="rotate"><div><span><a title="%s" class="tooltipsterInteract" href="crispor.py?batchId=%s&sortBy=%s">%s</a></span></div></th>' % (scoreDesc, batchId, scoreName, scoreLabel))
  2431. print('<th style="border-top:none"></th>')
  2432. print('<th style="border-top:none"></th>')
  2433. print("</tr>")
  2434. print('</thead>')
  2435. def scoreToColor(guideScore):
  2436. if guideScore is None:
  2437. color = ("#000000", "black")
  2438. elif guideScore > 50:
  2439. color = ("#32cd32", "green")
  2440. elif guideScore > 20:
  2441. color = ("#ffff00", "yellow")
  2442. elif guideScore==-1:
  2443. color = ("#000000", "black")
  2444. else:
  2445. color = ("#aa0114", "red")
  2446. return color
  2447. def hexToRgb(hexCode):
  2448. " convert hex color to RGB in UCSC format, https://stackoverflow.com/questions/29643352/converting-hex-to-rgb-value-in-python "
  2449. hexCode = hexCode.lstrip("#")
  2450. return ",".join(tuple(str(int(hexCode[i:i+2], 16)) for i in (0, 2 ,4)))
  2451. def makeOtBrowserLinks(otData, chrom, dbInfo, pamId):
  2452. " return a list with the html texts of the offtarget links "
  2453. links = []
  2454. i = 0
  2455. for otSeq, score, cfdScore, editDist, pos, gene, alnHtml, inLinkage in otData:
  2456. cssClasses = ["tooltipster"]
  2457. if not gene.startswith("exon:"):
  2458. cssClasses.append("notExon")
  2459. if pos.split(":")[0]!=chrom:
  2460. cssClasses.append("diffChrom")
  2461. if inLinkage:
  2462. cssClasses.append("inLinkage")
  2463. classStr = ""
  2464. if len(cssClasses)!=0:
  2465. classStr = ' class="%s"' % " ".join(cssClasses)
  2466. link = makeBrowserLink(dbInfo, pos, gene, alnHtml, ["tooltipster"])
  2467. editDist = str(editDist)
  2468. links.append( '''<div%(classStr)s>%(editDist)s:%(link)s</div>''' % locals() )
  2469. return links
  2470. def filterOts(otDatas, minScore):
  2471. " remove all offtargets with score < minScore "
  2472. newList = []
  2473. for otData in otDatas:
  2474. score = otData[1]
  2475. if score > minScore:
  2476. newList.append(otData)
  2477. return newList
  2478. def findOtCutoff(otData):
  2479. " try cutoffs 0.5, 1.0, 2.0, 3.0 until not more than 20 offtargets left "
  2480. for cutoff in [0.3, 0.5, 1.0, 2.0, 3.0, 10.0, 99.9]:
  2481. otData = filterOts(otData, cutoff)
  2482. if len(otData)<=30:
  2483. return otData, cutoff
  2484. if len(otData)>30:
  2485. return otData[:30], None
  2486. return otData, 1000
  2487. def printNote(s):
  2488. print('<div style="text-align:left; background-color: aliceblue; padding:5px; border: 1px solid black"><strong>Note:</strong>')
  2489. print(s)
  2490. print("</div>")
  2491. def printWarning(s):
  2492. print('<div style="text-align:left; background-color: #FFDDDD; padding:5px; border: 1px solid black"><strong>Warning:</strong>')
  2493. print(s)
  2494. print('</div>')
  2495. def printNoEffScoreFoundWarn(effScoresCount, pam):
  2496. if effScoresCount==0 and not pamIsCpf1(pam):
  2497. note = "No guide could be scored for efficiency. This happens when the input sequence is shorter than 100bp and there is no genome available to extend it or if there is simply not guide socring method. In the first case, please add flanking 50bp on both sides of the input sequence and submit this new, longer sequence. For the second case, you can contact me and suggest an efficiency scoring method, send me the published paper in this case."
  2498. printNote(note)
  2499. def showGuideTable(guideData, pam, otMatches, dbInfo, batchId, org, chrom, varHtmls):
  2500. " shows table of all PAM motif matches "
  2501. print("<br><div class='title'>Predicted guide sequences for PAMs</div>")
  2502. global scoreNames
  2503. if (cgiParams.get("showAllScores", "0")=="1"):
  2504. scoreNames = allScoreNames
  2505. showColumns = set()
  2506. # show the CFD guide score?
  2507. if pamIsSpCas9(pam):
  2508. showColumns.add("cfdGuideScore")
  2509. showPamWarning(pam)
  2510. showNoGenomeWarning(dbInfo)
  2511. printTableHead(pam, batchId, chrom, org, varHtmls, showColumns)
  2512. count = 0
  2513. effScoresCount = 0
  2514. showProxGcCol = ("proxGc" in scoreNames)
  2515. for guideRow in guideData:
  2516. guideScore, guideCfdScore, effScores, pamStart, guideStart, strand, pamId, guideSeq, \
  2517. pamSeq, otData, otDesc, last12Desc, mutEnzymes, ontargetDesc, repCount = guideRow
  2518. color = scoreToColor(guideScore)[0]
  2519. classStr = cssClassesFromSeq(guideSeq)
  2520. print('<tr id="%s" class="%s" style="border-left: 5px solid %s">' % (pamId, classStr, color))
  2521. # position and strand
  2522. #print '<td id="%s">' % pamId
  2523. print('<td>')
  2524. print('<a href="#list%s">' % (pamId))
  2525. print(str(pamStart+1)+" /")
  2526. if strand=="+":
  2527. print('fw')
  2528. else:
  2529. print('rev')
  2530. print('</a>')
  2531. print("</td>")
  2532. # sequence with variants and PCR primer link
  2533. print("<td>")
  2534. print("<small>")
  2535. # guide sequence + PAM sequence
  2536. if pamIsFirst:
  2537. fullGuideHtml = "<tt><i>"+pamSeq+"</i> " + guideSeq+"</tt>"
  2538. spacePos = len(pamSeq)
  2539. else:
  2540. fullGuideHtml = "<tt>"+guideSeq + " <i>" + pamSeq+"</i></tt>"
  2541. spacePos = len(guideSeq)
  2542. print(fullGuideHtml)
  2543. print("<br>")
  2544. # variant-string
  2545. if varHtmls is not None:
  2546. varFound = False
  2547. varStrs = []
  2548. guideHtmlStart = min(guideStart, pamStart)
  2549. guideHtmls = varHtmls[guideHtmlStart:guideHtmlStart+len(guideSeq)+len(pamSeq)]
  2550. if strand=="-":
  2551. guideHtmls = list(reversed(guideHtmls))
  2552. for i in range(len(guideHtmls)):
  2553. html = guideHtmls[i]
  2554. if html!=".":
  2555. varFound = True
  2556. if i==spacePos:
  2557. varStrs.append("&nbsp;")
  2558. varStrs.append(html)
  2559. print(("<tt style='color:#888888'>%s</tt><br>" % ("".join(varStrs))))
  2560. if "TTTT" in guideSeq.upper():
  2561. text = "This guide contains the sequence TTTT. It cannot be transcribed with a U6 or U3 promoter, as TTTT terminates the transcription."
  2562. htmlWarn(text)
  2563. print(' Not with U6/U3')
  2564. print("<br>")
  2565. if pam=="NGG":
  2566. grafType = crisporEffScores.getGrafType(guideSeq)
  2567. if grafType:
  2568. if grafType=="tt":
  2569. grafText = "The guide ends with TTC or TTT or contains only T and C in the last four nucleotides and more than 2 Ts or at least one TT and one T or C ('TT-motif'). These guides should be avoided in polymerase III (Pol III)-based gene editing experiments requiring high sgRNA expression levels."
  2570. elif grafType=="gcc":
  2571. grafText = "The guide ends with [AGT]GCC or GCCT ('GCC motif'). These sgRNAs appear to be inefficient irrespective of the delivery method and should thus be generally avoided."
  2572. text = "This guide contains one of the motifs described by <a target=_blank href='https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6352712/'>Graf et al, Cell Reports 2019</a>. %s " % grafText
  2573. htmlWarn(text)
  2574. print(' Inefficient')
  2575. print("<br>")
  2576. if gcContent(guideSeq)>0.75:
  2577. text = "This sequence has a GC content higher than 75%.<br>In the data of Tsai et al Nat Biotech 2015, the two guide sequences with a high GC content had almost as many off-targets as all other sequences combined. We do not recommend using guide sequences with such a high GC content."
  2578. htmlWarn(text)
  2579. print(' High GC content')
  2580. print("<br>")
  2581. if gcContent(guideSeq)<0.25:
  2582. text = "This sequence has a GC content lower than 25%.<br>In the data of Wang/Sabatini/Lander Science 2014, guides with a very low GC content had low cleavage efficiency."
  2583. htmlWarn(text)
  2584. print(' Low GC content<br>')
  2585. print("<br>")
  2586. if len(mutEnzymes)!=0:
  2587. print("<div style='margin-top: 3px'>Enzymes: <i>", end=' ')
  2588. print(", ".join([x.split("/")[0] for x,y,z in list(mutEnzymes.keys())]))
  2589. print("</i></div>")
  2590. scriptName = basename(__file__)
  2591. if otData!=None and repCount == 0:
  2592. print(('&nbsp;<a href="%s?batchId=%s&pamId=%s&pam=%s" target="_blank"><strong>Cloning / PCR primers</strong></a>' % (scriptName, batchId, urllib.parse.quote(str(pamId)), pam) ))
  2593. print("</small>")
  2594. print("</td>")
  2595. # off-target score, aka specificity score aka MIT score
  2596. if not pamIsCpf1(pam):
  2597. print("<td>")
  2598. if guideScore==None:
  2599. print("No matches")
  2600. else:
  2601. print("%d" % guideScore)
  2602. print("</td>")
  2603. # guide score based on CFD scores, aka guidescan score
  2604. if "cfdGuideScore" in showColumns:
  2605. print("<td>")
  2606. if guideCfdScore==None:
  2607. print("No matches")
  2608. else:
  2609. print("%d" % guideCfdScore)
  2610. print("</td>")
  2611. # eff scores
  2612. if effScores==None:
  2613. print('<td colspan="%d">Too close to end</td>' % len(scoreNames))
  2614. htmlHelp("The efficiency scores require some flanking sequence<br>This guide does not have enough flanking sequence in your input sequence and could not be extended as it was not found in the genome.<br>")
  2615. else:
  2616. for scoreName in scoreNames:
  2617. # out-of-frame and prox. gc need special treatment
  2618. if scoreName in ["oof", "proxGc"]:
  2619. continue
  2620. score = effScores.get(scoreName, None)
  2621. if score!=None:
  2622. effScoresCount += 1
  2623. if score==None:
  2624. print('''<td>--</td>''')
  2625. elif scoreName=="ssc":
  2626. # save some space
  2627. numStr = '%.1f' % (float(score))
  2628. print('''<td style="font-size:small">%s</td>''' % numStr)
  2629. elif scoreDigits.get(scoreName, 0)==0:
  2630. print('''<td>%d</td>''' % int(score))
  2631. else:
  2632. print('''<td>%0.1f</td>''' % (float(score)))
  2633. #print "<!-- %s -->" % seq30Mer
  2634. if showProxGcCol:
  2635. print("<td>")
  2636. # close GC > 4
  2637. finalGc = int(effScores.get("finalGc6", -1))
  2638. if finalGc==1:
  2639. print("+")
  2640. elif finalGc==0:
  2641. print("-")
  2642. else:
  2643. print("--")
  2644. # main motif is "NGG" and last nucleotides are GGNGG
  2645. if int(effScores.get("finalGg", 0))==1:
  2646. print("<br>")
  2647. print("<small>-GG</small>")
  2648. print("</td>")
  2649. if not baseEditor:
  2650. for mutScoreName in mutScoreNames:
  2651. print("<td>")
  2652. oofScore = str(effScores.get(mutScoreName, None))
  2653. if mutScoreName=="oof":
  2654. scoreDesc = "out-of-frame deletions"
  2655. else:
  2656. scoreDesc = "frameshift mutations"
  2657. if oofScore==None or oofScore=="None":
  2658. print("--")
  2659. else:
  2660. print("""<a href="%s?batchId=%s&pamId=%s&showMh=%s" target=_blank class="tooltipster" title="This score indicates how likely %s are. Click to show the induced deletions based on the micro-homology around the cleavage site.">%s</a>""" % (myName, batchId, urllib.parse.quote(pamId), mutScoreName, scoreDesc, oofScore))
  2661. #print """<br><br><small><a href="%s?batchId=%s&pamId=%s&showMh=1" target=_blank class="tooltipster">Micro-homology</a></small>""" % (myName, batchId, pamId)
  2662. print("</td>")
  2663. # mismatch description
  2664. print("<td>")
  2665. #otCount = sum([int(x) for x in otDesc.split("/")])
  2666. if otData==None:
  2667. # no genome match
  2668. print(otDesc)
  2669. htmlHelp("This exact sequence was not found in the genome.<br>If you have pasted a cDNA multi-exon sequence, note that sequences that overlap a splice site cannot be used as guide sequences. If you only have a cDNA sequence, please BLAST or BLAT your sequence first against the genome, then use the resulting exon from the genome for CRISPOR.<br>This warning also appears if you have selected the wrong or no genome.")
  2670. elif repCount > 0:
  2671. print ("Repeat")
  2672. htmlHelp("At <= 4 mismatches, %d alignments were found in the genome for this sequence, without looking at the PAM sequence around these alignments.<br>This guide is a repeated region, it is too unspecific.<br>Usually, CRISPR cannot be used to target repeats. Also, note that sequences that include long repeats will make the CRISPOR website slow. You can mask repeats with Ns to speed up the search." % repCount)
  2673. else:
  2674. print(otDesc)
  2675. print("<br>")
  2676. # mismatch description, last 12 bp
  2677. print('<small style="color:grey">'+last12Desc+"</small><br>")
  2678. otCount = len(otData)
  2679. print("<br><small>%d off-targets</small>" % otCount)
  2680. print("</td>")
  2681. # links to offtargets
  2682. print("<td><small>")
  2683. if otData!=None:
  2684. if len(otData)>500 and len(guideData)>1:
  2685. otData, cutoff = findOtCutoff(otData)
  2686. if cutoff==None:
  2687. print("More than 1000 off-targets, showing only top "+str(len(otData)))
  2688. else:
  2689. print("More than 500 off-targets, showing %d with score &gt;%0.1f " % (len(otData), cutoff))
  2690. htmlHelp("This guide sequence has a high number of off-targets, its use is discouraged.<br>To show all off-targets, paste only the guide sequence into the input sequence box.")
  2691. otLinks = makeOtBrowserLinks(otData, chrom, dbInfo, pamId)
  2692. print("\n".join(otLinks[:3]))
  2693. if len(otLinks)>3:
  2694. cssPamId = pamId.replace("-","minus").replace("+","plus") # +/-: not valid in css
  2695. cssPamId = cssPamId+"More"
  2696. print('<div id="%s" class="otMore" style="display:none; width:100%%">' % cssPamId)
  2697. print("\n".join(otLinks[3:]))
  2698. print('''<a style="float:right;text-decoration:underline" href="%s?batchId=%s&pamId=%s&otPrimers=1" id="%s">''' % (myName, batchId, urllib.parse.quote(pamId), cssPamId))
  2699. print('<strong>Off-target primers</strong></a>')
  2700. print('</div>')
  2701. print('''<a id="%sMoreLink" class="otMoreLink" onclick="showAllOts('%s')">''' % (cssPamId, cssPamId))
  2702. print('show all...</a>')
  2703. print('''<a id="%sLessLink" class="otLessLink" style="display:none" onclick="showLessOts('%s')">''' % (cssPamId, cssPamId))
  2704. print('show less...</a>')
  2705. print("</small></td>")
  2706. print("</tr>")
  2707. count = count+1
  2708. print("</table>")
  2709. printDownloadTableLinks(batchId, addTsv=True)
  2710. printNoEffScoreFoundWarn(effScoresCount, pam)
  2711. def linkLocalFiles(listFname):
  2712. """ write a <link> statement for each filename in listFname. Version them via mtime
  2713. (-> browser cache)
  2714. """
  2715. for fname in open(listFname).read().splitlines():
  2716. fname = fname.strip()
  2717. if not isfile(fname):
  2718. fname = join(HTMLDIR, fname)
  2719. if not isfile(fname):
  2720. print("missing: %s<br>" % fname)
  2721. continue
  2722. mTime = str(os.path.getmtime(fname)).split(".")[0] # seconds is enough
  2723. if fname.endswith(".css"):
  2724. #url = fname.replace("/var/www/", "http://tefor.net/")
  2725. print("<link rel='stylesheet' media='screen' type='text/css' href='%s%s?%s'/>" % (HTMLPREFIX, fname, mTime))
  2726. def printHeader(batchId, title):
  2727. " print the html header "
  2728. print('''<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">''')
  2729. print("<html><head>")
  2730. if title==None:
  2731. if batchName!="":
  2732. print("""<title>CRISPOR - %s</title>""" % batchName)
  2733. else:
  2734. print("""<title>CRISPOR</title>""")
  2735. else:
  2736. print("""<title>%s</title>""" % title)
  2737. print("""
  2738. <meta name='description' content='Design CRISPR guides with off-target and efficiency predictions, for more than 100 genomes.'/>
  2739. <meta http-equiv='Content-Type' content='text/html; charset=utf-8' />
  2740. <meta property='fb:admins' content='692090743' />
  2741. <meta name="google-site-verification" content="OV5GRHyp-xVaCc76rbCuFj-CIizy2Es0K3nN9FbIBig" />
  2742. <meta property='og:type' content='website' />
  2743. <meta property='og:url' content='http://crispor.gi.ucsc.edu/' />
  2744. <meta property='og:image' content='http://crispor.gi.ucsc.edu/image/CRISPOR.png' />
  2745. <script src="https://cdn.jsdelivr.net/npm/clipboard@2/dist/clipboard.min.js"></script>
  2746. """)
  2747. # load jquery from local copy, not from CDN, for offline use
  2748. print(("""<script src='%sjs/jquery.min.js'></script>
  2749. <script src='%sjs/jquery-ui.min.js'></script>
  2750. """ % (HTMLPREFIX, HTMLPREFIX)))
  2751. #print('<link rel="stylesheet" href="//fonts.googleapis.com/css?family=Roboto:300,300italic,700,700italic" />')
  2752. #print('<link rel="stylesheet" type="text/css" href="https://cdnjs.cloudflare.com/ajax/libs/normalize/5.0.0/normalize.min.css" />')
  2753. #print('<link rel="stylesheet" href="//cdn.rawgit.com/milligram/milligram/master/dist/milligram.min.css">')
  2754. linkLocalFiles("includes.txt")
  2755. print('<link rel="stylesheet" type="text/css" href="%sstyle/tooltipster.css" />' % HTMLPREFIX)
  2756. print('<link rel="stylesheet" type="text/css" href="%sstyle/tooltipster-shadow.css" />' % HTMLPREFIX)
  2757. print('<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/chosen/1.6.2/chosen.css" />')
  2758. print('<link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/source-code-pro@2.38.0/source-code-pro.css" />')
  2759. # the UFD combobox, https://code.google.com/p/ufd/wiki/Usage
  2760. # patched to allow mouse wheel
  2761. # https://code.google.com/p/ufd/issues/detail?id=86&q=mouse%20wheel
  2762. print('<script type="text/javascript" src="%sjs/jquery.ui.ufd.js"></script>' % HTMLPREFIX)
  2763. #print '<link rel="stylesheet" type="text/css" href="%sstyle/ufd-base.css" />' % HTMLPREFIX
  2764. print('<link rel="stylesheet" type="text/css" href="%sstyle/plain.css" />' % HTMLPREFIX)
  2765. print('<link rel="stylesheet" type="text/css" href="%sstyle/jquery-ui.css" />' % HTMLPREFIX)
  2766. print('<script type="text/javascript" src="js/jquery.tooltipster.min.js"></script>')
  2767. print('<script src="https://cdnjs.cloudflare.com/ajax/libs/chosen/1.6.2/chosen.jquery.min.js"></script>')
  2768. # override the main TEFOR css
  2769. print("""
  2770. <style>
  2771. select { font-size: 80%; }
  2772. body {
  2773. text-align: left;
  2774. /* float: left; */
  2775. }
  2776. p {
  2777. }
  2778. ul {
  2779. -webkit-margin-before: 0;
  2780. -webkit-margin-after: 0;
  2781. }
  2782. .editTable {
  2783. border: 1px solid black;
  2784. background-color: white;
  2785. }
  2786. mut {
  2787. color: blue;
  2788. background-color: yellow;
  2789. }
  2790. .editTable th {
  2791. background-color: #F0F0F0;
  2792. }
  2793. tt { font-size: 90% }
  2794. div.contentcentral { text-align: left; float: left}
  2795. /* for chosen.js */
  2796. .chosen-container { width: 600px }
  2797. .chosen-container .chosen-results li.active-result { float: left}
  2798. </style>
  2799. """)
  2800. # activate tooltipster
  2801. #theme: 'tooltipster-shadow',
  2802. # activate jqueryUI tooltips
  2803. print ("""
  2804. <script>
  2805. $(function () {
  2806. $(".tooltip").tooltip({
  2807. relative : true,
  2808. tooltipClass : "alignStyle",
  2809. content: function () {
  2810. return '<div style="width:300px">'+$(this).prop('title')+"</div>";
  2811. }
  2812. });
  2813. });
  2814. $(function () {
  2815. $(".tooltipAuto").tooltip({
  2816. contentAsHtml : true
  2817. });
  2818. });
  2819. </script>""")
  2820. # style of Jquery UI tooltips, default style is div.ui-tooltip
  2821. print("""<style>
  2822. .alignStyle {
  2823. background-color: #FFFFFF;
  2824. width: 350px;
  2825. max-width: 400px;
  2826. height: 110px;
  2827. position : absolute;
  2828. text-align: left;
  2829. border:1px solid #cccccc;
  2830. }
  2831. </style>""")
  2832. # style from https://css-tricks.com/rotated-table-column-headers/ to rotate table headers
  2833. print("""<style>
  2834. th.rotate {
  2835. /* Something you can count on */
  2836. /* height: 10px; */
  2837. white-space: nowrap;
  2838. }
  2839. th.rotate > div {
  2840. float:left;
  2841. white-space: nowrap;
  2842. position: relative;
  2843. border-style: none;
  2844. """)
  2845. # if we're showing all scores, we have very little space, so turn by
  2846. # 90degrees. otherwise, we can afford 45 degrees, which is easier to read.
  2847. #if cgiParams.get("showAllScores", "0")=="1":
  2848. print("""
  2849. -webkit-transform: rotate(-90);
  2850. -moz-transform: rotate(270deg);
  2851. -ms-transform: rotate(270deg);
  2852. -o-transform: rotate(270deg);
  2853. transform: rotate(270deg);
  2854. width: 25px;""")
  2855. #else:
  2856. #print("""
  2857. #-webkit-transform: rotate(-45deg);
  2858. #-moz-transform: rotate(315deg);
  2859. #-ms-transform: rotate(315deg);
  2860. #-o-transform: rotate(315deg);
  2861. #transform: rotate(315deg);
  2862. #width: 25px;""")
  2863. print("""
  2864. }
  2865. th.rotate > div > span {
  2866. /* border-bottom: 1px solid #ccc; */
  2867. padding: 0px 3px;
  2868. white-space: nowrap;
  2869. }
  2870. </style>""")
  2871. print("</head>")
  2872. print('<body id="wrapper">')
  2873. def firstFreeLine(lineMasks, y, start, end):
  2874. " recursively search for first free line to place a feature (start, end) "
  2875. #print "first free line called with y", y, "<br>"
  2876. if y>=len(lineMasks):
  2877. return None
  2878. lineMask = lineMasks[y]
  2879. for x in range(start, end):
  2880. #print "checking pos", x, "<br>"
  2881. if lineMask[x]!=0:
  2882. return firstFreeLine(lineMasks, y+1, start, end)
  2883. return y
  2884. #return None
  2885. def distrOnLines(seq, startDict, featLen, pam):
  2886. """ given a dict with start -> (start,end,name,strand) and a motif len, create lines of annotations such that
  2887. the motifs don't overlap on the lines
  2888. """
  2889. # max number of lines in y direction to draw
  2890. MAXLINES = 18
  2891. # amount of free space around each feature
  2892. SLOP = 2
  2893. # bitmask, one per line, 1 = we have a feature here, 0 = no feature here
  2894. lineMasks = []
  2895. for i in range(0, MAXLINES):
  2896. lineMasks.append( [0]* (len(seq)+10) )
  2897. # dict with lineCount (0...MAXLINES) -> list of (start, strand) tuples
  2898. ftsByLine = defaultdict(list)
  2899. maxY = 0
  2900. for start in sorted(startDict):
  2901. end = start+featLen
  2902. strand = startDict[start]
  2903. # Cannot use Unicode here: these symbols are not part of the
  2904. # monospace font on some platforms and therefore their width
  2905. # is not the same as the other characters
  2906. #arrNE = u'\u2197'
  2907. #arrSE = u'\u2198'
  2908. arrNE = '/'
  2909. arrSE = '\\'
  2910. #arrNE = u'\u2a3c' # hebrew
  2911. #arrSE = u'\ufb27' # math
  2912. ftSeq = seq[start:end]
  2913. if strand=="+":
  2914. if pamIsFirst:
  2915. if pamIsCas12max(pam): #modify cleavage sites to 14-16 (target strand) and 24nt (non-target strand) / hfCas12Max
  2916. label = '%s'%(ftSeq)+'.............%s%s%s.......%s' % (arrNE, arrNE, arrNE, arrSE)
  2917. else:
  2918. label = '%s'%(ftSeq)+'.................%s....%s' % (arrNE, arrSE)
  2919. startFt = start
  2920. endFt = start+len(label)
  2921. else:
  2922. #label = '%s..%s'%(seq[start-3].lower(), ftSeq)
  2923. #label = '---%s'%(ftSeq)
  2924. #label = '&#45;&#45;&#45;%s'%(ftSeq)
  2925. label = '&#8722;&#8722;&#8722;%s'%(ftSeq)
  2926. startFt = start - 3
  2927. endFt = end
  2928. else:
  2929. if pamIsFirst:
  2930. if pamIsCas12max(pam): #modify cleavage sites to 14-16 (target strand) and 24nt (non-target strand) / hfCas12Max
  2931. spc1 = "......."
  2932. spc2 = "............."
  2933. labelPrefix = '%s%s%s%s%s%s' % (arrSE, spc1, arrNE, arrNE, arrNE, spc2)
  2934. label = labelPrefix + ftSeq
  2935. else:
  2936. spc1 = "...."
  2937. spc2 = "................."
  2938. labelPrefix = '%s%s%s%s' % (arrSE, spc1, arrNE, spc2)
  2939. label = labelPrefix + ftSeq
  2940. startFt = start - len(labelPrefix)
  2941. endFt = startFt+len(label)
  2942. else:
  2943. #label = '%s..%s'%(ftSeq, seq[end+2].lower())
  2944. label = '%s&#45;&#45;&#45;'%(ftSeq)
  2945. startFt = start
  2946. endFt = end + 3
  2947. #print "feature", strand, start, startFt, endFt, SLOP,"<br>"
  2948. #print "mask", lineMasks[0][startFt:endFt], "<br>"
  2949. y = firstFreeLine(lineMasks, 0, startFt, endFt)
  2950. #print "free line: %s<br>" % y
  2951. if y==None:
  2952. errAbort("not enough space to plot features")
  2953. # fill the current mask
  2954. mask = lineMasks[y]
  2955. maskStart = max(startFt-SLOP, 0)
  2956. maskEnd = min(endFt+SLOP, len(seq))
  2957. #print "mask:", maskStart, maskEnd
  2958. for i in range(maskStart, maskEnd):
  2959. mask[i]=1
  2960. maxY = max(y, maxY)
  2961. pamId = "s%d%s" % (start, strand)
  2962. ft = (startFt, endFt, label, strand, pamId)
  2963. #print "labelLen: %d<br>" % len(label)
  2964. #print "ft: %s<br>" % repr(ft)
  2965. ftsByLine[y].append(ft )
  2966. return ftsByLine, maxY
  2967. def writePamFlank(seq, startDict, pam, faFname):
  2968. " write pam flanking sequences to fasta file, optionally with versions where each nucl is removed "
  2969. #print "writing pams to %s<br>" % faFname
  2970. faFh = open(faFname, "w")
  2971. for pamId, pamStart, guideStart, strand, flankSeq, pamSeq, pamPlusSeq in flankSeqIter(seq, startDict, len(pam), True):
  2972. faFh.write(">%s\n%s\n" % (pamId, flankSeq))
  2973. faFh.close()
  2974. def runCmd(cmd, ignoreExitCode=False, useShell=True):
  2975. " run shell command, check ret code, replaces BIN and SCRIPTS special variables "
  2976. if useShell:
  2977. cmd = cmd.replace("$BIN", binDir)
  2978. cmd = cmd.replace("$PYTHON", sys.executable)
  2979. cmd = cmd.replace("$SCRIPT", scriptDir)
  2980. cmd = "set -o pipefail; " + cmd
  2981. executable = "/bin/bash"
  2982. else:
  2983. cmd = [x.replace("$BIN", binDir).replace("$PYTHON", sys.executable).replace("$SCRIPT", scriptDir) for x in cmd]
  2984. executable=None
  2985. debug("Running %s" % cmd)
  2986. ret = subprocess.call(cmd, shell=useShell, executable=executable)
  2987. if ret!=0 and not ignoreExitCode:
  2988. if not useShell:
  2989. cmd = " ".join(cmd)
  2990. if commandLineMode:
  2991. logging.error("Error: could not run command %s." % cmd)
  2992. sys.exit(1)
  2993. else:
  2994. print("Server error: could not run command %s, error %d.<p>" % (cmd, ret))
  2995. print("please send us an email, we will fix this error as quickly as possible. %s " % contactEmail)
  2996. raise
  2997. sys.exit(0)
  2998. def isAltChrom(chrom):
  2999. """ return true is chrom name looks like it's not on the primary assembly. This is mostly relevant for hg38.
  3000. examples: chr6_*_alt (hg38)
  3001. """
  3002. return chrom.endswith("_alt")
  3003. def parseOfftargets(db, batchId, onTargetChrom=""):
  3004. """ parse a bed file with annotataed off target matches from overlapSelect,
  3005. has two name fields, one with the pam position/strand and one with the
  3006. overlapped segment
  3007. return as dict pamId -> editDist -> (chrom, start, end, seq, strand, segType, segName, totalAlnCount, isRep)
  3008. segType is "ex" "int" or "ig" (=intergenic)
  3009. if intergenic, geneNameStr is two genes, split by |
  3010. The isRep flag is true if BWA reported more than one alignment with X0+X1 but didn't report these with the
  3011. XA tag. It means that we can't get the alignments for this sequence from BWA (=repeats).
  3012. """
  3013. # edge case: target is on chr6_alt -> we remove a single off-target with 0 mismatches on chr6
  3014. targetIsAlt = isAltChrom(onTargetChrom)
  3015. # keep track of pamIds already handled for this edge case
  3016. skippedPams = set()
  3017. # ideally we would check if the offtarget falls into the chrom area that gave rise to the alt
  3018. # but that would mean parsing yet another non-small file and this case should be sufficiently rare
  3019. # to not bother 99% of users.
  3020. batchBase = join(batchDir, batchId)
  3021. bedFname = batchBase+".bed.gz"
  3022. # example input:
  3023. # chrIV 9864393 9864410 s41-|-|5|ACTTGACTG|0 chrIV 9864303 9864408 ex:K07F5.16
  3024. # chrIV 9864393 9864410 s41-|-|5|ACTGTAGCTAGCT|9999 chrIV 9864408 9864470 in:K07F5.16
  3025. debug("reading offtargets from %s" % bedFname)
  3026. # first sort into dict (pamId,chrom,start,end,editDist,strand)
  3027. # -> (segType, segName)
  3028. pamData = {}
  3029. #ifh = open(bedFname) # switched to gzip compression in Dec 2018, converted old files with bash script
  3030. try:
  3031. ifh = gzip.open(bedFname, "rt")
  3032. except FileNotFoundError:
  3033. print("Off-target results from this link were temporarily removed to save space. ")
  3034. linkUrl = "crispor.py?batchId="+batchId
  3035. print("Please go back to <a href='%s'>your job page</a>, which will rerun the job, then try the link again." % linkUrl)
  3036. exit(0)
  3037. maxOtLines = 500000
  3038. count = 0
  3039. for line in ifh:
  3040. fields = line.rstrip("\n").split("\t")
  3041. count +=1
  3042. if count > maxOtLines:
  3043. print("Error: More than %d off-targets. CRISPR has trouble with handling extremely unspecific inputs. Please email us at "
  3044. "%s and discuss. Your input sequence most likely includes a recent L1HS, SVA or similar repetitive elements. "
  3045. "It is hard to design guides for these, you can try to re-run your input sequence, but without the repetitive element. "%
  3046. (maxOtLines, contactEmail))
  3047. sys.exit(1)
  3048. chrom, start, end, name, segment = fields
  3049. logging.debug("off-target: %s" % name)
  3050. # hg38: ignore alternate chromosomes otherwise the
  3051. # regions on the main chroms look as if they could not be
  3052. # targeted at all with Cas9
  3053. if isAltChrom(chrom):
  3054. logging.debug("skipping off-target: on alt-chromosome")
  3055. continue
  3056. nameFields = name.split("|")
  3057. pamId, strand, editDist, seq = nameFields[:4]
  3058. #print pamId, strand, editDist, seq, chrom, start, end, name, segment, "<br>"
  3059. if targetIsAlt:
  3060. if editDist=='0' and onTargetChrom.split("_")[0]==chrom.split("_")[0] and not pamId in skippedPams:
  3061. logging.debug("altChrom edge case: target is on alt-chrom, skipping a single 0-mismatch off-target on primary chrom")
  3062. skippedPams.add(pamId)
  3063. continue
  3064. isRep = 0
  3065. totalAlnCount = 0
  3066. # for compatibility with old bed files, only parse these fields if they are present
  3067. # note: for some reason, the MIT hitScore was always written to these files
  3068. # However, it's not parsed here and never was. In order to not break the old files
  3069. # I kept it in the files, but am not reading it here
  3070. # these are the different formats until now:
  3071. # seqId+"|"+strand+"|"+editDist+"|"+seq+"|"+str(hitScore) # five fields
  3072. # seqId+"|"+strand+"|"+editDist+"|"+seq+"|"+x1Score+"|"+str(hitScore) # six fields
  3073. # (x1Score was roughly the alnCount, similar enough for practical purposes, fixed in 2019)
  3074. # seqId+"|"+strand+"|"+editDist+"|"+seq+"|"+alnCount+"|"+str(hitScore)+"|"+isRep # seven fields
  3075. if len(nameFields)>5:
  3076. totalAlnCount = int(nameFields[4])
  3077. if len(nameFields)>6:
  3078. isRep = bool(int(nameFields[6]))
  3079. editDist = int(editDist)
  3080. # some gene models include colons
  3081. if ":" in segment:
  3082. segType, segName = segment.split(":", maxsplit=1)
  3083. else:
  3084. segType, segName = "", segment
  3085. start, end = int(start), int(end)
  3086. otKey = (pamId, chrom, start, end, editDist, seq, strand, totalAlnCount, isRep)
  3087. # if an offtarget is in the PAR region, we keep only the chrY off-target
  3088. parNum = isInPar(db, chrom, start, end)
  3089. # keep only matches on chrX
  3090. if parNum is not None and chrom=="chrX":
  3091. logging.debug("off-target on PAR region, skipping")
  3092. continue
  3093. # if a offtarget overlaps an intron/exon or ig/exon boundary it will
  3094. # appear twice; in this case, we only keep the exon offtarget
  3095. if otKey in pamData and segType!="ex":
  3096. logging.debug("skipping off-target: ex/ig boundary edge case")
  3097. continue
  3098. pamData[otKey] = (segType, segName)
  3099. # index by pamId and edit distance
  3100. indexedOts = defaultdict(dict)
  3101. for otKey, otVal in pamData.items():
  3102. pamId, chrom, start, end, editDist, seq, strand, totalAlnCount, isRep = otKey
  3103. segType, segName = otVal
  3104. otTuple = (chrom, start, end, seq, strand, segType, segName, totalAlnCount, isRep)
  3105. indexedOts[pamId].setdefault(editDist, []).append( otTuple )
  3106. return indexedOts
  3107. class ConsQueue:
  3108. """ a pseudo job queue that does nothing but report progress to the console """
  3109. def startStep(self, batchId, desc, label):
  3110. logging.info("Progress %s - %s - %s" % (batchId, desc, label))
  3111. def annotateBedWithPos(inBed, outBed, genome):
  3112. """
  3113. given an input bed4 and an output bed filename, add an additional column 5 to the bed file
  3114. that is a descriptive text of the chromosome pos (e.g. chr1:1.23 Mbp).
  3115. """
  3116. ofh = gzip.open(outBed, "wt")
  3117. for line in open(inBed):
  3118. chrom, start = line.split("\t")[:2]
  3119. chrom = applyChromAlias(genome, chrom)
  3120. start = int(start)
  3121. if start>1000000:
  3122. startStr = "%.2f Mbp" % (float(start)/1000000)
  3123. else:
  3124. startStr = "%.2f Kbp" % (float(start)/1000)
  3125. desc = "%s %s" % (chrom, startStr)
  3126. ofh.write(line.rstrip("\n"))
  3127. ofh.write("\t")
  3128. ofh.write(desc)
  3129. ofh.write("\n")
  3130. ofh.close()
  3131. def findAllGuides(seq, pam):
  3132. startDict, endSet = findAllPams(seq, pam)
  3133. pamInfo = list(flankSeqIter(seq, startDict, len(pam), False))
  3134. return pamInfo
  3135. def extractMutScores(scoreDict, pamIds):
  3136. " make a list of the guide-related outcome scores in the order of pamIds "
  3137. res = []
  3138. for pamId in pamIds:
  3139. res.append(scoreDict[pamId][0])
  3140. return res
  3141. def calcSaveEffScores(batchId, seq, extSeq, pam, queue):
  3142. """ given a sequence and an extended sequence, get all potential guides
  3143. with pam, extend them to 100mers and score them with various eff. scores.
  3144. Return a
  3145. list of rows [headers, (guideSeq, 100mer, score1, score2, score3,...), ... ]
  3146. Also write the results to a database so they can be retrieved later.
  3147. extSeq can be None, if we were unable to extend the sequence
  3148. """
  3149. seq = seq.upper()
  3150. if extSeq:
  3151. extSeq = extSeq.upper()
  3152. pamInfo = findAllGuides(seq, pam)
  3153. pamIds = []
  3154. guides = []
  3155. longSeqs = []
  3156. for pamId, startPos, guideStart, strand, guideSeq, pamSeq, pamPlusSeq in pamInfo:
  3157. logging.debug("PAM ID: %s - guideSeq %s" % (pamId, guideSeq))
  3158. gStart, gEnd = pamStartToGuideRange(startPos, strand, len(pam))
  3159. longSeq = getExtSeq(seq, gStart, gEnd, strand, 50-GUIDELEN, 50, extSeq) # +-50 bp from the end of the guide
  3160. if longSeq!=None:
  3161. longSeqs.append(longSeq)
  3162. pamIds.append(pamId)
  3163. guides.append(guideSeq+pamSeq)
  3164. if len(longSeqs)>0 and doEffScoring:
  3165. enz = None
  3166. if pamIsCpf1(pam) and not pam=="NGTN":
  3167. enz = "cpf1"
  3168. elif pamIsSaCas9(pam):
  3169. enz = "sacas9"
  3170. # for spcas9, we use the extended list for the calculation
  3171. global scoreNames
  3172. if enz is None:
  3173. scoreNames = allScoreNames
  3174. effScores = crisporEffScores.calcAllScores(longSeqs, enzyme=enz, scoreNames=scoreNames)
  3175. # these are slow algorithms, so store the results for later
  3176. queue.startStep(batchId, "outcome", "Calculating editing outcomes")
  3177. mutScores = crisporEffScores.calcMutSeqs(pamIds, longSeqs, enz, scoreNames=mutScoreNames)
  3178. saveOutcomeData(batchId, mutScores)
  3179. # for output and sorting, it's easier to treat the outcome-derived scores like an efficiency score
  3180. for mutScoreName in mutScoreNames:
  3181. if mutScoreName in mutScores:
  3182. effScores[mutScoreName] = extractMutScores(mutScores[mutScoreName], pamIds)
  3183. # make sure the "N bug" reported by Alberto does never happen again:
  3184. # we must get back as many scores as we have sequences
  3185. for scoreName, scores in effScores.items():
  3186. if len(scores)!=len(longSeqs):
  3187. print("Internal error when calculating score %s" % scoreName)
  3188. assert(False)
  3189. else:
  3190. effScores = {}
  3191. activeScoreNames = list(effScores.keys())
  3192. # reformat to rows, write all scores to file
  3193. assert(len(pamIds)==len(guides)==len(longSeqs))
  3194. rows = []
  3195. for i, (guideId, guide, longSeq) in enumerate(zip(pamIds, guides, longSeqs)):
  3196. row = [guideId, guide, longSeq]
  3197. for scoreName in activeScoreNames:
  3198. scoreList = effScores[scoreName]
  3199. if len(scoreList) > 0:
  3200. row.append(scoreList[i])
  3201. else:
  3202. row.append("noScore?")
  3203. rows.append(row)
  3204. headerRow = ["guideId", "guide", "longSeq"]
  3205. headerRow.extend(activeScoreNames)
  3206. rows.insert(0, headerRow)
  3207. return rows
  3208. def writeRow(ofh, row):
  3209. " write list to file as tab-sep row "
  3210. row = [str(x) for x in row]
  3211. ofh.write("\t".join(row))
  3212. ofh.write("\n")
  3213. def createBatchEffScoreTable(batchId, queue):
  3214. """ annotate all potential guides with efficiency scores and write to file.
  3215. tab-sep file for easier debugging, no pickling
  3216. """
  3217. outFname = join(batchDir, batchId+".effScores.tab")
  3218. # Todo: why don't we get these from the caller as arguments instead of reading the batch?
  3219. batchInfo = readBatchAsDict(batchId)
  3220. seq = batchInfo["seq"]
  3221. extSeq = batchInfo.get("extSeq") # cannot always extend a sequence, e.g. when no perfect match
  3222. pam = batchInfo["pam"]
  3223. pam = setupPamInfo(pam)
  3224. seq = seq.upper()
  3225. if extSeq:
  3226. extSeq = extSeq.upper()
  3227. guideRows = calcSaveEffScores(batchId, seq, extSeq, pam, queue)
  3228. guideFh = open(outFname, "w")
  3229. for row in guideRows:
  3230. writeRow(guideFh, row)
  3231. guideFh.close()
  3232. logging.info("Wrote eff scores to %s" % guideFh.name)
  3233. def readEffScores(batchId):
  3234. " parse eff scores from tab sep file and return as dict pamId -> dict of scoreName -> value "
  3235. effScoreFname = join(batchDir, batchId)+".effScores.tab"
  3236. seqToScores = {}
  3237. if isfile(effScoreFname):
  3238. for row in lineFileNext(open(effScoreFname)):
  3239. scoreDict = {}
  3240. rowDict = row._asdict()
  3241. # the first three fields are the pamId, shortSeq, longSeq, they are not scores
  3242. allScoreNames = row._fields[3:]
  3243. for scoreName in allScoreNames:
  3244. score = rowDict[scoreName]
  3245. if score=="None":
  3246. score = "NA"
  3247. elif "." in score or "e" in score:
  3248. score = float(score)
  3249. else:
  3250. score = int(score)
  3251. scoreDict[scoreName] = score
  3252. seqToScores[row.guideId] = scoreDict
  3253. return seqToScores
  3254. def findOfftargetsBwa(queue, batchId, batchBase, faFname, genome, pamDesc, bedFname):
  3255. " align faFname to genome and create matchedBedFname "
  3256. matchesBedFname = batchBase+".matches.bed"
  3257. saFname = batchBase+".sa"
  3258. pam = setupPamInfo(pamDesc)
  3259. pamLen = len(pam)
  3260. genomeDir = genomesDir # make var local, see below
  3261. open(matchesBedFname, "w") # truncate to 0 size
  3262. # increase MAXOCC if there is only a single query, but only in CGI mode
  3263. #if len(parseFasta(open(faFname)))==1 and not commandLineMode:
  3264. #global MAXOCC
  3265. #global maxMMs
  3266. #MAXOCC=max(HIGH_MAXOCC, MAXOCC)
  3267. #maxMMs=HIGH_maxMMs
  3268. maxDiff = maxMMs
  3269. queue.startStep(batchId, "bwa", "Alignment of potential guides, mismatches <= %d" % maxDiff)
  3270. convertMsg = "Converting alignments"
  3271. seqLen = GUIDELEN
  3272. bwaM = MFAC*MAXOCC # -m is queue size in bwa
  3273. cmd = "$BIN/bwa aln -o 0 -m %(bwaM)s -n %(maxDiff)d -k %(maxDiff)d -N -l %(seqLen)d %(genomeDir)s/%(genome)s/%(genome)s.fa %(faFname)s > %(saFname)s" % locals()
  3274. runCmd(cmd)
  3275. queue.startStep(batchId, "saiToBed", convertMsg)
  3276. maxOcc = MAXOCC # create local var from global
  3277. # EXTRACTION OF POSITIONS + CONVERSION + SORT/CLIP
  3278. # the sorting should improve the twoBitToFa runtime
  3279. python = sys.executable
  3280. cmd = "$BIN/bwa samse -n %(maxOcc)d %(genomeDir)s/%(genome)s/%(genome)s.fa %(saFname)s %(faFname)s | $SCRIPT/xa2multi.pl | %(python)s $SCRIPT/samToBed %(pam)s %(seqLen)d | sort -k1,1 -k2,2n | $BIN/bedClip stdin %(genomeDir)s/%(genome)s/%(genome)s.sizes stdout >> %(matchesBedFname)s " % locals()
  3281. runCmd(cmd)
  3282. filtMatchesBedFname = batchBase+".filtMatches.bed"
  3283. queue.startStep(batchId, "filter", "Removing matches without a PAM motif")
  3284. altPats = ",".join(offtargetPams.get(pam, ["na"]))
  3285. bedFnameTmp = bedFname+".tmp"
  3286. altPamMinScore = str(ALTPAMMINSCORE)
  3287. shmFaFname = join("/dev/shm", genome+".fa")
  3288. # EXTRACTION OF SEQUENCES + ANNOTATION - big headache!!
  3289. # twoBitToFa was 15x slower than python's twobitreader, after markd's fix it is better
  3290. # but bedtools uses an fa.idx file and also mmap, so is a LOT faster
  3291. # arguments: guideSeq, mainPat, altPats, altScore, passTotalAlnCount
  3292. if isfile(shmFaFname):
  3293. logging.info("Using bedtools and genome fasta on ramdisk, %s" % shmFaFname)
  3294. cmd = "time bedtools getfasta -s -name -fi %(shmFaFname)s -bed %(matchesBedFname)s -fo /dev/stdout | $SCRIPT/filterFaToBed %(faFname)s %(pam)s %(altPats)s %(altPamMinScore)s > %(filtMatchesBedFname)s" % locals()
  3295. else:
  3296. cmd = "time $BIN/twoBitToFa %(genomeDir)s/%(genome)s/%(genome)s.2bit stdout -bed=%(matchesBedFname)s | %(python)s $SCRIPT/filterFaToBed %(faFname)s %(pam)s %(altPats)s %(altPamMinScore)s > %(filtMatchesBedFname)s" % locals()
  3297. #cmd = "$SCRIPT/twoBitToFaPython %(genomeDir)s/%(genome)s/%(genome)s.2bit %(matchesBedFname)s | $SCRIPT/filterFaToBed %(faFname)s %(pam)s %(altPats)s %(altPamMinScore)s %(maxOcc)d > %(filtMatchesBedFname)s" % locals()
  3298. runCmd(cmd)
  3299. segFname = "%(genomeDir)s/%(genome)s/%(genome)s.segments.bed" % locals()
  3300. # if we have gene model segments, annotate them, otherwise just use the chrom position
  3301. if isfile(segFname):
  3302. queue.startStep(batchId, "genes", "Annotating matches with genes")
  3303. cmd = "cat %(filtMatchesBedFname)s | $BIN/overlapSelect %(segFname)s stdin stdout -mergeOutput -selectFmt=bed -inFmt=bed | cut -f1,2,3,4,8 | gzip > %(bedFnameTmp)s " % locals()
  3304. runCmd(cmd)
  3305. else:
  3306. queue.startStep(batchId, "chromPos", "Annotating matches with chromosome position")
  3307. annotateBedWithPos(filtMatchesBedFname, bedFnameTmp, genome)
  3308. # make sure the final bed file is never in a half-written state,
  3309. # as it is our signal that the job is complete
  3310. shutil.move(bedFnameTmp, bedFname)
  3311. queue.startStep(batchId, "done", "Job completed")
  3312. # remove the temporary files
  3313. tempFnames = [saFname, matchesBedFname, filtMatchesBedFname]
  3314. if not DEBUG:
  3315. for tfn in tempFnames:
  3316. if isfile(tfn):
  3317. os.remove(tfn)
  3318. return bedFname
  3319. def makeVariants(seq):
  3320. " generate all possible variants of sequence at 1bp-distance"
  3321. seqs = []
  3322. for i in range(0, len(seq)):
  3323. for l in "ACTG":
  3324. if l==seq[i]:
  3325. continue
  3326. newSeq = seq[:i]+l+seq[i+1:]
  3327. seqs.append((i, seq[i], l, newSeq))
  3328. return seqs
  3329. def expandIupac(seq):
  3330. """ expand all IUPAC characters to nucleotides, returns list.
  3331. >>> expandIupac("NY")
  3332. ['GC', 'GT', 'AC', 'AT', 'TC', 'TT', 'CC', 'CT']
  3333. """
  3334. # http://stackoverflow.com/questions/27551921/how-to-extend-ambiguous-dna-sequence
  3335. d = {'A': 'A', 'C': 'C', 'B': 'CGT', 'D': 'AGT', 'G': 'G', \
  3336. 'H': 'ACT', 'K': 'GT', 'M': 'AC', 'N': 'GATC', 'S': 'CG', \
  3337. 'R': 'AG', 'T': 'T', 'W': 'AT', 'V': 'ACG', 'Y': 'CT', 'X': 'GATC'}
  3338. seqs = []
  3339. for i in product(*[d[j] for j in seq]):
  3340. seqs.append("".join(i))
  3341. return seqs
  3342. def writeBowtieSequences(inFaFname, outFname, pamPat):
  3343. """ write the sequence and one-bp-distant-sequences + all possible PAM sequences to outFname
  3344. Return dict querySeqId -> querySeq and a list of all
  3345. possible PAMs, as nucleotide sequences (not IUPAC-patterns)
  3346. """
  3347. ofh = open(outFname, "w")
  3348. outCount = 0
  3349. inCount = 0
  3350. guideSeqs = {} # 20mer guide sequences
  3351. qSeqs = {} # 23mer query sequences for bowtie, produced by expanding guide sequences
  3352. allPamSeqs = expandIupac(pamPat)
  3353. for seqId, seq in parseFastaAsList(open(inFaFname)):
  3354. inCount += 1
  3355. guideSeqs[seqId] = seq
  3356. for pamSeq in allPamSeqs:
  3357. # the input sequence + the PAM
  3358. newSeqId = "%s.%s" % (seqId, pamSeq)
  3359. newFullSeq = seq+pamSeq
  3360. ofh.write(">%s\n%s\n" % (newSeqId, newFullSeq))
  3361. qSeqs[newSeqId] = newFullSeq
  3362. # all one-bp mutations of the input sequence + the PAM
  3363. for nPos, fromNucl, toNucl, newSeq in makeVariants(seq):
  3364. newSeqId = "%s.%s.%d:%s>%s" % (seqId, pamSeq, nPos, fromNucl, toNucl)
  3365. newFullSeq = newSeq+pamSeq
  3366. ofh.write(">%s\n%s\n" % (newSeqId, newFullSeq))
  3367. qSeqs[newSeqId] = newFullSeq
  3368. outCount += 1
  3369. ofh.close()
  3370. logging.debug("Wrote %d variants+expandedPam of %d sequences to %s" % (outCount, inCount, outFname))
  3371. return guideSeqs, qSeqs, allPamSeqs
  3372. def applyModifStr(seq, modifStrs, strand):
  3373. """ bowtie: given a list of pos:toNucl>fromNucl and a seq, return the original seq.
  3374. position is 0-based
  3375. >>> applyModifStr("ACAATAAGACATAAACATATCGG", "14:T>A,21:A>G,22:C>G".split(","), "+")
  3376. 'ACAATAAGACATAATCATATCAC'
  3377. """
  3378. seq = list(seq)
  3379. for modifStr in modifStrs:
  3380. #logging.debug( modifStr)
  3381. pos, toFromNucl = modifStr.split(":")
  3382. fromNucl, toNucl = toFromNucl.split(">")
  3383. pos = int(pos)
  3384. if strand=="-":
  3385. fromNucl = revComp(fromNucl)
  3386. seq[pos] = fromNucl
  3387. return "".join(seq)
  3388. def parseRefout(tmpDir, guideSeqs, pamLen):
  3389. """ parse all .map file in tmpDir and return as list of chrom,start,end,strand,guideSeq,tSeq
  3390. """
  3391. fnames = glob.glob(join(tmpDir, "*.map"))
  3392. # while parsing, make sure we keep only the hit with the lowest number of mismatches
  3393. # to the guide. Saves time when parsing.
  3394. posToHit = {}
  3395. hitBestMismCount = {}
  3396. for fname in fnames:
  3397. for line in open(fname):
  3398. # s20+.17:A>G - chr8 26869044 CCAGCACGTGCAAGGCCGGCTTC IIIIIIIIIIIIIIIIIIIIIII 7 4:C>G,13:T>G,15:C>G
  3399. guideIdWithMod, strand, chrom, start, tSeq, weird, someScore, alnModifStr = \
  3400. line.rstrip("\n").split("\t")
  3401. guideId = guideIdWithMod.split(".")[0]
  3402. modifParts = alnModifStr.split(",")
  3403. if modifParts==['']:
  3404. modifParts = []
  3405. mismCount = len(modifParts)
  3406. hitId = (guideId, chrom, start, strand)
  3407. oldMismCount = hitBestMismCount.get(hitId, 9999)
  3408. if mismCount < oldMismCount:
  3409. hit = (mismCount, guideIdWithMod, strand, chrom, start, tSeq, modifParts)
  3410. posToHit[hitId] = hit
  3411. hitBestMismCount[hitId] = mismCount # thanks to github user mbsimonovic
  3412. ret = []
  3413. for guideId, hit in posToHit.items():
  3414. mismCount, guideIdWithMod, strand, chrom, start, tSeq, modifParts = hit
  3415. if strand=="-":
  3416. tSeq = revComp(tSeq)
  3417. guideId = guideIdWithMod.split(".")[0]
  3418. guideSeq = guideSeqs[guideId]
  3419. genomeSeq = applyModifStr(tSeq, modifParts, strand)
  3420. start = int(start)
  3421. bedRow = (guideId, chrom, start, start+GUIDELEN+pamLen, strand, guideSeq, genomeSeq)
  3422. ret.append( bedRow )
  3423. return ret
  3424. def getEditDist(str1, str2):
  3425. """ return edit distance between two strings of equal length
  3426. >>> getEditDist("HIHI", "HAHA")
  3427. 2
  3428. """
  3429. assert(len(str1)==len(str2))
  3430. str1 = str1.upper()
  3431. str2 = str2.upper()
  3432. editDist = 0
  3433. for c1, c2 in zip(str1, str2):
  3434. if c1!=c2:
  3435. editDist +=1
  3436. return editDist
  3437. def findOfftargetsBowtie(queue, batchId, batchBase, faFname, genome, pamPat, bedFname):
  3438. " align guides with pam in faFname to genome and write off-targets to bedFname "
  3439. tmpDir = batchBase+".bowtie.tmp"
  3440. os.mkdir(tmpDir)
  3441. # make sure this directory gets removed, no matter what
  3442. global tmpDirsDelExit
  3443. tmpDirsDelExit.append(tmpDir)
  3444. if not DEBUG:
  3445. atexit.register(delTmpDirs)
  3446. # write out the sequences for bowtie
  3447. queue.startStep(batchId, "seqPrep", "preparing sequences")
  3448. bwFaFname = abspath(join(tmpDir, "bowtieIn.fa"))
  3449. guideSeqs, qSeqs, allPamSeqs = writeBowtieSequences(faFname, bwFaFname, pamPat)
  3450. genomePath = abspath(join(genomesDir, genome, genome))
  3451. oldCwd = os.getcwd()
  3452. # run bowtie
  3453. queue.startStep(batchId, "bowtie", "aligning with bowtie")
  3454. os.chdir(tmpDir) # bowtie writes to hardcoded output filenames with --refout
  3455. # -v 3 = up to three mismatches
  3456. # -y = try hard
  3457. # -t = print time it took
  3458. # -k = output up to X alignments
  3459. # -m = do not output any hit if a read has more than X hits
  3460. # --max = write all reads that exceed -m to this file
  3461. # --refout = output in bowtie format, not SAM
  3462. # --maxbts=2000 maximum number of backtracks
  3463. # -p 4 = use four threads
  3464. # --mm = use mmap
  3465. maxOcc = MAXOCC # meaning in BWA: includes any PAM, in bowtie we have the PAM in the input sequence
  3466. cmd = "$BIN/bowtie -e 1000 %(genomePath)s -f %(bwFaFname)s -v 3 -y -t -k %(maxOcc)d -m %(maxOcc)d dummy --max tooManyHits.txt --mm --refout --maxbts=2000 -p 4" % locals()
  3467. runCmd(cmd)
  3468. os.chdir(oldCwd)
  3469. queue.startStep(batchId, "parse", "parsing alignments")
  3470. pamLen = len(pamPat)
  3471. hits = parseRefout(tmpDir, guideSeqs, pamLen)
  3472. queue.startStep(batchId, "scoreOts", "scoring off-targets")
  3473. # make the list of alternative PAM sequences
  3474. altPats = offtargetPams.get(pamPat, [])
  3475. altPamSeqs = []
  3476. for altPat in altPats:
  3477. altPamSeqs.extend(expandIupac(altPat))
  3478. # iterate over bowtie hits and write to a BED file with scores
  3479. # if the hit looks OK (right PAM + score is high enough)
  3480. tempBedPath = join(tmpDir, "bowtieHits.bed")
  3481. tempFh = open(tempBedPath, "w")
  3482. offTargets = {}
  3483. isSaCas9 = pamIsSaCas9(pamPat)
  3484. for guideIdWithMod, chrom, start, end, strand, _, tSeq in hits:
  3485. guideId = guideIdWithMod.split(".")[0]
  3486. guideSeq = guideSeqs[guideId]
  3487. genomePamSeq = tSeq[-pamLen:]
  3488. logging.debug( "PAM seq: %s of %s" % (genomePamSeq, tSeq))
  3489. if genomePamSeq in altPamSeqs:
  3490. minScore = ALTPAMMINSCORE
  3491. elif genomePamSeq in allPamSeqs:
  3492. minScore = MINSCORE
  3493. else:
  3494. logging.debug("Skipping off-target for %s: %s:%d-%d" % (guideId, chrom, start, end))
  3495. continue
  3496. logging.debug("off-target minScore = %f" % minScore )
  3497. # check if this match passes the off-target score limit
  3498. if pamIsCpf1(pamPat):
  3499. otScore = 0.0
  3500. else:
  3501. tSeqNoPam = tSeq[:-pamLen]
  3502. if isSaCas9:
  3503. otScore = calcSaHitScore(guideSeq, tSeqNoPam)
  3504. else:
  3505. otScore = calcHitScore(guideSeq, tSeqNoPam)
  3506. if otScore < minScore:
  3507. logging.debug("off-target not accepted")
  3508. continue
  3509. editDist = getEditDist(guideSeq, tSeqNoPam)
  3510. guideHitCount = 0
  3511. guideId = guideId.split(".")[0] # full guide ID looks like s33+.0:A>T
  3512. name = guideId+"|"+strand+"|"+str(editDist)+"|"+tSeq+"|"+str(guideHitCount)+"|"+str(otScore)
  3513. row = [chrom, str(start), str(end), name]
  3514. # this way of collecting the features will remove the duplicates
  3515. otKey = (chrom, start, end, strand, guideId)
  3516. logging.debug("off-target key is %s" % str(otKey))
  3517. offTargets[ otKey ] = row
  3518. for rowKey, row in offTargets.items():
  3519. tempFh.write("\t".join(row))
  3520. tempFh.write("\n")
  3521. tempFh.flush()
  3522. # create a tempfile which is moved over upon success
  3523. # makes sure we do not leave behind a half-written file if
  3524. # we crash later
  3525. tmpFd, tmpAnnotOffsPath = tempfile.mkstemp(dir=tmpDir, prefix="annotOfftargets")
  3526. tmpFh = open(tmpAnnotOffsPath, "w")
  3527. # get name of file with genome locus names
  3528. genomeDir = genomesDir # make local var
  3529. segFname = "%(genomeDir)s/%(genome)s/%(genome)s.segments.bed" % locals()
  3530. # annotate with genome locus names
  3531. cmd = "$BIN/overlapSelect %(segFname)s %(tempBedPath)s stdout -mergeOutput -selectFmt=bed -inFmt=bed | cut -f1,2,3,4,8 > %(tmpAnnotOffsPath)s" % locals()
  3532. runCmd(cmd)
  3533. shutil.move(tmpAnnotOffsPath, bedFname)
  3534. queue.startStep(batchId, "done", "Job completed")
  3535. if DEBUG:
  3536. logging.info("debug mode: Not deleting %s" % tmpDir)
  3537. else:
  3538. shutil.rmtree(tmpDir)
  3539. def processSubmission(faFname, genome, pamDesc, bedFname, batchBase, batchId, queue):
  3540. """ search fasta file against genome, filter for pam matches and write to bedFName
  3541. optionally write status updates to work queue. Remove faFname.
  3542. """
  3543. batchInfo = readBatchAsDict(batchId)
  3544. if genome=="noGenome":
  3545. posStr = "?"
  3546. elif "batchName" in batchInfo and batchInfo["batchName"].count(":")==2: # chrom:start-end:strand
  3547. posStr = batchInfo["batchName"]
  3548. else:
  3549. queue.startStep(batchId, "bwasw", "Searching genome for one 100% identical match to input sequence")
  3550. posStr = findPerfectMatch(batchId)
  3551. batchInfo["posStr"] = posStr
  3552. if posStr!="?":
  3553. # get a 100bp-extended version of the input seq
  3554. chrom, start, end, strand = parsePos(posStr)
  3555. extSeq = extendAndGetSeq(genome, chrom, start, end, strand, batchInfo["seq"])
  3556. if extSeq is None:
  3557. # this can only happen if there is a 100%-M match but small SNPs in it compared to the input sequence
  3558. # so the extension of the input fails.
  3559. # in this case, we also invalidate the position, as there was no perfect match and the user
  3560. # has to do something to fix it
  3561. batchInfo["posStr"] = "?"
  3562. else:
  3563. logging.debug("100pb-extended seq (len: %d) is: %s" % (len(extSeq), extSeq))
  3564. batchInfo["extSeq"] = extSeq
  3565. # must save the batch again, as otherwise display won't work, we need the position saved
  3566. writeBatchAsDict(batchInfo, batchId)
  3567. if doEffScoring:
  3568. queue.startStep(batchId, "effScores", "Calculating guide efficiency scores")
  3569. createBatchEffScoreTable(batchId, queue)
  3570. if genome=="noGenome":
  3571. # skip the off-target search entirely
  3572. open(bedFname, "w") # create a 0-byte file to signal job completion
  3573. queue.startStep(batchId, "done", "Job completed")
  3574. return
  3575. if useBowtie:
  3576. findOfftargetsBowtie(queue, batchId, batchBase, faFname, genome, pamDesc, bedFname)
  3577. else:
  3578. findOfftargetsBwa(queue, batchId, batchBase, faFname, genome, pamDesc, bedFname)
  3579. if not DEBUG:
  3580. os.remove(faFname)
  3581. return bedFname
  3582. def lineFileNext(fh):
  3583. """
  3584. parses tab-sep file with headers as field names
  3585. yields collection.namedtuples
  3586. strips "#"-prefix from header line
  3587. """
  3588. line1 = fh.readline()
  3589. while line1.startswith("##"):
  3590. line1 = fh.readline()
  3591. line1 = line1.strip("\n").strip("#")
  3592. headers = line1.split("\t")
  3593. Record = namedtuple('tsvRec', headers)
  3594. for line in fh:
  3595. line = line.rstrip("\n")
  3596. fields = line.split("\t")
  3597. try:
  3598. rec = Record(*fields)
  3599. except Exception as msg:
  3600. logging.error("Exception occured while parsing line, %s" % msg)
  3601. logging.error("Filename %s" % fh.name)
  3602. logging.error("Line was: %s" % repr(line))
  3603. logging.error("Does number of fields match headers?")
  3604. logging.error("Headers are: %s" % headers)
  3605. #raise Exception("wrong field count in line %s" % line)
  3606. continue
  3607. # convert fields to correct data type
  3608. yield rec
  3609. allGenomes = None
  3610. def readGenomes():
  3611. " return list of all genomes supported "
  3612. global allGenomes
  3613. if allGenomes:
  3614. return allGenomes
  3615. genomes = {}
  3616. myDir = dirname(__file__)
  3617. genomesDir = join(myDir, "genomes")
  3618. inFnames = []
  3619. globalFname = join(genomesDir, "genomeInfo.all.tab")
  3620. if isfile(globalFname):
  3621. inFnames = [globalFname]
  3622. else:
  3623. for subDir in os.listdir(genomesDir):
  3624. infoFname = join(genomesDir, subDir, "genomeInfo.tab")
  3625. if isfile(infoFname):
  3626. inFnames.append(infoFname)
  3627. for infoFname in inFnames:
  3628. for row in lineFileNext(open(infoFname)):
  3629. # add a note to identify UCSC genomes
  3630. if row.server.startswith("ucsc"):
  3631. addStr="UCSC "
  3632. else:
  3633. addStr = ""
  3634. genomes[row.name] = row.scientificName+" - "+row.genome+" - "+addStr+row.description
  3635. genomes = list(genomes.items())
  3636. genomes.sort(key=operator.itemgetter(1))
  3637. allGenomes = genomes
  3638. return allGenomes
  3639. def printOrgDropDown(lastorg, genomes):
  3640. " print the organism drop down box. "
  3641. print('<select id="genomeDropDown" class style="max-width:600px" name="org" tabindex="2">')
  3642. print('<option ')
  3643. if lastorg == "noGenome":
  3644. print('selected ')
  3645. print('value="noGenome">-- No Genome: no specificity, only cleavage efficiency scores (max. len 25kbp)</option>')
  3646. for db, desc in genomes:
  3647. print('<option ')
  3648. if db == lastorg :
  3649. print('selected ')
  3650. print('value="%s">%s</option>' % (db, desc))
  3651. print("</select>")
  3652. #print ('''
  3653. #<script type="text/javascript">
  3654. #$("#genomeDropDown").ufd({maxWidth:350, listWidthFixed:false});
  3655. #</script>''')
  3656. print ('''<br>''')
  3657. print ("""<script>
  3658. $("#genomeDropDown").chosen();
  3659. $(".chosen-choices li").css("background","red");
  3660. </script>
  3661. """)
  3662. def printPamDropDown(lastpam):
  3663. print('<select style="float:left" name="pam" tabindex="3">')
  3664. for key,value in pamDesc:
  3665. print('<option ')
  3666. if key == lastpam :
  3667. print('selected ')
  3668. print('value="%s">%s</option>' % (key, value))
  3669. print("</select>")
  3670. def printForm(params):
  3671. " print html input form "
  3672. scriptName = basename(__file__)
  3673. genomes = readGenomes()
  3674. haveHuman = False
  3675. for g in genomes:
  3676. if g[0]=="hg19":
  3677. haveHuman = True
  3678. # The returned cookie is available in the os.environ dictionary
  3679. cookies=http.cookies.SimpleCookie(os.environ.get('HTTP_COOKIE'))
  3680. if "lastorg" in cookies and "lastseq" in cookies and "lastpam" in cookies:
  3681. lastorg = cookies['lastorg'].value
  3682. lastseq = cookies['lastseq'].value
  3683. lastpam = cookies['lastpam'].value
  3684. else:
  3685. if not haveHuman:
  3686. global DEFAULTSEQ
  3687. global DEFAULTORG
  3688. DEFAULTSEQ = ALTSEQ
  3689. DEFAULTORG = ALTORG
  3690. lastorg = DEFAULTORG
  3691. lastseq = DEFAULTSEQ
  3692. lastpam = DEFAULTPAM
  3693. # SerialCloner is sending us the sequence via a HTTP get parameter
  3694. if "seq" in params:
  3695. lastseq = params["seq"]
  3696. if "org" in params:
  3697. lastorg = params["org"]
  3698. seqName = ""
  3699. if "seqName" in params:
  3700. seqName = params["seqName"]
  3701. printTeforBodyStart()
  3702. #print('''March 6 2023: Sorry, no CRISPOR on the new UCSC-based-server (with RS3 scores) today. Too many performance problems on the new server. We were able to renew the old server. Please use the <a href="http://37.187.154.234">old server</a> temporarily.''')
  3703. #sys.exit(0)
  3704. print("""
  3705. <form id="main-form" method="post" action="%s">
  3706. <div style="text-align:left; margin-left: 10px">
  3707. CRISPOR (<a href="https://academic.oup.com/nar/article/46/W1/W242/4995687">citation</a>) is a program that helps design, evaluate and clone guide sequences for the CRISPR/Cas9 system. <a target=_blank href="/manual/">CRISPOR Manual</a>
  3708. <br><i>July 2025: Added hasCas12Max and e-SpotOn. Also allowing old primer links to Crispor to work again. See <a href="doc/changes.html">Full list of changes</a></i><br>
  3709. </div>
  3710. <div class="windowstep subpanel" style="width:40%%;">
  3711. <div class="substep">
  3712. <div class="title">
  3713. Step 1
  3714. </div>
  3715. Planning a lentiviral gene knockout screen? Use <a href="crispor.py?libDesign=1">CRISPOR Batch</a><br>
  3716. Sequence name (optional): <input type="text" name="name" size="20" value="%s"><br>
  3717. Enter a single genomic sequence, &lt; %d bp, typically an exon
  3718. <img src="%simage/info-small.png" title="CRISPOR conserves the lowercase and uppercase format of your sequence, allowing to highlight sequence features of interest such as ATG or STOP codons.<br>Avoid using cDNA sequences as input, CRISPR guides that straddle splice sites are unlikely to work.<br>You can paste a single >23bp sequence and even multiple sequences, separated by N characters." class="tooltipster">
  3719. <br>
  3720. <small><a href="javascript:clearInput()">Clear Box</a> - </small>
  3721. <small><a href="javascript:resetToExample()">Reset to default</a></small>
  3722. </div>
  3723. <textarea tabindex="1" style="width:100%%" name="seq" rows="12"
  3724. placeholder="Paste here the genomic - not a cDNA - sequence of the exon you want to target. The sequence has to include the PAM site for your enzyme of interest, e.g. NGG. Maximum size %d bp. If you only have a cDNA, please BLAST or BLAT the cDNA first to find the right exon sequence for CRISPOR.">%s</textarea>
  3725. <small>Text case is preserved, e.g. you can mark ATGs with lowercase.<br>Instead of a sequence, you can paste a chromosome range, e.g. chr1:11,130,540-11,130,751</small>
  3726. </div>
  3727. <div class="windowstep subpanel" style="width:50%%">
  3728. <div class="substep" style="margin-bottom: 1px">
  3729. <div class="title" style="cursor:pointer;" onclick="$('#helpstep2').toggle('fast')">
  3730. Step 2
  3731. </div>
  3732. Select a genome
  3733. </div>
  3734. """% (scriptName, seqName, MAXSEQLEN, HTMLPREFIX, MAXSEQLEN, lastseq))
  3735. printOrgDropDown(lastorg, genomes)
  3736. print("""
  3737. <div id="trackHubNote" style="margin-bottom:5px">
  3738. <small>Note: pre-calculated exonic guides for this species are on the <a id='hgTracksLink' target=_blank href="">UCSC Genome Browser</a>.</small>
  3739. </div>
  3740. """)
  3741. print('<small style="float:left">We have %d genomes, but not yours? Search <a href="https://www.ncbi.nlm.nih.gov/assembly">NCBI assembly</a> and send a GCF_/GCA_ ID to <a href="mailto:%s">CRISPOR support</a>.</small>' % (len(genomes), contactEmail))
  3742. print("""
  3743. </div>
  3744. <div class="windowstep subpanel" style="width:50%%; height:158px">
  3745. <div class="substep">
  3746. <div class="title" style="cursor:pointer;" onclick="$('#helpstep3').toggle('fast')">
  3747. Step 3
  3748. <img src="%simage/info-small.png" title="The most common system uses the NGG PAM recognized by Cas9 from S. <i>pyogenes</i>. The VRER and VQR mutants were described by <a href='http://www.nature.com/nature/journal/vaop/ncurrent/abs/nature14592.html' target='_blank'>Kleinstiver et al</a>, Cas9-HF1 by <a href='https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4851738/'>Kleinstiver 2016</a>, eSpCas1.1 by <a href='https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4714946/'>Slaymaker 2016</a>, Cpf1 by <a href='http://www.cell.com/abstract/S0092-8674(15)01200-3'>Zetsche 2015</a>, SaCas9 by <a href='https://www.ncbi.nlm.nih.gov/pmc/articles/pmid/25830891/'>Ran 2015</a> and KKH-SaCas9 by <a href='https://www.ncbi.nlm.nih.gov/pmc/articles/pmid/26524662/'>Kleinstiver 2015</a>, modified As-Cpf1s by <a href='http://biorxiv.org/content/early/2016/12/04/091611'>Gao et al. 2017</a>." class="tooltipsterInteract">
  3749. </div>
  3750. Select a Protospacer Adjacent Motif (PAM)
  3751. </div>
  3752. """ % HTMLPREFIX)
  3753. printPamDropDown(lastpam)
  3754. print("""<br>See <a target=_blank href="manual/manual.html#enzymes">notes on enzymes</a> in the manual.<br>""")
  3755. print("""
  3756. <div style="width:40%; margin-top: 10px; margin-left:50px; text-align:center; display:block">
  3757. <input type="submit" name="submit" value="SUBMIT" tabindex="4"/>
  3758. </div>
  3759. </div>
  3760. """)
  3761. print("""
  3762. <script>
  3763. /* set the dropbox to hg19 and paste the example sequence into the input box. */
  3764. function resetToExample() {
  3765. $("textarea[name='seq']").val("%s");
  3766. $("#genomeDropDown").val("%s");
  3767. $("select[name='pam']").val("NGG");
  3768. }
  3769. /* clear the sequence input box */
  3770. function clearInput() {
  3771. $("textarea[name='seq']").val("");
  3772. }
  3773. </script>
  3774. <script>
  3775. /* hide the track hub note if genome is not hg19 */
  3776. ucscTrackDbs=['hg19', 'hg38', 'rn5', 'mm10', 'mm9', 'ci2', 'danRer7', 'sacCer3', 'dm6'];
  3777. function showHideHubNote() {
  3778. var valSel = $("#genomeDropDown").val();
  3779. if (jQuery.inArray(valSel, ucscTrackDbs)!=-1)
  3780. {
  3781. $("#trackHubNote").css('visibility', 'visible');
  3782. $("#hgTracksLink").attr("href", "http://genome.ucsc.edu/cgi-bin/hgTracks?db="+valSel+"&crispr=show");
  3783. }
  3784. else
  3785. $("#trackHubNote").css('visibility', 'hidden');
  3786. }
  3787. $("#genomeDropDown").on('change', showHideHubNote);
  3788. showHideHubNote();
  3789. </script>
  3790. </form>
  3791. """ % (DEFAULTSEQ, DEFAULTORG))
  3792. def readBatchAsDict(batchId):
  3793. " return contents of batch as a dictionary or None "
  3794. batchBase = join(batchDir, batchId)
  3795. jsonFname = batchBase+".json"
  3796. if isfile(jsonFname):
  3797. params = json.load(open(jsonFname))
  3798. else:
  3799. db = sqlite3.connect(batchArchive)
  3800. c = db.cursor()
  3801. c.execute("select data from jobArchive where id=?", (batchId,))
  3802. data = None
  3803. for row in c.fetchall():
  3804. data = row[0]
  3805. db.close()
  3806. if data is None:
  3807. return None
  3808. jsonStr = gzip.decompress(data).decode("utf8")
  3809. params = json.loads(jsonStr)
  3810. if "batchName" in params:
  3811. global batchName
  3812. batchName = params["batchName"]
  3813. return params
  3814. def writeBatchAsDict(batchInfo, batchId):
  3815. batchBase = join(batchDir, batchId)
  3816. tmpFname = batchBase+".json.tmp"
  3817. ofh = open(tmpFname, "w")
  3818. json.dump(batchInfo, ofh)
  3819. ofh.close()
  3820. jsonFname = batchBase+".json"
  3821. os.rename(tmpFname, jsonFname)
  3822. logging.debug("Wrote batch info to %s: %s" % (jsonFname, batchInfo))
  3823. def readBatchParams(batchId):
  3824. """ given a batchId, return the genome, the pam, the input sequence and the
  3825. chrom pos and extSeq, a 100bp-extended version of the input sequence.
  3826. Returns None for pos if not found. """
  3827. params = readBatchAsDict(batchId)
  3828. if params != None:
  3829. return params["seq"], params["org"], params["pam"], params.get("posStr"), params.get("extSeq")
  3830. # FROM HERE UP TO END OF FUNCTION: legacy cold for old batches pre-end-2016 (no json files back then)
  3831. # remove in 2017
  3832. batchBase = join(batchDir, batchId)
  3833. inputFaFname = batchBase+".input.fa"
  3834. if not isfile(inputFaFname):
  3835. errAbort('Could not find the batch %s. We cannot keep Crispor runs for more than '
  3836. 'a few months. Please resubmit your input sequence via'
  3837. ' <a href="crispor.py">the query input form</a>' % batchId)
  3838. ifh = open(inputFaFname, encoding="utf8")
  3839. ifhFields = ifh.readline().replace(">","").strip().split()
  3840. if len(ifhFields)==2:
  3841. genome, pamSeq = ifhFields
  3842. position = None
  3843. else:
  3844. genome, pamSeq, position = ifhFields
  3845. inSeq = ifh.readline().strip()
  3846. ifh.seek(0)
  3847. seqs = parseFasta(ifh)
  3848. ifh.close()
  3849. extSeq = None
  3850. if "extSeq" in seqs:
  3851. extSeq = seqs["extSeq"]
  3852. return inSeq, genome, pamSeq, position, extSeq
  3853. def gzipStr(s):
  3854. " compress a string with gzip and return "
  3855. out = StringIO()
  3856. with gzip.GzipFile(fileobj=out, mode="w") as f:
  3857. f.write(s)
  3858. return out.getvalue()
  3859. def gunzipStr(s):
  3860. " uncompress a string with gzip and return "
  3861. print(len(s), type(s), dir(s))
  3862. f = gzip.GzipFile(StringIO(s))
  3863. result = f.read()
  3864. f.close()
  3865. return result
  3866. def openDbm(dbFname, mode):
  3867. " some distributions don't include the dbm module anymore "
  3868. #import dbm.ndbm
  3869. #dbMod = dbm
  3870. #import dbm.gnu
  3871. #dbMod = gdbm
  3872. #import semidbm
  3873. # lmdbm is faster than everything else: https://pypi.org/project/lmdbm/
  3874. # though semidbm is not bad either
  3875. # Also see leveldb Wiki page
  3876. #import lmdb
  3877. #import dbm.dumb
  3878. from lmdbm import Lmdb
  3879. db = Lmdb.open(dbFname+".lmdb", mode)
  3880. if mode=="c":
  3881. # hack, I don't know how to set the permissions on the open call
  3882. cmd = ["chmod", "-R", "a+rw",dbFname+".lmdb"]
  3883. runCmd(cmd, useShell=False)
  3884. return db
  3885. def saveOutcomeData(batchId, data):
  3886. """ save outcome data of batch. data is a dictionary with key = score name """
  3887. batchBase = join(batchDir, batchId)
  3888. dbFname = batchBase
  3889. db = openDbm(dbFname, "c")
  3890. #conn = sqlite3.connect(dbFname, "w")
  3891. #c = conn.cursor()
  3892. #c.execute('''CREATE TABLE outcomes (id text PRIMARY KEY, data blob))''' % scoreName)
  3893. #c.commit()
  3894. for scoreName, data in data.items():
  3895. #c.execute("INSERT INTO outcomes values (?, ?)", (scoreName, gzipStr(json.dumps(data))))
  3896. db[scoreName] = zlib.compress(json.dumps(data).encode("utf8"))
  3897. db.close()
  3898. #c.commit()
  3899. def readOutcomeData(batchId, scoreName):
  3900. """ open outcome data of batch, key is score name """
  3901. batchBase = join(batchDir, batchId)
  3902. #conn = sqlite3.connect(dbFname, "r")
  3903. #c = conn.cursor()
  3904. #binData = c.execute("SELECT data FROM outcomes where id=?", scoreName)
  3905. #try:
  3906. # import dbm.ndbm
  3907. # db = dbm.ndbm.open(batchBase, "r") # dbm always adds .db to the file name
  3908. #except:
  3909. # # old batches on crispor.org are still using gdbm
  3910. #dbFname = batchBase+".dbm"
  3911. # import dbm.gnu
  3912. # db = dbm.gnu.open(dbFname, "r")
  3913. #import dbm.dumb
  3914. #db = dbm.dumb.open(dbFname, "r") # dbm always adds .db to the file name
  3915. db = openDbm(batchBase, "r")
  3916. dbObj = db[scoreName]
  3917. jsonStr = zlib.decompress(dbObj)
  3918. data = json.loads(jsonStr)
  3919. db.close()
  3920. return data
  3921. def findAllPams(seq, pam):
  3922. """ find all matches for PAM and return as dict startPos -> strand and a set
  3923. of end positions. The start positions for the negative strand are for the
  3924. rev-complemented PAM
  3925. """
  3926. seq = seq.upper()
  3927. startDict, endSet = findPams(seq, pam, "+", {}, set())
  3928. startDict, endSet = findPams(seq, revComp(pam), "-", startDict, endSet)
  3929. if pam in multiPams:
  3930. for pam2 in multiPams[pam]:
  3931. startDict, endSet = findPams(seq, pam2, "+", startDict, endSet)
  3932. startDict, endSet = findPams(seq, revComp(pam2), "-", startDict, endSet)
  3933. return startDict, endSet
  3934. def newBatch(batchName, seq, org, pam):
  3935. """ obtain a batch ID and write seq/org/pam to their files.
  3936. Return batchId.
  3937. """
  3938. batchId = makeTempBase(seq, org, pam, batchName)
  3939. batchData = {}
  3940. batchData["org"] = org
  3941. batchData["pam"] = pam
  3942. batchData["batchName"] = batchName
  3943. batchData["seq"] = seq
  3944. batchData["posStr"] = ""
  3945. writeBatchAsDict(batchData, batchId)
  3946. return batchId
  3947. def readDbInfo(org):
  3948. " return a dbInfo object with the columsn in the genomeInfo.tab file "
  3949. myDir = dirname(__file__)
  3950. genomesDir = join(myDir, "genomes")
  3951. infoFname = join(genomesDir, org, "genomeInfo.tab")
  3952. if not isfile(infoFname):
  3953. return None
  3954. dbInfo = next(lineFileNext(open(infoFname)))
  3955. return dbInfo
  3956. def printQueryNotFoundNote(dbInfo):
  3957. print("<div class='title'>Query sequence, not found in the selected genome, %s (%s)</div>" % (dbInfo.scientificName, dbInfo.name))
  3958. print("<div class='substep' style='border: 1px black solid; padding:5px; background-color: aliceblue'>")
  3959. print("<strong>Warning:</strong> The query sequence was not found in the selected genome.")
  3960. print("This can be a valid query, e.g. a GFP sequence.<br>")
  3961. print("If not, you might want to check if you selected the right genome for your query sequence.<br>")
  3962. print("Use a tool like <a target=_blank href='http://genome.ucsc.edu/cgi-bin/hgBlat'>BLAT</a> to check if the " \
  3963. "sequence really has a 100% identical match in the target genome.<p>")
  3964. print("When reading the list of guide sequences and off-targets below, bear in mind that in case that the input sequence is really in the genome and just has a few differences, the software will use the first found match as the on-target as it cannot distinguish 0-mismatch off-targets from 0-mismatch on-targets. In this case, the specificity scores of guide sequences are too low. In other words, some guides may be fine, the problem may just be that the on-target is shown as an off-target. <br>")
  3965. print("Because there is no flanking sequence available, the guides in your sequence that are within 50bp of the ends will have no efficiency scores. The efficiency scores will instead be shown as '--'. Include more flanking sequence > 50bp to obtain these scores.")
  3966. print("</div>")
  3967. def getOfftargets(seq, org, pamDesc, batchId, startDict, queue):
  3968. """ write guides to fasta and run bwa or use cached results.
  3969. Return name of the BED file with the matches or None if not yet available.
  3970. Write progress status updates to queue object.
  3971. """
  3972. pam = setupPamInfo(pamDesc)
  3973. assert('-' not in pam)
  3974. batchBase = join(batchDir, batchId)
  3975. otBedFname = batchBase+".bed.gz"
  3976. batchInfo = readBatchAsDict(batchId)
  3977. flagFile = batchBase+".running"
  3978. if isfile(flagFile):
  3979. errAbort("This sequence is still being processed. Please wait for ~20 seconds "
  3980. "and try again, e.g. by reloading this page. If you see this message for "
  3981. "more than 2-3 minutes, please send an email to %s. Thanks!" % contactEmail)
  3982. if not batchInfo or not isfile(otBedFname) or commandLineMode or not "posStr" in batchInfo or \
  3983. (batchInfo["posStr"]=="" and not batchInfo["org"]=="noGenome"): # pre-4.8 batches don't have a posStr at all
  3984. # write potential PAM sites to file
  3985. faFname = batchBase+".fa"
  3986. writePamFlank(seq, startDict, pam, faFname)
  3987. if commandLineMode:
  3988. processSubmission(faFname, org, pamDesc, otBedFname, batchBase, batchId, queue)
  3989. else:
  3990. q = JobQueue()
  3991. q.openSqlite()
  3992. ip = os.environ.get("REMOTE_ADDR", "noIp")
  3993. if ip=="195.176.112.240":
  3994. errAbort("IP address blocked.")
  3995. wasOk = q.addJob("search", batchId, "ip=%s,org=%s,pam=%s" % (ip, org, pamDesc))
  3996. if not wasOk:
  3997. print("CRISPOR job %s failed-running..." % batchId)
  3998. pass
  3999. q.close()
  4000. return None
  4001. return otBedFname
  4002. def showPamWarning(pam):
  4003. if pamIsCpf1(pam):
  4004. print('<div style="text-align:left; border: 1px solid; background-color: aliceblue; padding: 3px">')
  4005. print("<strong>Note:</strong> You are using the Cpf1 enzyme or related enzyme.")
  4006. print("While there is an efficiency score specificially for Cpf1, there is no off-target ranking algorithm available in the literature, to our knowledge. We use Hsu and CFD scores below for off-target ranking, but they were developed for spCas9. There is not enough data yet to support their usefulness for Cpf1. Contact us for more info if you need to rank Cpf1 off-targets for validation or if you have a dataset that could elucidate this question. We are showing out-of-frame scores, but they are based on micro-homology that assumes a spCas9 cut site, so most likely the out-of-frame scores are not accurate for the staggered cut of Cpf1 either.")
  4007. print('</div>')
  4008. #elif pamIsSaCas9(pam):
  4009. #print '<div style="text-align:left; border: 1px solid; background-color: aliceblue; padding: 3px">'
  4010. #print "<strong>Note:</strong> Your query is using a Cas9 from S. aureus.<br>"
  4011. #print "Please note that while the efficiency scoring was built for saCas9, the off-target ranking below and specificity scores are based on CFD/Hsu models, which were developed for spCas9. The ranking of off-targets could be very inaccurate. If you have a saCas9 off-target dataset, you can contact us for further info, we are only aware of the BLESS dataset by <a href='https://www.nature.com/articles/nature14299' target=_blank>Ran et al. 2015</a>.<br>As for out-of-frame and micro-homology, this model is also based on spCas9, but <a target=_blank href='https://www.nature.com/articles/nature14299'>Ran et al 2015</a> showed that the saCas9 cleavage pattern looks identical to spCas9's, so the OOF micro-homology model should work with saCas9."
  4012. #print '</div>'
  4013. elif not pamIsSpCas9(pam) and not pamIsSaCas9(pam):
  4014. print('<div style="text-align:left; border: 1px solid; background-color: aliceblue; padding: 3px">')
  4015. print("<strong>Warning:</strong> Your query involves a Cas9 that is not from S. Pyogenes and is also not Cpf1 nor saCas9.")
  4016. print("Please bear in mind that specificity and efficiency scores were designed using data with S. Pyogenes Cas9 and will very likely not be applicable to this particular Cas9. There is nothing we can do about this, we are unaware of a published dataset for this enzyme. If you know one, please contact us. Also contact us if you think another one of the existing scoring model would be more appropriate for this enzyme.<br>")
  4017. print('</div>')
  4018. if pam=="NNNNACA":
  4019. printNote("You selected the old version of the CjCas9 PAM. You may want to select the more recent "+
  4020. "PAMs from the menu on the first page, based on the study by "+
  4021. "<a target=_blank href='https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5473640/'>Kim et al 2017</a>.")
  4022. if pam=="NGN":
  4023. printNote("You have selected the NGN pam for xCas9. While this PAM is documented to work, "+
  4024. "if you read the paper in detail, you will notice that the editing efficiency is much lower. "+
  4025. "For optimal efficiency, consider going back and switching to the 'high-efficiency' xCas9 PAM.")
  4026. if pam=="NGK":
  4027. printNote("You have selected the most efficient PAM for xCas9. You can also select the more general/flexible"+
  4028. " NGN PAM from the menu when you submit your job. If you read the xCas9 paper in detail, you will find "+
  4029. "that NGN is not as efficient though.")
  4030. if pam=="TTTN":
  4031. printWarning("You selected TTTN as the PAM for Cpf1. " +
  4032. "This is not the best PAM. The actual PAM is TTTV, as shown in Fig. 2a of " +
  4033. "<a href='https://www.ncbi.nlm.nih.gov/pubmed/27992409'>Kim HK et al. Nat Meth 2017</a>.<br>")
  4034. def showNoGe

crispor.py at commit 3aa198f, under other · at the source

Overview

Authors: Florine Roussange1, Jacqueline Gide2, Johana Tournois3, Michel Cailleret1, Anne Boland4, Christophe Battail4,5, Jean-François Deleuze4, Hélène Polvèche3, Didier Auboeuf6, Knut Brockmann7, Edor Kabashi8, Anca Marian8, Lina El Kassar3, Sophie Blondel2, François Salachas9, Gaëlle Bruneteau9,10, Marc Peschanski2, Cécile Martinat1,2, Sandrine Baghdoyan1
  1. Université Paris Saclay, Université d’Evry, Inserm, IStem, UMR861, 91100 Corbeil-Essonnes, France
  2. IStem, CECS, 91100 Corbeil-Essonnes, France
  3. IStem, CECS, the research and innovation team, 91100 Corbeil-Essonnes, France
  4. Université Paris-Saclay, CEA, Centre National de Recherche en Génomique Humaine (CNRGH), 91057 Evry, France
  5. Université Grenoble Alpes, Inserm, CEA, UA13, BGE, 38000 Grenoble, France
  6. Ecole Normale Supérieure de Lyon, Inserm, U1293, CNRS, UMR 5239, Université Claude Bernard Lyon 1, Laboratory of Biology and Modelling of the Cell, 46 allée d'Italie 69364 Lyon, France
  7. Department of Pediatrics and Adolescent Medicine, University Medical Center, Göttingen, Germany
  8. Laboratory of Translational Research for Neurological Disorders, Imagine Institute, Université de Paris, INSERM, UMR 1163, 75015 Paris, France
  9. Sorbonne Université, Institut du Cerveau - Paris Brain Institute - ICM, APHP, Inserm, CNRS, Département de Neurologie, Centre SLA de Paris, Hôpital Pitié-Salpêtrière, 75013 Paris, France
  10. Alliance on Clinical Trials for ALS-MND (ACT4ALS-MND), Neuroscience Clinical Investigation Center, Paris Brain Institute, 75013 Paris, France
Journal: Stem cell reports, volume 21, issue 7, article 102977
Dates: received 25 November 2025; accepted 28 May 2026; published online 25 June 2026; in print July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1016/j.stemcr.2026.102977 · PMID 42349423 · PMCID PMC13385447 · OpenAlex W7165873553
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), human (organism), zebrafish (organism), other condition (population), cellular / molecular (subfield)
Methods: Statistics
Keywords: pluripotent stem cells, drug repositioning, amyotrophic lateral sclerosis, motor neurons, autophagy, SQSTM1, repurposing, personalized therapies, gene expression mapping
MeSH: Amyotrophic Lateral Sclerosis*, Pluripotent Stem Cells*, Sequestosome-1 Protein*, Animals, Cell Line, Fibroblasts, Gene Expression Profiling, Gene Expression Regulation, Humans, Induced Pluripotent Stem Cells, Motor Neurons, Zebrafish (* major topic)
Topic: Amyotrophic Lateral Sclerosis Research (Neurology, Medicine), according to OpenAlex
Funding: European Organisation for Rare Diseases; Association française contre les myopathies; French National Research Agency
Citations: cited by 1 paper (Europe PMC); 54 references in the paper
Research resources: Rabbit monoclonal anti-NANOG RRID:AB_10559205, SSEA3 RRID:AB_10564070, Rabbit monoclonal anti-ACTB RRID:AB_1850027, Goat polyclonal anti-ISLET1 RRID:AB_2126324, TUJ1/TUBB3 RRID:AB_2313773, anti-LC3B antibody RRID:AB_3750600, Mouse monoclonal anti-LAMP2 RRID:AB_470709, Mouse monoclonal anti-Oct3/4 RRID:AB_628051, Rabbit polyclonal anti-LC3B RRID:AB_669581, Rabbit polyclonal anti-PLP RRID:AB_776593, Mouse monoclonal anti-SQSTM1 RRID:AB_945626

Abstract

The classical paradigm of drug screening often faces significant limitations due to the challenges associated with identifying molecular or cellular read-outs that are relevant to specific genetic diseases. To remedy this, an alternative approach of reverse phenotypic mapping was tested: Compounds were evaluated for their effects on gene expression and alternative splicing in a healthy cell model, and the resulting data were matched to molecular signatures of diseases. A subset of 50 drugs was tested on mesenchymal stem cells derived from a human pluripotent stem cell line. Over half of the compounds altered gene expression, many affecting pathways linked to monogenic diseases. One hit, increased SQSTM1 expression induced by prazosin, was further validated in FTD/ALS type 3 models caused by SQSTM1 haploinsufficiency, including patient-derived fibroblasts, SQSTM1-depleted hiPSC-derived motor neurons, and a zebrafish model. Extending this paradigm could involve testing diverse cell types and larger drug libraries.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 2 matches between paragraphs and lines of code.

maximilianh/crisporWebsite

License: other
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 3aa198f1371bb15733acd98aaa097aad052cdbe3, 26 September 2026
Languages: C/C++ (278), C (275), JavaScript (98), Python (86), C++ (67), Perl (33), Java (8), Shell (8), R (4), Stan (1)
Size: 3,120 files, 858 scripts
Software Heritage: not archived
Found in: the resources table
Holds: README, license file, environment (requirements.txt, docker/Dockerfile, bin/Azimuth-2.0/setup.cfg, bin/Azimuth-2.0/setup.py, bin/src/lindel/setup.py, bin/src/ViennaRNA-2.1.9/interfaces/Python/setup.py), tests, documentation
Not found: CITATION.cff, continuous integration
Tools: NumPy (45 files), pandas (24 files), SciPy (19 files), scikit-learn (15 files), Biopython (11 files), Matplotlib (11 files), Keras (2 files), Stan (2 files), limma (1 file), reshape2 (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
860 files

usadellab/trimmomatic

License: GPL-3.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: ef98d6252abeae80cfee36acf9e0e1055097da0b, 3 July 2026
Languages: Java (150)
Size: 178 files, 150 scripts
Software Heritage: not archived
Found in: the resources table
Holds: README, license file, environment (Dockerfile), tests, continuous integration
Not found: CITATION.cff, documentation
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
152 files

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 1,008 scripts, each with its path and the digest of its content;
  • 2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data and code availability

Bulk RNA-seq data have been deposited at GEO: GSE261648 (https://ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE261648) accession number and are publicly available as of the date of publication.

Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 2, 28 September 2026

  • Authors: added Cécile Martinat (0000-0002-5234-1064); Sandrine Baghdoyan (0000-0002-5904-7394); removed Cécile Martinat; Sandrine Baghdoyan

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 19 authors, 9 keywords, 12 MeSH terms, 3 funders, 54 references, 11 RRIDs.

Cite

This paper

Roussange, F., Gide, J., Tournois, J., Cailleret, M., Boland, A., Battail, C., Deleuze, J.-F., Polvèche, H., Auboeuf, D., Brockmann, K., Kabashi, E., Marian, A., El Kassar, L., Blondel, S., Salachas, F., Bruneteau, G., Peschanski, M., Martinat, C., & Baghdoyan, S. (2026). Integrative analysis of drug-gene signatures in human pluripotent stem cells reveals prazosin as a novel SQSTM1 regulator for ALS therapeutics. Stem cell reports, 21(7), 102977. https://doi.org/10.1016/j.stemcr.2026.102977

BibTeX

@article{roussange2026integrative,
author = {Roussange, Florine and Gide, Jacqueline and Tournois, Johana and Cailleret, Michel and Boland, Anne and Battail, Christophe and Deleuze, Jean-François and Polvèche, Hélène and Auboeuf, Didier and Brockmann, Knut and Kabashi, Edor and Marian, Anca and El Kassar, Lina and Blondel, Sophie and Salachas, François and Bruneteau, Gaëlle and Peschanski, Marc and Martinat, Cécile and Baghdoyan, Sandrine},
title = {{Integrative analysis of drug-gene signatures in human pluripotent stem cells reveals prazosin as a novel SQSTM1 regulator for ALS therapeutics}},
journal = {Stem cell reports},
year = {2026},
month = jun,
volume = {21},
number = {7},
pages = {102977},
publisher = {Elsevier},
issn = {2213-6711},
doi = {10.1016/j.stemcr.2026.102977},
url = {https://doi.org/10.1016/j.stemcr.2026.102977},
pmid = {42349423},
pmcid = {PMC13385447}
}

RIS

TY - JOUR
AU - Roussange, Florine
AU - Gide, Jacqueline
AU - Tournois, Johana
AU - Cailleret, Michel
AU - Boland, Anne
AU - Battail, Christophe
AU - Deleuze, Jean-François
AU - Polvèche, Hélène
AU - Auboeuf, Didier
AU - Brockmann, Knut
AU - Kabashi, Edor
AU - Marian, Anca
AU - El Kassar, Lina
AU - Blondel, Sophie
AU - Salachas, François
AU - Bruneteau, Gaëlle
AU - Peschanski, Marc
AU - Martinat, Cécile
AU - Baghdoyan, Sandrine
TI - Integrative analysis of drug-gene signatures in human pluripotent stem cells reveals prazosin as a novel SQSTM1 regulator for ALS therapeutics
T2 - Stem cell reports
J2 - Stem Cell Reports
PY - 2026
DA - 2026/06/25
VL - 21
IS - 7
SP - 102977
SN - 2213-6711
PB - Elsevier
DO - 10.1016/j.stemcr.2026.102977
UR - https://doi.org/10.1016/j.stemcr.2026.102977
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.stemcr.2026.102977",
"type": "article-journal",
"title": "Integrative analysis of drug-gene signatures in human pluripotent stem cells reveals prazosin as a novel SQSTM1 regulator for ALS therapeutics",
"container-title": "Stem cell reports",
"author": [
{
"family": "Roussange",
"given": "Florine"
},
{
"family": "Gide",
"given": "Jacqueline"
},
{
"family": "Tournois",
"given": "Johana"
},
{
"family": "Cailleret",
"given": "Michel"
},
{
"family": "Boland",
"given": "Anne"
},
{
"family": "Battail",
"given": "Christophe"
},
{
"family": "Deleuze",
"given": "Jean-François"
},
{
"family": "Polvèche",
"given": "Hélène"
},
{
"family": "Auboeuf",
"given": "Didier"
},
{
"family": "Brockmann",
"given": "Knut"
},
{
"family": "Kabashi",
"given": "Edor"
},
{
"family": "Marian",
"given": "Anca"
},
{
"family": "El Kassar",
"given": "Lina"
},
{
"family": "Blondel",
"given": "Sophie"
},
{
"family": "Salachas",
"given": "François"
},
{
"family": "Bruneteau",
"given": "Gaëlle"
},
{
"family": "Peschanski",
"given": "Marc"
},
{
"family": "Martinat",
"given": "Cécile"
},
{
"family": "Baghdoyan",
"given": "Sandrine"
}
],
"container-title-short": "Stem Cell Reports",
"volume": "21",
"issue": "7",
"page": "102977",
"DOI": "10.1016/j.stemcr.2026.102977",
"PMID": "42349423",
"PMCID": "PMC13385447",
"ISSN": "2213-6711",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.stemcr.2026.102977",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
25
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: Stan, Biopython, limma, 7 other tools, zebrafish, cellular / molecular
[2] doi:10.1038/s41467-026-76675-1 [code]
Long-read proteogenomic atlas of human neuronal differentiation reveals isoform diversity informing neurodevelopmental risk mechanisms.
Journal: Nature communications
In common: Biopython, limma, reshape2, 5 other tools, genetics / omics, 2 references
[3] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: limma, Keras, reshape2, 5 other tools, genetics / omics, other condition, cellular / molecular, 1 reference
[4] doi:10.1038/s44318-026-00818-9 [code]
FAM134B-mediated ER-phagy degrades APP and suppresses Alzheimer's disease pathology.
Journal: The EMBO journal
In common: Biopython, limma, reshape2, 5 other tools, cellular / molecular
[5] doi:10.1186/s13293-026-00927-4 [code]
Gene regulatory network analysis identifies dysregulation of hypoxia pathways as contributing to glioblastoma treatment resistance in females.
Journal: Biology of sex differences
In common: Stan, limma, reshape2, 4 other tools, other condition, cellular / molecular
[6] doi:10.1038/s41592-026-03057-2 [code]
CREsted: modeling genomic and synthetic cell-type-specific enhancers across tissues and species.
Journal: Nature methods
In common: Biopython, Keras, scikit-learn, 4 other tools, zebrafish, genetics / omics
[7] doi:10.1038/s43587-026-01207-x [code]
A microprotein atlas of the human frontal cortex in Alzheimer's disease.
Journal: Nature aging
In common: Biopython, limma, reshape2, 4 other tools, genetics / omics
[8] doi:10.1126/sciadv.aed2952 [code]
Activation of transposable elements is linked to a region- and cell type-specific interferon response in Parkinson's disease.
Journal: Science advances
In common: Biopython, limma, reshape2, 4 other tools, cellular / molecular
[9] doi:10.1038/s41467-026-75700-7 [code]
Gene regulatory innovations from transposable elements in primate cerebellum development.
Journal: Nature communications
In common: Biopython, Keras, scikit-learn, 4 other tools, genetics / omics, cellular / molecular
[10] doi:10.1038/s41467-026-71391-2 [code]
Accelerating Leigh syndrome drug discovery through deep learning screening in brain organoids.
Journal: Nature communications
In common: scikit-learn, pandas, SciPy, 2 other tools, 3 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.