1 Convert text files to binary (bod) format

mkdir -p ${wd}/BodFiles

# Cortical files
for hemi in lh rh
do
for moda in area thickness
do
${bind}/qsubshcom " "${bind}"/osca --tefile "${wd}"/FS/"${hemi}"."${moda}".fwhm0.UKB.txt --make-bod --no-fid --out "${wd}"/BodFiles/"${hemi}"."${moda}".fwhm0.UKB --thread-num 5 " 5 30G ${hemi}_${moda}_intoBod 10:00:00 "-acct=UQ-IMB-CNSG" 
done
done 

# Subcortical
for hemi in lh rh
do
for moda in thick LogJacs
do
${bind}/qsubshcom " "${bind}"/osca --efile "${wd}"/FS/"${hemi}"."${moda}".transposed.txt --make-bod --no-fid --out "${wd}"/BodFiles/"${hemi}"."${moda}".fwhm0.UKB.vertexQC " 1 8G ${hemi}_${moda}_intoBod 10:00:00 "-acct=UQ-IMB-CNSG" 
done
done 

2 Exclude vertices outside of cortex

The VerticesToExclude files can be found in the atlas folder of the GitHub repository. Another more data-driven way is to exclude vertices based on their standard deviation (vertices outside of cortex have 0 values for most-to-all subjects) using the OSCA option –sd-min (https://cnsgenomics.com/software/osca/#ManageBODfile(s)).

for hemi in lh rh
do
for moda in thickness area
do
${bind}/qsubshcom " "${bind}"/osca --befile "${wd}"/BodFiles/"${hemi}"."${moda}".fwhm0.UKB --exclude-probe "${wd}"/VerticesToExclude_"${hemi}"_"${moda}"_Cortex.txt --make-bod --out "${wd}"/BodFiles/"${hemi}"."${moda}".fwhm0.UKB.vertexQC " 1 10G ${hemi}_${moda}_excludeVertices 10:00:00 "-acct=UQ-IMB-CNSG" 
done
done 

3 Calculate Brain Relatedness Matrices (BRM)

for hemi in lh rh
do
for moda in thickness area thick LogJacs
do
${bind}/qsubshcom " "${bind}"/osca --befile "${wd}"/BodFiles/"${hemi}"."${moda}".fwhm0.UKB.vertexQC --make-orm-bin --out "${wd}"/BodFiles/"${hemi}"."${moda}".fwhm0.UKB.BRM --thread-num 5 " 5 30G ${hemi}_${moda}_BRM_calculation 10:00:00 "-acct=UQ-IMB-CNSG" 
done
done

4 Extract mean and variance of each vertex

for hemi in lh rh
do
for moda in thickness area thick LogJacs
do
${bind}/qsubshcom " "${bind}"/osca --befile "${wd}"/BodFiles/"${hemi}"."${moda}".fwhm0.UKB.vertexQC --get-mean --get-variance --out "${wd}"/BodFiles/"${hemi}"."${moda}" " 1 8G ${hemi}_${moda}_mean_var 48:00:00 "-acct=UQ-IMB-CNSG" 
done
done 

5 Combine files across measurement types (for BLUP analyses)

for hemi in lh rh
do
for moda in area thickness thick LogJacs
do
${bind}/qsubshcom " "${bind}"/osca --befile "${wd}"/BodFiles/"${hemi}"."${moda}".fwhm0.UKB.vertexQC --std-probe --make-bod --out "${wd}"/BodFiles/${hemi}.${moda}.fwhm0.UKB.BLUP.STD " 1 8G ${hemi}_${moda}_STD_probes 02:00:00 "-acct=UQ-IMB-CNSG" 
done
done

${bind}/qsubshcom " "${bind}"/osca --befile-flist "${wd}"/mybod.flist --make-bod --out "${wd}"/BodFiles/UKB_allVertices " 1 8G BLUP_allVertices 02:00:00 "-acct=UQ-IMB-CNSG" 

6 Merge BRMs

You will need a text file mergeBRM_All.txt which contains the path/names of the BRM you want to merge

${bind}/qsubshcom " "${bind}"/osca --merge-orm "$wd"/mergeBRM_All.txt --make-orm --out "${wd}"/BodFiles/allVerticesBRM " 1 8G CombineBRMs 02:00:00 "-acct=UQ-IMB-CNSG" 

7 Perform QC using the BRMs

This section uses a R script to open the BRMs, produce plots and calculate BRM off-diagonal cut-off to use in QC.

# Launch interactive session (syntax depends on your HPC, here using a PBSpro)
qsub -I -X -l select=1:ncpus=1:mem=30GB -l walltime=05:00:00 -A UQ-IMB-CNSG

wd="working/directory"
cd $wd
module load R
R

8 Example of running the R script on two BRMs

source("RR_4.0_BRM_QC_functions.R")

# Get variance of diagonal and off-diagonal elements
vals=NULL
IDs2ExDiag=NULL
  for (BRM in c("allVerticesBRM", "allVerticesBRM_QC")){
  brm=asBRM(ReadORMBin(paste0( moda)))
  vals=rbind(vals, c(mean(diag(brm)),var(diag(brm)), mean(brm[upper.tri(brm)]),var(brm[upper.tri(brm)]) , sqrt(var(brm[upper.tri(brm)]))*5 , min(brm[upper.tri(brm)]), max(brm[upper.tri(brm)]) ) )
  print(moda)
  png(paste0("HistogramBRM_", BRM, ".png"), width = 20, height = 10, units = "cm", res = 400)
  par(mfrow=c(1,2))
  hist(diag(brm), main = paste0( BRM, " BRM diagonal"), breaks=500 )
  hist(brm[upper.tri(brm)], main = paste0(BRM, " BRM diagonal"), breaks=5000 )
  dev.off()

}
rownames(vals)=c("allVerticesBRM", "allVerticesBRM_QC")
colnames(vals)=c("Mean - diag", "Var - diag", "Mean - off diag", "Var - off diag", "BRM cutoff", "min off-diag", "max off-diag")
print(vals)

# FYI : BRM description from our "replication UKB sample, N about 5K"
# I found that those values were rather stable across samples

#          Mean - diag  Var - diag Mean - off diag Var - off diag BRM cutoff     min off-diag max off-diag
#allVerticesBRM      1.0000000 0.020123953   -0.0002192982   0.0017715966  0.2104517  -0.4211293    0.6721468
#allVerticesBRM_QC   0.9654818 0.009724773    0.0005558369   0.0008534802  0.1460719  -0.2206390    0.2027178

#See also QC_BRM_outliers() function for QC used in the discovery sample

9 Enforce cutoff on off-diagonal elements

Note that in the replication sample we used a cutoff a 5SD from mean but this results in a quite stringent QC.

${bind}/qsubshcom " "${bind}"/osca --reml --orm "${wd}"/allVerticesBRM --orm-cutoff 0.21 --orm-cutoff-2sides --make-orm --out "${wd}"/allVerticesBRM_QC " 1 4G BRM_QC 4:00:00 "-acct=UQ-IMB-CNSG"