!------------------------------------------------------------------------- ! Discover drug-adverse event interactions (part 2) ! ! In part 2 we are simply deconvolving for single drugs using all reports. ! ! (c) 2010 Nicholas P. Tatonetti, Stanford University ! nick.tatonetti@stanford.edu !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! This script needs to be run from the directory which contains all of the data. ! For example, ! Orochi: /Volumes/UsersHD/Users/nick/Stanford/AltmanLab/aers/fortran_data/ ! Abada: /Users/nick/Stanford/AltmanLab/aers/fortran_data/ ! ! Move to that directory and compile this script on Orochi as follows: ! ifort ~/Dropbox/aers_src/fortran/discover_ae_p2.f90 -O3 -o discover_ae_p2 ! ! Usage example, this will run the first 100 drugs and store the results to disk: ! echo "1,100" | ./discover_ae_p2 > results_drug_event_1,100.txt ! ! Here is an example which starts 6 parallel processes. ! for (( i=1; i<=12; i+=2)); do echo $i,2 | ./discover_ae_p2 > results/results_drug_event_$i+2.txt & done ! ! This bash command is handy for testing out splits: ! for (( i=1; i<=787; i+=140)); do echo $i,140; done ! ! Here is an example where we run 4 threads on Orochi and 5 on Abada ! Orochi: for (( i=1; i<=490; i+=245)); do echo $i,245 | ./discover_ae_p2 > results/results_drug_event_$i+245.txt & done ! ! Abada: for (( i=491; i<=1567; i+=784)); do echo $i,783 | ./discover_ae_p2 > results/results_drug_event_$i+784.txt & done ! ! CLUSTER ! for ((i=1; i<962; i+= 100)); do bsub -e output/$i.txt -o results/$i.txt ./discover_ae_p2 $i 100; done !------------------------------------------------------------------------- ! PRECISION - select the basic precision for the computations. ! module precision integer,parameter::R=kind(1.0) ! IEEE single precision or similar integer,parameter::D=kind(1.0d0) ! IEEE double precision or similar integer,parameter::I64=selected_int_kind(18) ! 64-bit int on most machines end module precision module functions use precision contains subroutine vi_max(vector, length, max) integer,intent(in)::length integer,intent(in)::vector(length) integer,intent(inout)::max max = 0 do i=1,length if (i.eq.1) then max = vector(i) else if (vector(i).gt.max) then max = vector(i) endif endif enddo end subroutine vi_max subroutine next_rand(rand, randvec, randptr) integer,intent(inout)::randptr,rand integer,intent(in)::randvec(1100000) rand = randvec(randptr) randptr = randptr + 1 if (randptr.gt.1100000) then randptr = 1 endif end subroutine next_rand subroutine vector_mean(mu, vector, n) integer,intent(in)::n real(D),intent(out)::mu integer,intent(in)::vector(n) integer::i mu = 0.0 do i=1,n mu = mu + real(vector(i)) enddo mu = mu / real(n) end subroutine vector_mean subroutine vector_mean_real(mu, vector, n) integer,intent(in)::n real(D),intent(out)::mu real(D),intent(in)::vector(n) integer::i mu = 0.0 do i=1,n mu = mu + vector(i) enddo mu = mu / real(n) end subroutine vector_mean_real subroutine vector_std_real(std, mu, vector, n) integer,intent(in)::n real(D),intent(in)::mu real(D),intent(out)::std real(D),intent(in)::vector(n) integer::i std = 0.0 do i=1,n std = std + (vector(i) - mu)*(vector(i) - mu) enddo std = sqrt( std / real(n-1) ) end subroutine vector_std_real subroutine vector_sample(sample, size, vector, n, randvec, randptr) integer,intent(in)::size,n integer,intent(in)::vector(n) integer,intent(inout)::sample(size), randvec(1100000), randptr integer::i,randint do i=1,size call next_rand(randint, randvec, randptr) sample(i) = vector(mod(randint,n)+1) enddo end subroutine vector_sample subroutine bootstrap_statistics(mu, variance, vector, vector_length, randvec, randptr) real(D),intent(inout)::mu, variance integer,intent(in)::vector_length integer,intent(inout)::randptr integer,intent(inout)::randvec(1100000),vector(vector_length) real(D)::estimates(1000) integer::rand, sample_size, i integer,dimension(:),allocatable::sample sample_size = vector_length / 100 if (sample_size.lt.10) then sample_size = 10 endif if (sample_size.gt.10000) then sample_size = 10000 endif allocate(sample(sample_size)) sample = 0 do i=1,1000 call vector_sample( sample, sample_size, vector, vector_length, randvec, randptr ) call vector_mean( estimates(i), sample, sample_size) enddo call vector_mean_real( mu, estimates, 1000) call vector_std_real( variance, mu, estimates, 1000) end subroutine bootstrap_statistics end module functions program discover use precision use functions implicit none integer::i,j,k,l,m integer::start_index, step_size, stop_index ! This is initialized from the stdin to enable easy multi-processing. integer::num_stitch, num_reports, num_events, num_inds, min_bg_reports integer::max_reportcount, max_eventcount, max_indreportcount real(D)::temp_real integer,dimension(:),allocatable::stitch_reportcount, report_umlscount, ind_reportcount integer,dimension(:),allocatable::temp_vector, temp_vector2 integer,dimension(:,:),allocatable::report2umls ! Binary matrix that stores which events are on which reports integer,dimension(:,:),allocatable::stitch2report ! Binary matrix that stores which drugs are on which reports. integer,dimension(:,:),allocatable::ind2report ! Binary matrix that stores which indications are on which reports. real(D),dimension(4)::results ! A vector which stores fg_mean, fg_sd, fg_n, bg_mean, bg_sd, bg_n for the current interation. real(D),dimension(:,:),allocatable::stitch_stitch, stitch_ind real(D),dimension(:),allocatable::stitch_stitch_threshold, stitch_ind_threshold integer::randvec(1100000),randptr,randint character*10::argument ! ================= INITIALIZATION ==================== ! Initialize the data matrices and vectors. if (iargc().ne.2) then write(0,*) "Wrong number of arguments provided." endif call getarg(1, argument) read(argument, *) start_index call getarg(2, argument) read(argument, *) step_size stop_index = start_index + step_size - 1 write(0,*) "========= Discover Drug-AE: ", start_index, " to ~", stop_index, "==========" open(unit = 13, file = "random_numbers.vec") randptr = 726 do i=1,1100000 read(13,*) randvec(i) enddo open(unit = 10, file= "emr_entity_counts.csv") read(10, *), num_stitch, num_reports, num_events, num_inds if (num_stitch.lt.stop_index) then stop_index = num_stitch write(0,*) "Adjusting stop_index to ", stop_index endif write(0,*) "Found ", num_stitch, " drugs, ", num_reports, " reports, " write(0,*) "and ", num_inds, " indications, ", num_events, " adverse events." allocate(stitch_stitch(num_stitch, num_stitch)) allocate(stitch_stitch_threshold(num_stitch)) allocate(stitch_ind_threshold(num_stitch)) allocate(stitch_ind(num_stitch, num_inds)) ! Load up how many reports we have for each drug. open(unit = 10, file= "emr_stitch_reportcount.txt") allocate(stitch_reportcount(num_stitch)) max_reportcount = 0 do i=1,num_stitch read(10, *) stitch_reportcount(i) if (stitch_reportcount(i).gt.max_reportcount) then max_reportcount = stitch_reportcount(i) endif enddo ! Load up the stitch2report open(unit = 11, file= "emr_stitch_reportlist.csv") allocate( stitch2report( num_stitch, max_reportcount ) ) stitch2report = 0 do i=1,num_stitch allocate( temp_vector( stitch_reportcount(i) ) ) read(11, *) temp_vector do j=1,stitch_reportcount(i) stitch2report(i, j) = temp_vector(j) enddo deallocate( temp_vector ) enddo write(0,*) "Loaded stitch2report matrix (", num_stitch, ",", max_reportcount, ")" ! Load up how many events we have for each report. open(unit = 10, file= "emr_report_umlscount.txt") allocate(report_umlscount(num_reports)) max_eventcount = 0 do i=1,num_reports read(10, *) report_umlscount(i) if (report_umlscount(i).gt.max_eventcount) then max_eventcount = report_umlscount(i) endif enddo ! Load up the event data. open(unit = 11, file= "emr_report_umlslist.csv") allocate( report2umls( num_reports, max_eventcount ) ) report2umls = 0 do i=1,num_reports allocate( temp_vector( report_umlscount(i) ) ) read(11, *) temp_vector do j=1,report_umlscount(i) report2umls(i, j) = temp_vector(j) enddo deallocate( temp_vector ) enddo write(0,*) "Loaded report2umls matrix (", num_reports, ",", max_eventcount, ")" ! Load up the indicaiton report data. open(unit = 10, file= "emr_ind_reportcount.txt") allocate(ind_reportcount(num_inds)) max_indreportcount = 0 do i=1,num_inds read(10, *) ind_reportcount(i) if (ind_reportcount(i).gt.max_indreportcount) then max_indreportcount = ind_reportcount(i) endif enddo open(unit = 11, file= "emr_ind_reportlist.csv") allocate( ind2report( num_inds, max_indreportcount ) ) ind2report = 0 do i=1,num_inds allocate( temp_vector( ind_reportcount(i) ) ) read(11, *) temp_vector do j=1,ind_reportcount(i) ind2report(i, j) = temp_vector(j) enddo deallocate( temp_vector ) enddo write(0,*) "Loaded ind2report matrix (", num_inds, ",", max_indreportcount, ")" ! Load up the drug-drug correlation significance matrix. open(unit = 11, file= "emr_stitch_stitch_est.csv") ! open(unit = 11, file= "aers_stitch_stitch_corr.csv") stitch_stitch = 0 do i=1,num_stitch read(11, *) stitch_stitch(i,:) enddo open(unit = 11, file= "emr_stitch_ind_est.csv") ! open(unit = 11, file= "aers_stitch_ind_corr.csv") stitch_ind = 0 do i=1,num_stitch read(11, *) stitch_ind(i,:) enddo ! =============== IDENTIFY THE DRUG-DRUG BG THRESHOLD ================== write(0,*) "Identifying the drug-drug correlation thresholds..." stitch_stitch_threshold = -1 allocate(temp_vector(num_reports)) do i=1,num_stitch write(0,*) "Working on drug ", i if (stitch_reportcount(i).ge.100) then min_bg_reports = stitch_reportcount(i) ! "e" ! write(0,*) "Finding drug-drug background with minimum of ", min_bg_reports, " reports." ! temp_real = 1200.0 temp_real = 1.0 temp_vector = 0 ! If the drug is on every report then it clogs up the system if (min_bg_reports.gt.(num_reports/2)) then stitch_stitch_threshold(i) = 0.0 else do while ((sum(temp_vector).lt.min_bg_reports).and.(temp_real.lt.10000)) ! temp_real = temp_real / 2.0 temp_real = temp_real + 1.0 do j=1,num_stitch if (stitch_stitch(i,j).ge.(1.0/temp_real)) then ! For each report that this other drug is on... do k=1,stitch_reportcount(j) ! Set the temp_vector position for the report to 1. temp_vector(stitch2report(j,k)) = 1 enddo endif enddo ! write(0,*) "Found ", sum(temp_vector), " before removing reports for this drugs." ! Set any of the report positions of this drug back to 0. do k=1,stitch_reportcount(i) temp_vector(stitch2report(i,k)) = 0 enddo ! write(0,*) "Found ", sum(temp_vector), " after removing reports for this drugs." enddo stitch_stitch_threshold(i) = (1.0/temp_real) endif endif enddo deallocate(temp_vector) ! =============== IDENTIFY THE DRUG-INDICATION BG THRESHOLD ================== write(0,*) "Identifying the drug-ind correlation thresholds..." stitch_ind_threshold = -1 allocate(temp_vector(num_reports)) do i=1,num_stitch if (stitch_reportcount(i).ge.100) then write(0,*) "Working on drug ", i min_bg_reports = stitch_reportcount(i) ! write(0,*) "Finding drug-indication background with minimum of ", min_bg_reports, " reports." ! temp_real = 1200 ! use this for t_statistic cutoffs temp_real = 1.0 temp_vector = 0 ! If the drug is on every report then it clogs up the system if (min_bg_reports.gt.(num_reports/3)) then stitch_ind_threshold(i) = 0.0 else do while ((sum(temp_vector).lt.min_bg_reports).and.(temp_real.lt.1000)) ! temp_real = temp_real / 2.0 ! use for t_statistic cutoffs temp_real = temp_real + 1.0 ! use for estimate cutoffs do j=1,num_inds if (stitch_ind(i,j).ge.(1.0/temp_real)) then ! For each report that this indication is on... do k=1,ind_reportcount(j) ! set the temp_vector position for the report to 1. temp_vector(ind2report(j,k)) = 1 enddo endif enddo !write(0,*) "Found ", sum(temp_vector), " before removing reports for this drugs." ! Set any of the report positions of this drug back to 0. do k=1,stitch_reportcount(i) temp_vector(stitch2report(i,k)) = 0 enddo !write(0,*) "Found ", sum(temp_vector), " after removing reports for this drugs." enddo stitch_ind_threshold(i) = (1.0/temp_real) endif endif enddo deallocate(temp_vector) ! =============== BOOTSTRAP THOSE STATISTICS ================== results = 0 do i=start_index,stop_index write(0,*) 'Working on drug ', i-start_index+1, " of ", step_size if (stitch_reportcount(i).ge.100) then do j=1,num_events results = 0 ! Given a drug, i, and an adverse event, j... ! 1. Construct a vector of reports for this drug. allocate( temp_vector( stitch_reportcount(i) ) ) temp_vector = 0 do k=1,stitch_reportcount(i) ! 2. Set vector element equal to 1 if the event is on this report ! "k" is now the index of the report id we want to check for this event do l=1,report_umlscount(stitch2report(i,k)) if (report2umls(stitch2report(i,k),l).eq.j) then temp_vector(k) = 1 endif enddo enddo if (sum(temp_vector).gt.0) then ! 3. We only need to do the computation on drugs which have this adverse event. call bootstrap_statistics(results(1), results(2), temp_vector, stitch_reportcount(i), randvec, randptr) deallocate( temp_vector ) allocate( temp_vector( num_reports ) ) temp_vector = 0 ! Turn "on" the reports for the correlated drugs. do k=1,num_stitch if (stitch_stitch(i,k).ge.stitch_stitch_threshold(i)) then do l=1,stitch_reportcount(k) temp_vector( stitch2report(k,l) ) = 1 enddo endif enddo ! Turn "on" the reports for the correlated indications. do k=1,num_inds if (stitch_ind(i,k).ge.stitch_ind_threshold(i)) then do l=1,ind_reportcount(k) temp_vector( ind2report(k,l) ) = 1 enddo endif enddo ! Turn "off" the reports for the current drug ("e" and "f" versions of algorithm) do k=1,stitch_reportcount(i) temp_vector(stitch2report(i,k)) = 0 enddo allocate( temp_vector2( sum(temp_vector) ) ) temp_vector2 = 0 m = 0 do k=1,num_reports if (temp_vector(k).eq.1) then m = m + 1 do l=1,report_umlscount(k) if (report2umls(k,l).eq.j) then temp_vector2(m) = 1 endif enddo endif enddo call bootstrap_statistics( results(3), results(4), temp_vector2, sum(temp_vector), randvec, randptr) write(6,"(I10,I10,ES20.10,ES20.10,ES20.10,ES20.10,I10,I10,ES20.10,ES20.10)") i,j,results,stitch_reportcount(i),sum(temp_vector),stitch_stitch_threshold(i),stitch_ind_threshold(i) temp_vector2 = 1 deallocate( temp_vector2 ) endif deallocate( temp_vector ) enddo endif enddo ! ========= DEBUGGING STUFF ========== ! do i=1,num_stitch ! k = 0 ! do j=1,max_reportcount ! if (stitch2report(i,j).ne.0) then ! k = k + 1 ! endif ! enddo ! if (k.ne.stitch_reportcount(i)) then ! write(0,*) "Error in reading stitch2report data for stitch at index ", i ! endif ! enddo ! ! do i=1,num_reports ! k = 0 ! do j=1,max_eventcount ! if (report2umls(i,j).ne.0) then ! k = k + 1 ! endif ! enddo ! if (k.ne.report_umlscount(i)) then ! write(0,*) "Error in reading stitch2report data for stitch at index ", i ! endif ! enddo ! ! write(0,*) stitch_reportcount(1:10) ! We're done. deallocate(stitch_reportcount) end program discover