dna/dnadamage1/scandamages.C

0001 // -------------------------------------------------------------------
0002 // $Id: scandamages.C
0003 // -------------------------------------------------------------------
0004 // 13/10/2023 - Le Tuan Anh renamed and improved
0005 // *********************************************************************
0006 // To execute this macro under ROOT after your simulation ended,
0007 //   1 - launch ROOT v6.xx (usually type 'root' at your machine's prompt)
0008 //   2 - type '.X analysis.C' at the ROOT session prompt
0009 // *********************************************************************
0010
0011 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0012
0013 struct Damage;
0014 struct Cluster;
0015 void CreateSDDfile(std::map<int,std::vector<Damage>>,long int);
0016 void ExtractDirectAndIndirectDamages(std::map<int,std::vector<Damage>>&,long int&,double,double);
0017 void ClassifyDamages(std::map<int,std::vector<Damage>>&,int,
0018                     std::map<int,vector<Cluster>>&,int verbose =0);
0019
0020 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0021 // play with scandamages() by changing firect/indirect criterion selection or Clusterdistance
0022
0023 void scandamages()
0024 {
0025     ifstream fs("output_t0.root");
0026     if ( fs.good() ) {
0027       system ("rm -rf output.root");
0028       system ("hadd output.root output_*.root");
0029     }
0030
0031     double fEnergyThresholdForDirectSelection = 17.5; // eV
0032     double fProbabilityForIndirectionSelection = 0.40; // 40 %
0033     int fClusterdistance = 10; // bp, minium distance that separate two clusters
0034
0035     std::map<int,std::vector<Damage>> fDamagesEventMap;
0036     std::map<int,vector<Cluster>> fClustersMap;
0037     long int totalDamages = 0;
0038
0039     ExtractDirectAndIndirectDamages(fDamagesEventMap,totalDamages,
0040                                     fEnergyThresholdForDirectSelection,
0041                                     fProbabilityForIndirectionSelection);
0042
0043     ClassifyDamages(fDamagesEventMap,fClusterdistance,fClustersMap);
0044     // print results in each event for checking:
0045     //ClassifyDamages(fDamagesEventMap,fClusterdistance,fClustersMap,1);
0046
0047     CreateSDDfile(fDamagesEventMap,totalDamages);
0048 }
0049
0050 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0051
0052 struct Damage
0053 {
0054     Damage(int a, int b, int c=-1): strand(a),copyNb(b), cause(c){}
0055     int strand;
0056     int copyNb;
0057     int cause; // Unknown = -1, Direct = 0, Indirect = 1
0058     bool operator<(const Damage& r) const {return copyNb<r.copyNb;} // for sorting
0059 };
0060
0061 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0062
0063 struct Cluster // to score SSB or DSB
0064 {
0065     std::vector<Damage> fDamagesVector;
0066     bool HasLeftStrand=false, HasRightStrand=false;
0067     int firstCopyNb;
0068 };
0069
0070 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0071
0072 void
0073 ClassifyDamages(std::map<int,std::vector<Damage>> &fDamages,int fClusterditance,
0074                 std::map<int,vector<Cluster>>& fClustersMap,int verbose)
0075 {
0076     long int dirSB = 0, indirSB = 0;
0077     long int nDSB = 0, nsDSB = 0, ncDSB=0, nSSB = 0;
0078     std::cout<<"-----------------------------------------------"<<std::endl;
0079     if (verbose>0)std::cout<<" ------> Print few events for checking <------"<<std::endl;
0080     if (verbose>0)std::cout<<"\t\tCopyNb\tStrand\tCause"<<std::endl;
0081     for (auto &[eventID, dmgVector] : fDamages) {
0082         std::sort(dmgVector.begin(),dmgVector.end());
0083         if (verbose>0) std::cout<<"Event #"<<eventID<<":  "<<std::endl;
0084         unsigned long fSBcountThisEvent = 0;
0085         long copyNBofPrevDamage = 0;
0086         long firstCopyNb;
0087         std::map<int,Cluster> clustersInthisEvent;
0088         for (auto dmg : dmgVector) {
0089             if (dmg.cause == 0) dirSB++;
0090             if (dmg.cause == 1) indirSB++;
0091             if (verbose>0) std::cout<<"\t\t"<<dmg.copyNb<<"\t"
0092                                     <<dmg.strand<<"\t"<<dmg.cause<<"\t"<<std::endl;
0093             if (clustersInthisEvent.size()==0 ||
0094                 ((dmg.copyNb - copyNBofPrevDamage) >=  fClusterditance)) {
0095                 Cluster newCluster;
0096                 newCluster.fDamagesVector.push_back(dmg);
0097                 copyNBofPrevDamage = dmg.copyNb;
0098                 firstCopyNb = dmg.copyNb;
0099                 if (dmg.strand == 0) newCluster.HasLeftStrand = true;
0100                 if (dmg.strand == 1) newCluster.HasRightStrand = true;
0101                 clustersInthisEvent.insert({firstCopyNb,newCluster});
0102             } else {
0103                 clustersInthisEvent[firstCopyNb].fDamagesVector.push_back(dmg);
0104                 copyNBofPrevDamage = dmg.copyNb;
0105                 if (dmg.strand == 0 ) clustersInthisEvent[firstCopyNb].HasLeftStrand = true;
0106                 if (dmg.strand == 1 ) clustersInthisEvent[firstCopyNb].HasRightStrand = true;
0107             }
0108         }
0109
0110         for (auto [firstCpNb,acluster] : clustersInthisEvent) {
0111             if (fClustersMap.find(eventID) == fClustersMap.end()) {
0112                 std::vector<Cluster> sclv{acluster};
0113                 fClustersMap.insert({eventID,sclv});
0114             } else {
0115                 fClustersMap[eventID].push_back(acluster);
0116             }
0117         }
0118         clustersInthisEvent.clear();
0119     }
0120     if (verbose>0)std::cout<<" \n------> Checking SSB and DSB:"<<std::endl;
0121     if (verbose>0)std::cout<<"\t\t#SSB\t#DSB\t#cDSB\t#sDSB"<<std::endl;
0122
0123     for (auto [evt,clustervector] : fClustersMap) {
0124         long int nDSBinThisEvt = 0, nSSBinThisEvt = 0;
0125         long int nsDSBinThisEvt = 0, ncDSBinThisEvt = 0;
0126         if (verbose>0) std::cout<<"Event #"<<evt<<":  "<<std::endl; // maximum 5 event for cheking
0127         for (auto sb : clustervector) {
0128             if (sb.HasLeftStrand && sb.HasRightStrand) {
0129                 nDSBinThisEvt++;
0130                 if (sb.fDamagesVector.size()>2) ncDSBinThisEvt++;
0131                 else nsDSBinThisEvt++;
0132             }
0133             else nSSBinThisEvt++;
0134         }
0135         nSSB += nSSBinThisEvt;
0136         nDSB += nDSBinThisEvt;
0137         ncDSB += ncDSBinThisEvt;
0138         nsDSB += nsDSBinThisEvt;
0139         if (verbose>0) std::cout<<"\t\t"<<nSSBinThisEvt<<"\t"<<nDSBinThisEvt
0140                                 <<"\t"<<ncDSBinThisEvt<<"\t"<<nsDSBinThisEvt<<std::endl;
0141     }
0142     std::cout<<"-----------------------------------------------"<<std::endl;
0143     std::cout<<"\tSummary of results:"
0144                 <<"\n #Total SBs = "<<(dirSB+indirSB)
0145                 <<"\n #dirSB     = "<<dirSB <<" \t;\t #indirSB = "<<indirSB
0146                 <<"\n #SSB       = "<<nSSB  <<" \t;\t #DSB     = "<<nDSB
0147                 <<"\n #sDSB      = "<<nsDSB <<" \t;\t #cDSB    = "<<ncDSB<<std::endl;
0148 }
0149
0150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0151
0152 void ExtractDirectAndIndirectDamages(std::map<int,std::vector<Damage>> &fDamagesEventMap,
0153 long int &totalDamages,double Ethresh,double ChemPro)
0154 {
0155     TRandom rdomN;
0156
0157     vector<Molecule> fMolecules = molecule();
0158     using Table = map<int,array<map<int,double>, 2>>;
0159     double xphy, yphy, zphy;
0160     double xChe, yChe, zChe;
0161     double xrad, yrad, zrad;
0162
0163     map<int,int> CopyNumTable;
0164     Table physTable;
0165
0166     double EnergyDeposit;
0167     double kineticEnergyDifference;
0168     int flagVolume;
0169     double copyN;
0170     int EventIDp;
0171     int EventIDc;
0172
0173     TFile f("output.root");
0174
0175 // Load the file, the directory and the ntuple
0176     TDirectoryFile* d = dynamic_cast<TDirectoryFile*>(f.Get("ntuple") );
0177 //analyse Physical stage
0178     TTree* tPhys = dynamic_cast<TTree*> (d->Get("ntuple_1") );
0179     tPhys->SetBranchAddress("x",&xphy);
0180     tPhys->SetBranchAddress("y",&yphy);
0181     tPhys->SetBranchAddress("z",&zphy);
0182     tPhys->SetBranchAddress("edep",&EnergyDeposit);
0183     tPhys->SetBranchAddress("diffKin",&kineticEnergyDifference);
0184     tPhys->SetBranchAddress("volumeName",&flagVolume);
0185     tPhys->SetBranchAddress("CopyNumber",&copyN);
0186     tPhys->SetBranchAddress("EventID",&EventIDp);
0187
0188     auto entryPhyNumber = tPhys->GetEntries();
0189
0190     bool addCopyN = true;
0191     for(decltype(entryPhyNumber) i = 0; i < entryPhyNumber; ++i)
0192     {
0193         tPhys->GetEntry(i);
0194 //avoid histone
0195         if(flagVolume == DNAVolumeType::histone)
0196         {
0197             continue;
0198         }
0199 //Determine the strand
0200         int strand = -1;
0201         bool isLeftStrand = DNAVolumeType::deoxyribose1 == flagVolume ||
0202                             DNAVolumeType::phosphate1 == flagVolume ||
0203                             DNAVolumeType::deoxyribose1_water == flagVolume ||
0204                             DNAVolumeType::phosphate1_water == flagVolume;
0205
0206         bool isRightStrand = DNAVolumeType::deoxyribose2 == flagVolume ||
0207                              DNAVolumeType::phosphate2 == flagVolume ||
0208                              DNAVolumeType::deoxyribose2_water== flagVolume ||
0209                              DNAVolumeType::phosphate2_water == flagVolume;
0210
0211
0212         if( isLeftStrand )
0213         {
0214             strand = 0;
0215         }
0216         else if( isRightStrand )
0217         {
0218             strand = 1;
0219         }
0220         else
0221         {
0222             //in case molecules are not assigned "strand"
0223             continue;
0224         }
0225
0226 //Determine the name
0227         bool isDeoxyribode = flagVolume == DNAVolumeType::deoxyribose1 ||
0228                              flagVolume == DNAVolumeType::deoxyribose2 ||
0229                              flagVolume == DNAVolumeType::deoxyribose1_water ||
0230                              flagVolume == DNAVolumeType::deoxyribose2_water;
0231
0232         bool  isPhosphate =  flagVolume == DNAVolumeType::phosphate1 ||
0233                              flagVolume == DNAVolumeType::phosphate2 ||
0234                              flagVolume == DNAVolumeType::phosphate1_water||
0235                              flagVolume == DNAVolumeType::phosphate2_water;
0236
0237         if(isDeoxyribode || isPhosphate)
0238         {
0239             physTable[EventIDp][strand][copyN] += EnergyDeposit;
0240         }
0241
0242         if(physTable[EventIDp][strand][copyN] < Ethresh)
0243         {
0244             continue;
0245         }
0246
0247         if(CopyNumTable.empty())
0248         {
0249             CopyNumTable.insert(pair<int,int>(copyN,strand));
0250             if (fDamagesEventMap.find(EventIDp) == fDamagesEventMap.end()) {
0251                 std::vector<Damage> admg{Damage(strand,copyN,0)};
0252                 fDamagesEventMap.insert({EventIDp,admg});
0253             } else {
0254                 fDamagesEventMap[EventIDp].push_back(Damage(strand,copyN,0));
0255             }
0256             totalDamages++;
0257         }
0258         else
0259         {
0260             addCopyN = true;
0261         }
0262
0263         auto itCopyNumTable = CopyNumTable.find(copyN);
0264         if (itCopyNumTable != CopyNumTable.end())
0265         {
0266             if (itCopyNumTable->second == strand)
0267             {
0268                 addCopyN = false;
0269             }
0270         }
0271
0272         if(addCopyN)
0273         {
0274             CopyNumTable.insert(pair<int,int>(copyN,strand));
0275             if (fDamagesEventMap.find(EventIDp) == fDamagesEventMap.end()) {
0276                 std::vector<Damage> admg{Damage(strand,copyN,0)};
0277                 fDamagesEventMap.insert({EventIDp,admg});
0278             } else {
0279                 fDamagesEventMap[EventIDp].push_back(Damage(strand,copyN,0));
0280             }
0281             totalDamages++;
0282         }
0283     }
0284 //Chemistry analyse
0285     TTree* tChem = dynamic_cast<TTree*> (d->Get("ntuple_2") );
0286     tChem->SetBranchAddress("x",&xrad);
0287     tChem->SetBranchAddress("y",&yrad);
0288     tChem->SetBranchAddress("z",&zrad);
0289     tChem->SetBranchAddress("EventID",&EventIDc);
0290
0291     auto entryNChem = tChem->GetEntries() ;
0292
0293     for(decltype(entryNChem) i=0; i < entryNChem; ++i)
0294     {
0295         tChem->GetEntry(i);
0296         ThreeVector<double> DNAElement(xrad,yrad,zrad);
0297
0298         for(const auto& moleculeIt : fMolecules)
0299         {
0300             if(moleculeIt.fPosition == DNAElement)
0301             {
0302                 int strand = -1;
0303                 if(moleculeIt.fName == "deoxyribose1")
0304                 {
0305                     strand = 0;
0306                 }
0307                 else if(moleculeIt.fName == "deoxyribose2")
0308                 {
0309                     strand = 1;
0310                 }
0311                 else
0312                 {
0313                     string msg = "Unknown DNA component";
0314                     throw std::runtime_error(msg);
0315                 }
0316
0317                 if(rdomN.Rndm() > ChemPro)//probability of reaction make damages
0318                 {
0319                     continue;
0320                 }
0321
0322                 if (!CopyNumTable.empty())
0323                 {
0324                     addCopyN = true;
0325                 }
0326
0327                 auto itCopyNumTable = CopyNumTable.find(moleculeIt.fCopyNumber);
0328
0329                 if (itCopyNumTable != CopyNumTable.end())
0330                 {
0331                     if (itCopyNumTable->second == strand)
0332                     {
0333                         addCopyN = false;
0334                     }
0335                 }
0336                 if(addCopyN)
0337                 {
0338                     if (fDamagesEventMap.find(EventIDc) == fDamagesEventMap.end()) {
0339                         std::vector<Damage> admg{Damage(strand,moleculeIt.fCopyNumber,1)};
0340                         fDamagesEventMap.insert({EventIDc,admg});
0341                     } else {
0342                         fDamagesEventMap[EventIDc].push_back(
0343                                                 Damage(strand,moleculeIt.fCopyNumber,1));
0344                     }
0345                     totalDamages++;
0346                 }
0347             }
0348         }
0349     }
0350 }
0351
0352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0353
0354 void CreateSDDfile(std::map<int,std::vector<Damage>> fDamagesEventMap,long int totalDamages)
0355 {
0356     ofstream ofile("SDDFormat.txt");
0357     if(ofile.is_open())
0358     {
0359         ofile<<'\n'
0360             <<"SDD version, SDDv1.0;\n"
0361             <<"Software, chromatin fiber Model;\n"
0362             <<"Author contact, Carmen Villagrasa,carmen.villagrasa@irsn.fr, "
0363               "02/04/2018, Melyn et al.,Sc. Rep. 7 (2017)11923;\n"
0364             <<"Simulation Details, DNA damages from direct and "
0365               "indirect effects;\n"
0366             <<"Source, Monoenergetic cylinder-parallel proton beam uniformly "
0367               "in 5 nm radius from left side of the cube"
0368             <<" exposing chromatin fiber. Energy: 5 keV;\n"
0369             <<"Source type, 1;\n"
0370             <<"Incident particles, 2212;\n"
0371             <<"Mean particle energy, 5 keV;\n"
0372             <<"Energy distribution, M, 0;\n"
0373             <<"Particle fraction, 1.0;\n"
0374             <<"Dose or fluence, 0, 0;\n"
0375             <<"Dose rate, 0.0;\n"
0376             <<"Irradiation target, Simple spherical chromatin fiber model"
0377               " in a voxel 40 nm;\n"
0378             <<"Volumes, 0,20,20,20,-20,-20,-20,2,15,15,20,-15,-15,20;\n"
0379             <<"Chromosome sizes, ;\n"
0380             <<"DNA Density, ;"
0381             <<"Cell Cycle Phase, G0/G1;\n"
0382             <<"DNA Structure, 4, 1;\n"
0383             <<"In vitro / in vivo, 0;\n"
0384             <<"Proliferation status, 1;\n"
0385             <<"Microenvironment, 27, 0.0;\n"
0386             <<"Damage definition, 1, 1, 10, -1, 17.5;\n"
0387             <<"Time, 2.5ns;\n"
0388             <<"Damage and primary count, "<<totalDamages<<", 0;\n"
0389             <<"Data entries, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0;\n"
0390             <<"Data field explaination, Field 1: [2]-eventID, Field 3: [4]-strand, "
0391             <<"Field 4:damage position (copynb), Field 5: Cause (direct: [1]=0) "
0392             <<" (indirect: [1]=1), Field 6: Damage types (Base:[1]>0) (Backbone: [2]>0);\n"
0393             <<"\n"
0394             <<"***EndOfHeader***;\n"
0395             <<endl;
0396 //For detail, please refer to Schuemann, J., et al. (2019). A New Standard DNA
0397 //Damage (SDD) Data Format. Radiation Research, 191(1), 76.
0398
0399     bool Primaryflag = true;
0400     for (auto const [eventID,dmgvector] : fDamagesEventMap)
0401     {
0402         Primaryflag = true; // update Primaryflag for new primary
0403         for (auto const dmg:dmgvector) {
0404             //Field 1 Calassification
0405             int newEvtFlag = 0; // = 2 if new event;
0406             if (Primaryflag) {
0407                 newEvtFlag = 2;
0408                 Primaryflag = false;
0409             }
0410             ofile<<newEvtFlag<<", "<<eventID<<"; ";
0411             //Field 3 Chromosome IDs and strand
0412             ofile<<0<<", "<<0<<", "<<0<<", "<<dmg.strand<<"; ";
0413             //Field 4, Chromosome position
0414             ofile<<dmg.copyNb<<"; ";
0415             //Field 5, Cause: Unknown = -1, Direct = 0, Indirect = 1
0416             ofile<<dmg.cause<<", "<<0<<", "<<0<<"; ";
0417             //Field 6, Damage types: only Backbone is considered here
0418             ofile<<0<<", "<<1<<", "<<0<<"; ";
0419             ofile<<"\n";
0420         }
0421     }
0422
0423     ofile<<endl;
0424     }
0425
0426     ofile.close();
0427 }
0428
0429 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....