dna/moleculardna/createSDD.py

0001 #!/usr/bin/env python
0002 # coding: utf-8
0003
0004 # In[1]:
0005
0006
0007 # To run, change the links and definitions in this script and
0008 # type in the terminal containing the file
0009 # $python3 createSDD.py
0010
0011 # Author contact, Konstantinos Chatzipapas, chatzipa@cenbg.in2p3.fr, 14/03/2022
0012
0013
0014 # In[2]:
0015
0016
0017 import ROOT as root
0018 import numpy as np
0019
0020
0021 # In[3]:
0022
0023
0024 # Define filenames
0025 outputFile = "SDD_minFormatMolecularDNA.txt"
0026 iMacFile   = "human_cell.mac"
0027 iRootFile  = "molecular-dna.root"
0028
0029
0030 # In[4]:
0031
0032
0033 # Definitions
0034 # They may be adapted in the code, if such information exists in the simulation definition file
0035
0036 # Define if it is a simulation of a single-track irradiation (0), a delivered dose (1) or a fluence (2)
0037 doseFluence = 1
0038 amountDoseFluence = 0
0039
0040 # Define dose rate
0041 doseRate = float(0)
0042
0043 # Define the number of chromosomes and the length in MBp
0044 nbChromo = 1
0045 chromoLength = 6405.886128
0046
0047 # Define DNA density (BPs/um3)
0048 dnaDensity = 1.2*1e07
0049
0050 # Define cell cycle phase (use G0 (1), G1 (2), S (3), G2 (4) and M (5), together with percentage, e.g [3,0.7] indicates a cell 70% of the way through S phase)
0051 cellCyclePhase = "0, 0.0"
0052
0053 # Define DNA structure :  whole nucleus (0), a heterochromatin region (1),
0054                         # euchromatin region (2),
0055                         # a mixed (heterochromatin and euchromatin) region (3),
0056                         # single DNA fiber (4), DNA wrapped around a single histone (5),
0057                         # DNA plasmid (6) or a simple circular (7) or straight (8) DNA section.
0058 dnaStructure = 4
0059 nakedWet = 1            # naked (0), wet (1)
0060
0061 # Define if the real experiment was in-vitro (0) or in-vivo (1)
0062 vitroVivo = 0
0063
0064 # Define the proliferation status, quiescent (0) or proliferating (1)
0065 proliferation = 0
0066
0067 # Define the microenvironment
0068 temperature = 27  # degrees
0069 oxygen = 0.0      # molarity
0070
0071 # Damage definition
0072 directIndirect = 1   # direct effects only (0) or including chemistry (1)
0073 bpNm = 0             # BPs (0) or in nm (1)
0074 bpThres = 10.0       # the distance in BPs or nm between backbone lesions that are considered DSBs
0075 baseLesions = -1     # This value then determines the distance (in BP or nm) beyond the outer backbone damages where base damages are also stored in the same site (float).
0076
0077 # Default, if nothing exists in mac file
0078 sourceTropus = "iso"
0079
0080
0081 # In[5]:
0082
0083
0084 # Initialize several parameters
0085 title    = "No Title Included"
0086 incidentParticles = " "
0087 particle = " "
0088 energy   = " "
0089 energyUnits = " "
0090 sourceShape = " "
0091 sourceTropus = " "
0092 sourceX = "0"
0093 sourceY = "0"
0094 sourceZ = "0"
0095 sourceUnits = "nm"
0096 targetType  = " "
0097 targetShape = " "
0098 targetRadius = 0
0099 targetX  = 0
0100 targetY  = 0
0101 targetZ  = 0
0102 r3 = 0
0103 worldSize = 0
0104 directDamageL = " "
0105 directDamageU = " "
0106 time = " "
0107 timeUnit = " "
0108
0109
0110
0111 # In[6]:
0112
0113
0114 with open(iMacFile, 'r') as infile:
0115     for line in infile:
0116         spl = line.split(" ")
0117         if spl[0] == "###":
0118             title = spl[1]
0119             #print(title)
0120         if spl[0] =="/run/beamOn":
0121             incidentParticles = int(spl[1])
0122             #print(incidentParticles)
0123         if spl[0] =="/gps/particle":
0124             spl2 = spl[1].split("\n")
0125             particle = spl2[0]
0126             #print(particle)
0127         if spl[0] =="/gps/energy":
0128             sourceType = 1
0129             energyDistribution = "M, 0"
0130             particleFraction = float(1)
0131             energy = float(spl[1])
0132             spl2 = spl[2].split("\n")
0133             energyUnits = spl2[0]
0134             #print(energy,energyUnits)
0135         if spl[0] =="/gps/pos/shape":
0136             spl2 = spl[1].split("\n")
0137             sourceShape = spl2[0]
0138             #print(sourceShape)
0139         if spl[0] =="/gps/ang/type":
0140             spl2 = spl[1].split("\n")
0141             sourceTropus = spl2[0]
0142             #print("type",sourceTropus)
0143         if spl[0] =="/gps/pos/centre":
0144             sourceX = spl[1]
0145             sourceY = spl[2]
0146             sourceZ = spl[3]
0147             spl2 = spl[4].split("\n")
0148             sourceUnits = spl2[0]
0149         if spl[0] =="/chromosome/add":
0150             targetType = spl[1]
0151             targetShape = spl[2]
0152             if targetShape == "sphere":
0153                 #print (spl)
0154                 if spl[4]=="nm\n":
0155                     targetRadius = float(spl[3])/1000
0156                     r3 = float(spl[3])*1e-09
0157                     #print ("r3=", r3)
0158                 elif spl[4]=="um\n":
0159                     targetRadius = float(spl[3])
0160                     r3 = float(spl[3])*1e-06
0161                     #print ("r3=", r3)
0162             if targetShape == "ellipse":
0163                 #print (spl)
0164                 if spl[9]=="nm":
0165                     targetX = float(spl[3])/1000
0166                     targetY = float(spl[4])/1000
0167                     targetZ = float(spl[5])/1000
0168                     r3 = float(spl[3])*1e-09*float(spl[4])*1e-09*float(spl[5])*1e-09
0169                     #print ("r3=", r3)
0170                 elif spl[9]=="um":
0171                     targetX = float(spl[3])
0172                     targetY = float(spl[4])
0173                     targetZ = float(spl[5])
0174                     r3 = float(spl[3])*1e-06*float(spl[4])*1e-06*float(spl[5])*1e-06
0175                     #print ("r3=", r3)
0176             #print(targetType,targetShape,targetRadius,targetX,targetY,targetZ)
0177         if spl[0] =="/world/worldSize":
0178             worldSize = float(spl[1])/2
0179             if spl[2] =="nm\n":
0180                 worldSize = worldSize/1000
0181             #print(worldSize)
0182         if spl[0] =="/dnadamage/directDamageLower":
0183             directDamageL = spl[1]
0184             #print(directDamageL)
0185         if spl[0] =="/dnadamage/directDamageUpper":
0186             directDamageU = spl[1]
0187             #print(directDamageU)
0188         if spl[0] =="/scheduler/endTime":
0189             time = spl[1]
0190             spl2 = spl[2].split("\n")
0191             timeUnit = spl2[0]
0192             #print(time, timeUnit)
0193
0194
0195 # In[7]:
0196
0197
0198 # Creating SDD file
0199 with open(outputFile, 'w') as ofile:
0200 # Starting with header part of SDD file
0201     ofile.write("\nSDD version, SDDv1.0;\n")
0202     ofile.write("Software, MolecularDNA;\n")
0203     ofile.write("Author contact, Konstantinos Chatzipapas, chatzipa@cenbg.in2p3.fr, "
0204                 "14/03/2022;\n")
0205     ofile.write("***Important information*********************************************\n"
0206                 "To provide some extra information on the quantification of DSB, "
0207                 "DSB+ and DSB++, the last column of data section, includes the values "
0208                 "of 4, 5 that correspond to DSB+ and DSB++ respectively;\n"
0209                 "*********************************************************************\n")
0210     #Adjust description
0211     ofile.write("Simulation Details, "+title+". DNA damages from direct and indirect effects;\n")
0212     ofile.write("Source, Monoenergetic "+sourceShape+" "+particle+" "+sourceTropus+"tropic,")
0213     ofile.write(" centered at "+sourceX+" "+sourceY+" "+sourceZ+" "+sourceUnits)
0214     ofile.write(" of a cell nucleus, containing a free DNA segment. Energy: ")
0215     ofile.write(str(energy)+" "+energyUnits+";\n")
0216     ######
0217     ofile.write("Source type, "+str(sourceType)+";\n")   # Needs improvement
0218     ofile.write("Incident particles, "+str(incidentParticles)+";\n")
0219     ofile.write("Mean particle energy, "+str(energy)+" MeV;\n")
0220     ofile.write("Energy distribution, "+energyDistribution+";\n")   # Needs improvement
0221     ofile.write("Particle fraction, "+str(particleFraction)+";\n")   # Needs improvement
0222
0223     f = root.TFile(iRootFile)
0224
0225     eVtoJ = 1.60218e-19
0226
0227     mass = 997 * 4 * 3.141592 * r3 / 3     # density * 4/3 * pi * r3
0228
0229     acc_edep = 0
0230     dose = 0
0231     nbEntries = 0
0232     ffTree = f.Get("tuples/chromosome_hits")
0233     nbEntries += ffTree.GetEntries()
0234     for ev in ffTree:
0235         acc_edep += (ev.e_chromosome_kev + ev.e_dna_kev) *1e3  # eV
0236
0237
0238     # Calculate the absorbed dose
0239     amountDoseFluence = acc_edep * eVtoJ / mass    # Dose in Gy
0240
0241
0242
0243     ofile.write("Dose or fluence, "+str(doseFluence)+", "+str(amountDoseFluence)+";\n")   # Needs improvement
0244     ofile.write("Dose rate, "+str(doseRate)+";\n")
0245     ofile.write("Irradiation target, Simple free DNA fragment in a "+targetShape+" with ")
0246     if targetShape == "sphere":
0247         ofile.write("radius "+str(targetRadius)+" um;\n")
0248     else: #targetShape == "ellipse":
0249         ofile.write("dimensions "+str(targetX)+" "+str(targetY)+" "+str(targetZ)+" um;\n")
0250     ofile.write("Volumes, 0,"+str(worldSize)+","+str(worldSize)+","+str(worldSize)+","+str(-worldSize)+","+str(-worldSize)+","+str(-worldSize)+",")
0251     if targetShape == "sphere":
0252         ofile.write(" 1,"+str(targetRadius)+","+str(targetRadius)+","+str(targetRadius)+";\n")
0253     elif targetShape == "ellipse":
0254         ofile.write(" 1,"+str(targetX)+","+str(targetY)+","+str(targetZ)+";\n")
0255     else:
0256         ofile.write(" 0,"+str(targetX)+","+str(targetY)+","+str(targetZ)+","+str(+targetX)+","+str(+targetY)+","+str(+targetZ)+";\n")
0257     ofile.write("Chromosome sizes, "+str(nbChromo)+", "+str(chromoLength)+";\n")
0258     ofile.write("DNA Density, "+str(dnaDensity)+";\n")
0259     ofile.write("Cell Cycle Phase, "+cellCyclePhase+";\n")
0260     ofile.write("DNA Structure, "+str(dnaStructure)+", "+str(nakedWet)+";\n")
0261     ofile.write("In vitro / in vivo, "+str(vitroVivo)+";\n")
0262     ofile.write("Proliferation status, "+str(proliferation)+";\n")
0263     ofile.write("Microenvironment, "+str(temperature)+", "+str(oxygen)+";\n")
0264     ofile.write("Damage definition, "+str(directIndirect)+", "+str(bpNm)+", "+str(bpThres)+", "+str(baseLesions)+", "+str(directDamageL)+";\n")
0265     ofile.write("Time, "+str(time)+" "+timeUnit+";\n")
0266
0267     number = 0
0268     gTree = f.Get("tuples/primary_source")
0269     number += gTree.GetEntries()
0270
0271     nEntries = 0
0272     fTree = f.Get("tuples/damage")
0273     nEntries += fTree.GetEntries()
0274     #print (nEntries)
0275
0276     ofile.write("Damage and primary count, "+str(nEntries)+", "+str(number)+";\n")
0277     ofile.write("Data entries, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0;\n")
0278     ofile.write("Data was generated in minimal output format and used as an example. ")
0279     ofile.write("Please modify description to match different conditions;\n")
0280     ofile.write("\n***EndOfHeader***;\n\n")
0281
0282 #############################################################################################
0283 # Writing the data part of SDD file
0284     currentEvent = 0
0285     Primaryflag = True
0286     Eventflag = False
0287
0288     for entry in fTree:
0289         if entry.Event!=currentEvent:
0290             Eventflag = True
0291         currentEvent = entry.Event
0292         #print(entry.SourceClassification)
0293         DSB = 0
0294         if (entry.TypeClassification == "DSB"):
0295             if (entry.SourceClassification == "DSBd"):
0296                 DSB = 1
0297             elif (entry.SourceClassification == "DSBi"):
0298                 DSB = 2
0299             elif (entry.SourceClassification == "DSBh" or entry.SourceClassification == "DSBm"):
0300                 DSB = 3
0301                 #print (entry.TypeClassification, entry.SourceClassification)
0302         elif (entry.TypeClassification == "DSB+"):
0303             DSB = 4
0304         elif (entry.TypeClassification == "DSB++"):
0305             DSB = 5
0306             #if (entry.SourceClassification == "DSBh" or entry.SourceClassification == "DSBm"):
0307                 #print (entry.TypeClassification, entry.SourceClassification)
0308         else:
0309             DSB = 0
0310
0311         #print(DSB)
0312         if Primaryflag:
0313             #print("True")
0314             ofile.write("2, "+str(entry.Event)+"; ")
0315             ofile.write(str(entry.Position_x_um)+", "+str(entry.Position_y_um)+", "+str(entry.Position_z_um)+"; ")
0316             ofile.write(str(entry.BaseDamage)+", "+str(entry.StrandDamage)+", "+str(DSB)+";\n")
0317             Primaryflag = False
0318             Eventflag = False
0319         else:
0320             if Eventflag:
0321                 ofile.write("1, "+str(entry.Event)+"; ")
0322                 ofile.write(str(entry.Position_x_um)+", "+str(entry.Position_y_um)+", "+str(entry.Position_z_um)+"; ")
0323                 ofile.write(str(entry.BaseDamage)+", "+str(entry.StrandDamage)+", "+str(DSB)+";\n")
0324                 Eventflag = False
0325             else:
0326                 ofile.write("0, "+str(entry.Event)+"; ")
0327                 ofile.write(str(entry.Position_x_um)+", "+str(entry.Position_y_um)+", "+str(entry.Position_z_um)+"; ")
0328                 ofile.write(str(entry.BaseDamage)+", "+str(entry.StrandDamage)+", "+str(DSB)+";\n")
0329
0330
0331 # In[8]:
0332
0333 print("Output file: ", outputFile)
0334
0335
0336