!------------------------------------------------------------------------- ! Discover drug-adverse event interactions (part 4) ! ! In part 4 we are simply deconvolving for pairs 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_p4.f90 -O3 -o discover_ae_p4 ! ! Usage example, this will run the first 100 drugs and store the results to disk: ! echo "1,100" | ./discover_ae_p4 > 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_p4 > 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: ! echo "1,36" | ./discover_ae_p4_me > results/results_me_drug_drug_event_1+36.txt & ! echo "36,40" | ./discover_ae_p4_me > results/results_me_drug_drug_event_36+40.txt & ! echo "76,45" | ./discover_ae_p4_me > results/results_me_drug_drug_event_76+45.txt & ! echo "121,50" | ./discover_ae_p4_me > results/results_me_drug_drug_event_121+50.txt & ! ! Abada: ! echo "171,50" | ./discover_ae_p4_me > results/results_me_drug_drug_event_171+50.txt & ! echo "221,60" | ./discover_ae_p4_me > results/results_me_drug_drug_event_221+60.txt & ! echo "281,80" | ./discover_ae_p4_me > results/results_me_drug_drug_event_281+80.txt & ! echo "361,120" | ./discover_ae_p4_me > results/results_me_drug_drug_event_361+120.txt & ! echo "481,120" | ./discover_ae_p4_me > results/results_me_drug_drug_event_481+120.txt & ! ! Here is some ipython code to monitor the progress of these scripts while they run. ! def monitor_progress(total_estimate = 1.5e6): ! x = !ssh nick@orochi wc -l /Volumes/UsersHD/Users/nick/Stanford/AltmanLab/aers/fortran_data/results/*.txt ! y = !ssh nick@abada wc -l /Users/nick/Stanford/AltmanLab/aers/fortran_data/results/*.txt ! complete = int(x[-1].split()[-2]) + int(y[-1].split()[-2]) ! print >> sys.stdout, "Finished %d of %d (%.5f perc complete)" % (complete, total_estimate, complete/total_estimate) !------------------------------------------------------------------------- ! 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,i2 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 integer::max_reportcount, max_eventcount, num_reports_shared real(D)::temp_real integer,dimension(:),allocatable::stitch_reportcount, report_umlscount integer,dimension(:),allocatable::shared_reports 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. 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 real(D),dimension(:),allocatable::stitch_stitch_threshold integer::randvec(1100000),randptr,randint ! ================= INITIALIZATION ==================== ! Initialize the data matrices and vectors. read(5,*) start_index, 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= "medeffect_entity_counts.csv") read(10, *), num_stitch, num_reports, num_events 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_events, " adverse events." allocate(stitch_stitch(num_stitch, num_stitch)) allocate(stitch_stitch_threshold(num_stitch)) ! Load up how many reports we have for each drug. open(unit = 10, file= "medeffect_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= "medeffect_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= "medeffect_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= "medeffect_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 drug-drug correlation significance matrix. open(unit = 11, file= "medeffect_stitch_stitch_corr.csv") stitch_stitch = 0 do i=1,num_stitch read(11, *) stitch_stitch(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 if (stitch_reportcount(i).ge.50) then temp_real = 600.0 temp_vector = 0 do while ((sum(temp_vector).lt.50).and.(temp_real.gt.5)) temp_real = temp_real / 2.0 do j=1,num_stitch if (stitch_stitch(i,j).ge.temp_real) then do k=1,stitch_reportcount(j) temp_vector(stitch2report(j,k)) = 1 enddo endif enddo enddo if (sum(temp_vector).ge.50) then stitch_stitch_threshold(i) = temp_real endif endif enddo deallocate(temp_vector) ! write(0,*) stitch_stitch_threshold(1:5) ! =============== BOOTSTRAP THOSE STATISTICS ================== results = 0 do i=start_index,stop_index ! We only built bacgrounds for drugs which were on at least 50 reports, so ! we check for that now. if (stitch_reportcount(i).ge.50) then do i2=i+1,num_stitch ! We check that drug i2 also has >= 50 reports. if (stitch_reportcount(i2).ge.50) then ! Determine if we want to proceed with this pair (ie. share > 20 reports) allocate(temp_vector(num_reports)) temp_vector = 0 do k=1,stitch_reportcount(i) temp_vector( stitch2report(i,k) ) = temp_vector( stitch2report(i,k) ) + 1 enddo do k=1,stitch_reportcount(i2) temp_vector( stitch2report(i2,k) ) = temp_vector( stitch2report(i2,k) ) + 1 enddo num_reports_shared = 0 do k=1,num_reports if (temp_vector(k).eq.2) then num_reports_shared = num_reports_shared + 1 endif enddo allocate(shared_reports(num_reports_shared)) shared_reports = -1 l = 1 do k=1,num_reports if (temp_vector(k).eq.2) then shared_reports(l) = k l = l + 1 endif enddo deallocate(temp_vector) if (num_reports_shared.ge.20) then write(0,*) 'Working on drug pair: (', i, ', ', i2, '), fg reports: ', num_reports_shared 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 pair. !! Remember 'num_reports_shared' is holding the number of reports this pair shares. allocate( temp_vector( num_reports_shared ) ) temp_vector = 0 do k=1,num_reports_shared ! 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(shared_reports(k)) if (report2umls(shared_reports(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, num_reports_shared, 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 the "other" (ie. k) drug passes the muster for either drug in the we include it's reports. if ((stitch_stitch(i,k).ge.stitch_stitch_threshold(i)).or.(stitch_stitch(i2,k).ge.stitch_stitch_threshold(i2))) then do l=1,stitch_reportcount(k) temp_vector( stitch2report(k,l) ) = 1 enddo endif 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,I10,ES20.10,ES20.10,ES20.10,ES20.10,I10,I10)") i,i2,j,results,num_reports_shared,sum(temp_vector) temp_vector2 = 1 deallocate( temp_vector2 ) endif ! End check if this drug pair - adverse event combination has any reports. deallocate( temp_vector ) enddo endif ! End check that drug pair (i, i2) share at least 20 reports. deallocate(shared_reports) endif ! End check for drug i2 to have >= 50 reports. enddo endif ! End check for drug i to have >= 50 reports. 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