Template for GTEx example on cluster

original data

Please read your original data as a \(N \times P\) matrix in R.

In this case you can create a folder the restore the data as .rds file

mkdir testflashr
cd testflashr
cp gtexEQTL_zscore.rds GTEX/testflashr/

now you have the data matrix in gtexEQTL_zscore.rds

to restore the result and track the errors, you can crest folders

mkdir output
mkdir ourlog

creat a .R file

library(ashr)
library(flashr)
load("./gtexEQTL_zscore.rds")
Y = t(zscore)
data = flash_set_data(Y)
f_greedy = flash_add_greedy(data,Kmax=60)
f_greedy_bf = flash_backfit(data,f_greedy)
#f_greedy = flash_add_greedy(data,Kmax=60,var_type = "by_column",ash_param=list(method = "fdr"))
#f_greedy_bf = flash_backfit(data,f_greedy,var_type = "by_column",ash_param=list(method = "fdr"))
saveRDS(f_greedy,file = "./output/gflashvarcol.rds")
saveRDS(f_greedy_bf,file = "./output/bflashvarcol.rds")

here 60 is much larger than the sample size. our method doesn’t restrict \(K < \min(P,N)\). But for this case, 60 is enough. We call this file as flashwrapper.R

creat a .sbatch file

#!/bin/bash

#SBATCH --job-name=arrayJob
#SBATCH --output=/home/weidong/HG/flash/data/GTEX/testflashr/outlog/arrayJob_%A_%a.out
#SBATCH --error=/home/weidong/HG/flash/data/GTEX/testflashr/outlog/arrayJob_%A_%a.err
#SBATCH --array=1
#SBATCH --time=30:00:00
#SBATCH --partition=mstephens
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=1
#SBATCH --mem-per-cpu=8G


######################
# Begin work section #
######################

# Print this sub-job's task ID
cd /home/weidong/HG/flash/data/GTEX/testflashr
Rscript --verbose flashwrapper.R

you can substitute the folder path /home/weidong/HG/flash/data/GTEX/testflashr. and we call this .sbatch file as runflash.sbatch

sfa wrapper

this example is on PPS cluster

SFA.wrapper = function(X, K){
  N = dim(X)[1]
  P = dim(X)[2]

  if(file.exists("SFAout")){
    unlink("SFAout", recursive= T)
  }
  system("mkdir SFAout")

  write.table(X,file="dscsfa.txt",row.names=F,col.names=F)
  # sfa command line on cluster
  # ~/flash/simulation/methods/sfa/src/sfa -gen ./dscsfamix.txt -g 600 -k 1 -n 200 -iter 100 -rand 999 -o ./SFAout/sfa

  # sqrt(mean(( loadings %*% factor - test$L_true%*%t(test$F_true))^2))/sqrt(mean((test$L_true%*%t(test$F_true))^2))
  system(paste("~/flash/simulation/methods/sfa/src/sfa", "-gen", "./dscsfa.txt", "-g", N, "-k", K,
               "-n", P, "-iter", 100, "-rand", 999, "-o", "./SFAout/sfa", sep=" "))

  if(file.info("./SFAout/sfa_F.out")$size == 1){
    return(list(L = 0, F = 0))
  } else{
    Fhat=read.table("./SFAout/sfa_F.out")
    lambda=read.table("./SFAout/sfa_lambda.out")
    Fhat=as.matrix(Fhat)
    lambda=as.matrix(lambda)
    return(list(L = lambda, F = Fhat))
  }

}

