Merge pull request #90 from justinwang97/master

update indirect marker effect MCMC sampling for GWAS

src/1.JWAS/src/JWAS.jl

@@ -338,7 +338,15 @@ function runMCMC(mme::MME,df;
  versioninfo()
  printstyled("\n\nThe analysis has finished. Results are saved in the returned ",bold=true)
  printstyled("variable and text files. MCMC samples are saved in text files.\n\n\n",bold=true)
+
+ # make MCMC samples for indirect marker effect
+ if causal_structure != false
+ JWAS.generate_indirect_marker_effect_sample(mme.lhsVec,output_folder,causal_structure,"structure_coefficient_MCMC_samples.txt")
+
+ end
+
  return mme.output
+
 end
 ################################################################################
 # Print out Model or MCMC information
@@ -397,6 +405,9 @@ function describe(model::MME;data=false)
  if model.MCMCinfo != false && model.MCMCinfo.printout_model_info == true
  getMCMCinfo(model)
  end
+
+
+
 end
 
 """
@@ -506,4 +517,7 @@ function getMCMCinfo(mme)
  @printf("%-30s %20.3f\n","marker effect variances:",mme.df.marker)
  end
  @printf("\n\n\n")
+
+
+
 end

src/1.JWAS/src/MCMC/MCMC_BayesianAlphabet.jl

@@ -316,7 +316,7 @@ function MCMC_BayesianAlphabet(mme,df)
  # 4. Causal Relationships among Phenotypes (Structure Equation Model)
  ########################################################################
  if is_multi_trait && causal_structure != false
- sample4λ = get_Λ(Y,mme.R,ycorr,Λy,mme.ySparse,causal_structure) #no missing phenotypes
+ sample4λ,sample4λ_vec = get_Λ(Y,mme.R,ycorr,Λy,mme.ySparse,causal_structure) #no missing phenotypes
  end
  ########################################################################
  # 5. Latent Traits
@@ -353,7 +353,7 @@ function MCMC_BayesianAlphabet(mme,df)
  #mean and variance of posterior distribution
  output_MCMC_samples(mme,mme.R,(mme.pedTrmVec!=0 ? inv(mme.Gi) : false),outfile)
  if causal_structure != false
- writedlm(causal_structure_outfile,sample4λ',',')
+ writedlm(causal_structure_outfile,sample4λ_vec',',')
  end
  end
  ########################################################################
@@ -388,6 +388,14 @@ function MCMC_BayesianAlphabet(mme,df)
  output_folder*"/MCMC_samples_genetic_variance(REML)"*"_"*Mi.name*".txt")
  end
  end
+
+
+ ############################################################################
+ # Compute the indirect marker effects
+ ############################################################################
+ output_folder
+
+
  output=output_result(mme,output_folder,
  mme.solMean,mme.meanVare,
  mme.pedTrmVec!=0 ? mme.G0Mean : false,

src/1.JWAS/src/structure_equation_model/SEM.jl

@@ -57,6 +57,9 @@ function SEM_setup(wArray,causal_structure,mme)
  Λy = kron(Λ,sparse(1.0I,nobs,nobs))*mme.ySparse
  causal_structure_filename = "structure_coefficient_MCMC_samples.txt"
  causal_structure_outfile = open(causal_structure_filename,"w") #write MCMC samples for Λ to a txt file
+
+
+
  return Y,Λy,causal_structure_outfile
 end
 # Get Y for all individuals ordered as individuals within traits (fully simultaneous model)
@@ -138,12 +141,14 @@ function get_Λ(Y,R,Λycorr,Λy,y,causal_structure)
  # Λ = [1 -λ12
  # -λ21 1]
  λ = rand(MvNormal(mu,var))
- Λ = I - tranform_lambda(λ,causal_structure)
 
+ causal_matrix = tranform_lambda(λ,causal_structure)
+ Λ = I - causal_matrix
+ λ_vec = vec(causal_matrix)
 
  Λy[:] = kron(Λ,sparse(1.0I,nind,nind))*y
  Λycorr[:] = ycorr - y + Λy #add new Λy
- return λ
+ return λ,λ_vec
 end
 
 function tranform_lambda(lambda,causal_structure)
@@ -156,3 +161,80 @@ function tranform_lambda(lambda,causal_structure)
  end
  return Lambda
 end
+
+
+# generate the MCMC samples for indirect marker effect
+# In each iteration, we want to form the following equation
+#
+# [0 0 0 [marker1_y1 marker2_y1 ... [0 0
+# λ12 0 0 * marker1_y2 marker2_y2 ... = λ12*marker1_y1 λ12*marker2_y1
+# λ13 0 0] marker1_y3 marker2_y3 ...] λ13*marker1_y1 λ12*marker2_y1 ]
+# = Λ * [α1,α2,...] = indirect effect
+
+function generate_indirect_marker_effect_sample(pheno_vec,output_folder,causal_structure,structure_coefficient_path)
+
+
+ number_traits = size(causal_structure,1) # the row number of causal structure matrix is number of traits
+ trait_vec = string.(pheno_vec) # transform the symbol to string
+ λ_file = CSV.read(structure_coefficient_path, DataFrame,header = false)
+
+ direct_effect_sample = Dict()
+ io_diction = Dict()
+ # read the direct effect sample file and create indirect sample file for each trait
+ for i in 1:number_traits
+ direct_file_name = output_folder*"/MCMC_samples_marker_effects_genotypes_"*trait_vec[i]*".txt"
+ direct_effect_sample[trait_vec[i]] = CSV.read(direct_file_name,DataFrame,header = true)
+
+ indirect_file_name = output_folder*"/MCMC_samples_indirect_marker_effects_genotypes_"*trait_vec[i]*".txt"
+ io_diction[trait_vec[i]] = indirect_file_name
+ end
+
+ number_sample = size(direct_effect_sample[trait_vec[1]],1)
+ number_marker = size(direct_effect_sample[trait_vec[1]],2)
+
+ marker_header = permutedims(names(direct_effect_sample[trait_vec[1]]))
+
+ # write marker header for each indirect effect file
+ for i in 1:number_traits
+ open(io_diction[trait_vec[i]], "a") do io
+ writedlm(io, marker_header ,", ")
+ end;
+
+ end
+
+ # compute the indirect effect in each sample and write to the target file
+ for i in 1:number_sample
+ direct_effect_current_sample = zeros(number_traits,number_marker)
+ λ_sample = Vector(λ_file[i,1:end])
+ for j in 1:number_traits
+ current_effect = direct_effect_sample[trait_vec[j]][i,1:end]
+ current_effect = Vector(current_effect)
+
+ direct_effect_current_sample[j,1:end] = current_effect'
+
+ end
+
+ causal_matrix = reshape(λ_sample ,number_traits,number_traits)
+ indirect_effect_current_sample = compute_indirect_effect(causal_matrix,direct_effect_current_sample )
+
+ for i in 1:number_traits
+ indirect_effect_this_trait = indirect_effect_current_sample[i,1:end]
+ open(io_diction[trait_vec[i]], "a") do io
+ writedlm(io, indirect_effect_this_trait' ,", ")
+ end;
+
+ end
+ end
+
+
+end
+
+# compute the indirect marker effect based on the formula in the SEM paper
+function compute_indirect_effect(Λ,marker_effects )
+ number_traits,number_marker = size(marker_effects)
+ result = zeros(number_traits,number_marker)
+ for i in 1:number_traits-1
+ result += Λ^i * marker_effects
+ end
+ return (result)
+end