dna/moleculardna/human_cell_chromosomes.py

0001 # %%
0002 # To run, change the links and definitions in this script and
0003 # type in the terminal containing the file
0004 # $python3 chromosAnalysis.py
0005
0006 # python version tested 3.10.12
0007 # developed in visual studio code (jupyter notebook)
0008 # Author contact, Konstantinos Chatzipapas, konstantinos.chatzipapas@cern.ch, 08/12/2023
0009
0010 # %%
0011 import ROOT as root
0012 import numpy as np
0013 import matplotlib.pyplot as plt
0014
0015 # %%
0016 # Define filenames
0017 iRootFile = "molecular-dna.root"
0018 print("Loaded file: ",iRootFile )
0019
0020 # %%
0021 # Initial parameters (need be inserted manually, if modifications are made in the geometry)
0022 eVtoJ = 1.60218e-19
0023 r3 = 10575e-9 * 3450e-9 * 10575e-9   # a * b * c
0024 mass = 997 * 4 * 3.141592 * r3 / 3   # waterDensity * 4/3 * pi * r3 in kg
0025
0026 # %%
0027 # Length of chromosomes in base pairs (bp) (need be inserted manually, if geometry file is modified)
0028 lc =np.zeros((24,1))
0029 lc[1]  = (335601410)
0030 lc[2]  = (424230628)
0031 lc[3]  = (330358245)
0032 lc[4]  = (214431219)
0033 lc[5]  = (225844963)
0034 lc[6]  = (187075736)
0035 lc[7]  = (104619404)
0036 lc[8]  = (523147774)
0037 lc[9]  = (286308941)
0038 lc[10] = (308642731)
0039 lc[11] = (281141907)
0040 lc[12] = (286680368)
0041 lc[13] = (297671931)
0042 lc[14] = (242439583)
0043 lc[15] = (517729277)
0044 lc[16] = (363900456)
0045 lc[17] = (407437606)
0046 lc[18] = (352426730)
0047 lc[19] = (165113373)
0048 lc[20] = (137667917)
0049 lc[21] = (132362463)
0050 lc[22] = (159987865)
0051 lc[23] = (99027882)
0052
0053 # Calculation of cumulative length and total length of the cell DNA
0054 sum_lc=0
0055 slc =np.zeros((24,1))
0056 for j in range(len(lc)):
0057     sum_lc += lc[j]
0058     slc[j] = sum_lc
0059 print("Total length of DNA in the cell(bp): ", sum_lc)
0060 #print(slc) #Length of each chromosome
0061
0062 # %%
0063 # Calculate the damage yield
0064 number = 0
0065 f = root.TFile(iRootFile)
0066 gTree = f.Get("tuples/primary_source")
0067 number += gTree.GetEntries()
0068
0069 nEntries = 0
0070 fTree = f.Get("tuples/damage")
0071 nEntries += fTree.GetEntries()
0072 #print (nEntries)
0073
0074 c_DSBBPID = {}
0075 for i in range(1,len(lc)):
0076     c_DSBBPID[i] = []
0077
0078
0079 SSB  = 0
0080 ttSSB = 0
0081 ttDSB = 0
0082
0083 c_SSB =np.zeros(24)
0084 c_DSB =np.zeros(24)
0085 cc_SSB=np.zeros(24)
0086 cc_DSB=np.zeros(24)
0087
0088 c_SSBp  = np.zeros(24)
0089 c_DSBp  = np.zeros(24)
0090 c_SSBpp = np.zeros(24)
0091 c_DSBpp = np.zeros(24)
0092
0093 DSBd = np.zeros(24)
0094 DSBi = np.zeros(24)
0095 DSBm = np.zeros(24)
0096 DSBh = np.zeros(24)
0097
0098 SSBd = np.zeros(24)
0099 SSBi = np.zeros(24)
0100 SSBm = np.zeros(24)
0101
0102 # Calculate total number of SSB, DSB, as well as initial break positions for the calculation of fragments
0103 DD = 0
0104 for e in fTree:
0105     if (e.TypeClassification=="DSB" or e.TypeClassification=="DSB+" or e.TypeClassification=="DSB++"):
0106         DD += 1
0107         ttDSB += e.DirectBreaks + e.IndirectBreaks
0108         for i in range(1,len(lc)):
0109             if (e.BasePair>=slc[i-1] and e.BasePair<slc[i]):
0110                 cc_DSB[i] += 1    #e.DirectBreaks + e.IndirectBreaks
0111                 c_DSBBPID[i].append((e.Event, e.BasePair, e.TypeClassification))
0112
0113
0114     if (e.TypeClassification=="SSB" or e.TypeClassification=="SSB+" or e.TypeClassification=="2SSB"):
0115         SSB  += 1
0116         ttSSB += e.DirectBreaks + e.IndirectBreaks
0117         for i in range(len(lc)):
0118             if (e.BasePair>=slc[i-1] and e.BasePair<slc[i]):
0119                 cc_SSB[i] += 1    #e.DirectBreaks + e.IndirectBreaks
0120
0121 #print(DD,SSB)  #Used for testing
0122
0123 # Calculate damage complexity for each chromosome
0124 for e in fTree:
0125     for i in range(1,len(lc)):
0126         if (e.BasePair>=slc[i-1] and e.BasePair<slc[i]):
0127             if e.TypeClassification=="DSB":
0128                 c_DSB[i] += 1
0129             elif e.TypeClassification=="DSB+":
0130                 c_DSBp[i] += 1
0131             elif e.TypeClassification=="DSB++":
0132                 c_DSBpp[i] += 1
0133             elif e.TypeClassification=="SSB":
0134                 c_SSB[i] += 1
0135             elif e.TypeClassification=="SSB+":
0136                 c_SSBp[i] += 1
0137             elif e.TypeClassification=="2SSB":
0138                 c_SSBpp[i] += 1
0139
0140 for e in fTree:
0141     for i in range(1,len(lc)):
0142         if (e.BasePair>=slc[i-1] and e.BasePair<slc[i]):
0143             if e.SourceClassification=="DSBd":
0144                 DSBd[i] += 1
0145             elif e.SourceClassification=="DSBi":
0146                 DSBi[i] += 1
0147             elif e.SourceClassification=="DSBm":
0148                 DSBm[i] += 1
0149             elif e.SourceClassification=="DSBh":
0150                 DSBh[i] += 1
0151             elif e.SourceClassification=="SSBd":
0152                 SSBd[i] += 1
0153             elif e.SourceClassification=="SSBi":
0154                 SSBi[i] += 1
0155             elif e.SourceClassification=="SSBm":
0156                 SSBm[i] += 1
0157
0158 # Damage classification, this is not needed
0159 repDSB = 0
0160 irrDSB =0
0161 tDSB = 0
0162
0163 gTree = f.Get("tuples/classification")
0164 for e in gTree:
0165     repDSB += e.DSB
0166     irrDSB += e.DSBp + e.DSBpp
0167     tDSB += e.DSB + e.DSBp + 2*e.DSBpp
0168
0169
0170 tSSB = 0
0171 ffTree = f.Get("tuples/source")
0172 for e in ffTree:
0173     tSSB += e.SSBd + e.SSBi + e.SSBm
0174
0175 totalDSB = repDSB + irrDSB
0176
0177 #print (ttDSB,tDSB,totalDSB,repDSB,irrDSB,"  ",SSB,ttSSB,tSSB,tSSB/ttDSB, tSSB/i)  #Used for testing
0178
0179 acc_edep=0
0180 hTree= f.Get("tuples/chromosome_hits")
0181 for e in hTree:
0182     acc_edep += (e.e_chromosome_kev + e.e_dna_kev )*1e3
0183
0184 dose = acc_edep * eVtoJ / mass
0185 #print("Accumulated deposited energy in the cell(MeV): ", acc_edep)
0186 #print("Absorbed Dose to the cell(Gy): ", dose)
0187
0188 # %%
0189 #print(c_DSBBPID)
0190 #print(c_DSBBPID[22])
0191
0192 # %%
0193 # Printing damage classification for each chromosome
0194 cSSBsum   =0
0195 cSSBpsum  =0
0196 cSSBppsum =0
0197 cDSBsum   =0
0198 cDSBpsum  =0
0199 cDSBppsum =0
0200 ccSSBsum  =0
0201 ccDSBsum  =0
0202
0203
0204 for i in range(1,len(lc)):
0205     print("\nchromosome",i," total DSB/Gy/GBp: ", cc_DSB[i]/(lc[i]/1e+9)/dose," --> DSB:", c_DSB[i]/(lc[i]/1e+9)/dose, "DSB+:", c_DSBp[i]/(lc[i]/1e+9)/dose, "DSB++:", c_DSBpp[i]/(lc[i]/1e+9)/dose)
0206     print(  "chromosome",i," total SSB/Gy/GBp: ", cc_SSB[i]/(lc[i]/1e+9)/dose," --> SSB:", c_SSB[i]/(lc[i]/1e+9)/dose, "SSB+:", c_SSBp[i]/(lc[i]/1e+9)/dose, "SSB++:", c_SSBpp[i]/(lc[i]/1e+9)/dose)
0207     print(  "chromosome",i," DSBd: ", DSBd[i]/(lc[i]/1e+9)/dose,"DSBi:", DSBi[i]/(lc[i]/1e+9)/dose, "DSBm:", DSBm[i]/(lc[i]/1e+9)/dose, "DSBh:", DSBh[i]/(lc[i]/1e+9)/dose)
0208     print(  "chromosome",i," SSBd: ", SSBd[i]/(lc[i]/1e+9)/dose,"SSBi:", SSBi[i]/(lc[i]/1e+9)/dose, "DSBm:", SSBm[i]/(lc[i]/1e+9)/dose)
0209     #print("chromosome",i," total DSB/Gy/GBp: ", cc_DSB[i]/(lc[i]/1e+9)/dose, "  \tchromosome",i," total SSB/Gy/GBp: ", cc_SSB[i]/(lc[i]/1e+9)/dose)
0210     cDSBsum += c_DSB[i]+c_DSBp[i]+c_DSBpp[i]
0211     cSSBsum += c_SSB[i]+c_SSBp[i]+c_SSBpp[i]
0212     cSSBpsum  += c_SSBp[i]
0213     cSSBppsum += c_SSBpp[i]
0214     cDSBpsum  += c_DSBp[i]
0215     cDSBppsum += c_DSBpp[i]
0216
0217     ccDSBsum += cc_DSB[i]
0218     ccSSBsum += cc_SSB[i]
0219 #print("total DSB/Gy/GBp:",ccDSBsum/(sum_lc/1e+9)/dose, "\ttotal SSB/Gy/GBp:",ccSSBsum/(sum_lc/1e+9)/dose)
0220 print("\nIf the geometry is considered as a whole nucleus:")
0221 #print("total DSB/Gy/GBp:",cDSBsum/(sum_lc/1e+9)/dose, "\ttotal SSB/Gy/GBp:",cSSBsum/(sum_lc/1e+9)/dose)
0222 print("total DSB/Gy/GBp:",cDSBsum/(sum_lc/1e+9)/dose, " --> DSB:", (cDSBsum-cDSBpsum-cDSBppsum)/(sum_lc/1e+9)/dose, "DSB+:", cDSBpsum/(sum_lc/1e+9)/dose, "DSB++:", cDSBppsum/(sum_lc/1e+9)/dose)
0223 print("total SSB/Gy/GBp:",cSSBsum/(sum_lc/1e+9)/dose, " --> SSB:", (cSSBsum-cSSBpsum-cSSBppsum)/(sum_lc/1e+9)/dose, "SSB+:", cSSBpsum/(sum_lc/1e+9)/dose, "SDSB++:", cSSBppsum/(sum_lc/1e+9)/dose)
0224 print("SSB/DSB ratio: ", cSSBsum/cDSBsum)
0225 print("Total dose (Gy): ", dose)
0226 #print(cDSBsum,ccDSBsum,cSSBsum,ccSSBsum)
0227
0228 print("\nThank you for using molecularDNA and supporting the Geant4-DNA collaboration.!")