https://github.com/vaskarageorg/SCA_MR
Raw File
Tip revision: 33a6b14e3083b370e7f0a209b18b45a26c5529f1 authored by Vasilios Karageorgiou on 25 May 2023, 18:42:06 UTC
Update README.md
Tip revision: 33a6b14
correl_exps.R
correl_expos=function(nX,nG,nbl1,n_ind,sum_stats=T){
  maf= runif(n = nG,min = 0.01,max = 0.5)
  
  G = matrix(nrow=n_ind,ncol = nG)
  for (nh in 1:ncol(G)){ G[,nh] = matrix(rbinom(n_ind, size = 2, prob = maf[nh]),ncol = 1)}
  
  U = matrix(rnorm(n = n_ind, 0, 5))
  
  epsilon_X=matrix(rnorm(n = n_ind*nX,0,4),ncol = nX)
  epsilon_y=matrix(rnorm(n = n_ind*1,0,3),ncol = 1)
  
  n_blocks=nbl1;bre1 = seq(1, nG, nG %/% n_blocks);bre1 = append(x = bre1,values = nG)
  
  list_causal=list()
  for (jl in 2:length(bre1)){ list_causal[[jl-1]] = seq(bre1[jl-1],bre1[jl]) }
  onelist1=which(unlist(lapply(list_causal,length))==1)
  if (length(onelist1)!=0){list_causal = list_causal[-onelist1]}
  
  gamma_SNP_X=matrix(nrow = nG,ncol = nX)
  gamma_SNP_X_mu=seq(4,5,length.out = length(list_causal))
  
  X_per_block1=seq(1, nX, nX %/% n_blocks)
  X_per_block1 = append(x = X_per_block1,values = nX)
  list_causalX=list()
  for (jl in 2:length(X_per_block1)){ list_causalX[[jl-1]] = seq(X_per_block1[jl-1],X_per_block1[jl]) }
  onelist1=which(unlist(lapply(list_causalX,length))==1)
  if (length(onelist1)!=0){list_causalX = list_causalX[-onelist1]}
  
  while (length(list_causalX) > length(list_causal)){
    list_causalX[[length(list_causalX)-1]] = append(list_causalX[[length(list_causalX)-1]],
                                                    list_causalX[[length(list_causalX)]])
    list_causalX[[length(list_causalX)]] =NULL
  }
  
  for (jm in 1:length(list_causalX)){ 
    gamma_SNP_X[,list_causalX[[jm]]] = matrix(rnorm(n = nG*length(list_causalX[[jm]]), 
                                                    mean = gamma_SNP_X_mu[[jm]], sd = 1),
                                              ncol = length(list_causalX[[jm]]))}
  
  X = matrix(nrow = n_ind,ncol = nX)
  for (lk in 1:length(list_causalX)){
    X[,list_causalX[[lk]]] = G[,list_causal[[lk]]] %*% gamma_SNP_X[list_causal[[lk]],list_causalX[[lk]]] + epsilon_X[,list_causalX[[lk]]]}
  
  X=apply(X,2, function(x) x+U)
  
  if (sum_stats==T){
    #Summary Statistics
    beta_matrix = SE_matrix = pval_matrix = data.frame(matrix(nrow = ncol(G), ncol = ncol(X)))
    pv_df=data.frame(matrix(nrow = ncol(G),ncol = ncol(X)))
    pvg=paste('p_val_GX',1:ncol(X),sep='')
    betag=paste('beta_GX',1:ncol(X),sep='')
    SEg=paste('SE_GX',1:ncol(X),sep='')
    
    for (k in 1:ncol(X)){
      lm_obj1=summary(lm(X[,k] ~ .,data = data.frame(G)))$coefficients
      p_val_temp = lm_obj1[-1,4]
      pv_df[,k] =p_val_temp; colnames(pv_df)[k]=pvg[k]
      beta_matrix[,k]=lm_obj1[-1,1]
      colnames(beta_matrix)[k]=betag[k]
      SE_matrix[,k]=lm_obj1[-1,2]
      colnames(SE_matrix)[k]=SEg[k]
    }
    rownames(beta_matrix)=rownames(SE_matrix)=rownames(pv_df)=paste('SNP',1:dim(G)[2],sep = '')
    b1=beta_matrix;se1=SE_matrix
    return(list(X=X,U=U,G=G,betaX=b1,seX=se1))
  }
  return(list(X=X,U=U,G=G))}

back to top