SFAmix.wrapper = function(X,K){
  N = dim(X)[1]
  P = dim(X)[2]
  
  write.table(X,file="dscsfamix.txt",row.names=F,col.names=F)
  print(class(X))
  print(dim(X))
  
  if(file.exists("SFAmixout")){
    unlink("SFAmixout", recursive= T)
  }
  system("mkdir SFAmixout")
  
  # ~/mvash/mvsim/SFAmix/SFAmix --nf K --y dscrsfamix.txt --out SFAmixout --sep space
  # this is run on the PPS cluster
  system(paste("~/mvash/mvsim/SFAmix/SFAmix","--nf",K,"--y","dscsfamix.txt","--out",
               "SFAmixout","--sep","space",sep=" "))
  
  alpha=read.table("./SFAmixout/PSI")
  
  if(file.info("./SFAmixout/EX")$size == 1){
    return(list(L = 0, F = 0))
  } else{
    Fhat=read.table("./SFAmixout/EX")
    lambda=read.table("./SFAmixout/LAM")
    Psi=as.vector(alpha)
    Fhat=as.matrix(Fhat)
    lambda=as.matrix(lambda)
    P=dim(lambda)[2]
    n=dim(Fhat)[1]
    return(list(L = Fhat, F = lambda))
  }
  
}

this is sfawrapper.R

creat run_sfa.R file

load("~/flash/dataanalysis/GTeX/EQTLzscore/gtexEQTL_zscore.rds")
setwd("~/flash/dataanalysis/GTeX/EQTLzscore/")
source("./sfawrapper.R")
source("./SFAmixwrapper.R")
Y = zscore
N = dim(Y)[1]
P = dim(Y)[2]
gsfa = SFA.wrapper(Y, 30)
saveRDS(gsfa, file = "~/flash/dataanalysis/GTeX/EQTLzscore/sfaGTExzscore.rds")


Y = t(zscore)
gsfamix = SFAmix.wrapper(Y,30)
saveRDS(gsfamix, file = "~/flash/dataanalysis/GTeX/EQTLzscore/sfamixGTExzscore.rds")

and run it

R CMD BATCH run_sfa.R

> gsfa = SFA.wrapper(Y, 30)
 Setting G to be 16069
 Setting N to be 44
Trying to read in matrix with G=16069 and N =44
Opening matrix ./dscsfa.txt with 16069 rows and 44 columns
Opening matrix ./dscsfa.txt with 16069 rows and 44 columns
read in matrix 16069 by 44
initializing sfa...
G = 16069
N = 44
K = 30

SFA K=3

g3 = readRDS('../data/GTExdata/SFAres/sfaK_3.rds')
par(mfrow = c(1,3),mar=c(5.1,4.1,4.1,2.1)-1.9)
for(k in 1:3){
  barplot(g3$F[k,],main = paste("factor",k))
}

plot SFA result K = 30

sfares = readRDS("../data/GTExdata/SFAres/sfaGTExzscore.rds")
# the F is K by N matrix! not N by K matrix
par(mfrow = c(2,3),mar=c(5.1,4.1,4.1,2.1)-1.9)
for(k in 1:30){
  barplot(sfares$F[k,],main = paste("factor",k))
}

plot of SFA result with t(Y) K = 10

If you would like to use t(Y) you should use -vn in the sfa command line in order to get \((\hat{\sigma}_1^2,\cdots,\hat{\sigma}_j^2,\cdots)\) is a 16069 vector

you need add -mg to get the mean vector with length of 44

load("../data/GTExdata/gtexEQTL_zscore.rds")
gt10 = readRDS('../data/GTExdata/SFAres/sfaTK_10m.rds')
ssy = sum(zscore^2)
pve = sapply(seq(1:10),function(x){sum((gt10$LF$L[,x] %*% t(gt10$LF$F[x,]))^2)})
pve = pve/ssy

mu = sapply(seq(1:44),function(x){gt10$mu[[x]]})
par(mfrow = c(2,3),mar=c(5.1,4.1,4.1,2.1)-1.9)
barplot(mu,main = "mean vector")
for(k in 1:10){
  barplot(gt10$LF$L[,k],main = paste("factor",k,"pve",round(pve[k],3)))
}

gsl: lu.c:147: ERROR: matrix is singular
Default GSL error handler invoked.

plot SFAmix result

sfamixres = readRDS("../data/GTExdata/SFAres/sfamixGTExzscore.rds")
dim(sfamixres$L)

[1] 44  1

dim(sfamixres$F)

[1]     1 16069

this is only a rank one matrix. I tried two time. two runs give rank one matrix.