1 Apply smoothing and coarses mesh

This can be integrated in the MRI images processing script (RR_2), it relies on FreeSurfer scripts and files created as part of the FreeSurfer processing (FSresults folder)

for batch in {1..5}
do
${bind}/qsubshcom " echo \${TASK_ID}.\$(sed -n \"\${TASK_ID}{p;q}\" "${medici}"/T1/raw/batch"${batch}"/T1_fetch_batch"${batch}".lis) |;
file=\$(sed -n \"\${TASK_ID}{p;q}\" "${medici}"/T1/raw/batch"${batch}"/T1_fetch_batch"${batch}".lis) |;
ID=\${file:0:7} |;
echo \${file} |;
echo \${ID} |;
for hemi in lh rh |;
do |;
for moda in area thickness |;
do |;
for fwhm in 5 10 15 20 25 |; 
do |;
echo \${hemi}.\${moda}.\${fwhm} |;
"${fsdir}"/bin/mris_convert -c "${wd}"/FS/FSresults/batch"${batch}"/\${ID}/\${hemi}.\${moda}.fwhm\${fwhm}.fsaverage.mgh "${fsdir}"/subjects/fsaverage/surf/\${hemi}.orig "${wd}"/FS/FSresults/batch"${batch}"/\${ID}/\${hemi}.\${moda}.fwhm\${fwhm}.fsaverage.asc |;
done |;
for fsav in fsaverage3 fsaverage4 fsaverage5 fsaverage6 |;
do |;
echo \${hemi}.\${moda}.\${fsav} |;
"${fsdir}"/bin/mri_surf2surf --s fsaverage --hemi \${hemi} --sval "${wd}"/FS/FSresults/batch"${batch}"/\${ID}/\${hemi}.\${moda}.fwhm0.fsaverage.mgh --trgsubject \${fsav} --tval "${wd}"/FS/FSresults/batch"${batch}"/\${ID}/\${hemi}.\${moda}.fwhm0.\${fsav}.mgh |;
"${fsdir}"/bin/mris_convert -c "${wd}"/FS/FSresults/batch"${batch}"/\${ID}/\${hemi}.\${moda}.fwhm0.\${fsav}.mgh "${fsdir}"/subjects/\${fsav}/surf/\${hemi}.orig "${wd}"/FS/FSresults/batch"${batch}"/\${ID}/\${hemi}.\${moda}.fwhm0.\${fsav}.asc |;
done |;
done |;
done |;
" 1 4G ExtractASC_${batch} 48:00:00 "-array=1-500 -acct=UQ-IMB-CNSG" 
done

2 Create the tables of vertex-wise data

Similar to code in RR_2, except that we initialise the vertex names using file the first individual of each batch. I have included a check that processing completed for this subject, if not the job fails immediately and you will have to provide manually an ID.
This is better than the previous version (RR_2 that required a manual input) but still not ideal

for batch in {1..5}
do
for hemi in lh rh
do
for moda in thickness area
do
for fwhm in 5 10 15 20 25
do
echo ${hemi}.${moda}
${bind}/qsubshcom " firstID=\$(ls "${wd}"/FS/FSresults/batch"${batch}"/ | sort -n | tail -1) |;
echo \${firstID} |;
ls "${wd}"/FS/FSresults/batch"${batch}"/\${firstID}/*"${hemi}"."${moda}".fwhm"${fwhm}".fsaverage*.asc |;
if [ -f "${wd}"/FS/FSresults/batch"${batch}"/\${firstID}/"${hemi}"."${moda}".fwhm"${fwhm}".fsaverage.asc ] |; 
then |;
echo \"vertexnum\" > "${wd}"/FS/"${hemi}"."${moda}".batch"${batch}".fwhm"${fwhm}".UKB.txt |; 
awk '{print \$1}' "${wd}"/FS/FSresults/batch"${batch}"/\${firstID}/"${hemi}"."${moda}".fwhm"${fwhm}".fsaverage.asc >> "${wd}"/FS/"${hemi}"."${moda}".batch"${batch}".fwhm"${fwhm}".UKB.txt |;
for ID in \$(awk -F\",\" \"NR>0 {print \$1}\" "${medici}"/T1/raw/batch"${batch}"/T1_fetch_batch"${batch}".lis) |;
do |;
ID=\${ID:0:7} |;
echo \${ID} |;
if [ -f "${wd}"/FS/FSresults/batch"${batch}"/\${ID}/"${hemi}"."${moda}".fwhm"${fwhm}".fsaverage.asc ] |;
then |;
echo \${ID} > "${wd}"/FS/"${hemi}"."${moda}"."${batch}".fwhm"${fwhm}".temp.lta |;
awk '{print \$5}' "${wd}"/FS/FSresults/batch"${batch}"/\${ID}/"${hemi}"."${moda}".fwhm"${fwhm}".fsaverage.asc >> "${wd}"/FS/"${hemi}"."${moda}"."${batch}".fwhm"${fwhm}".temp.lta |;
paste "${wd}"/FS/"${hemi}"."${moda}".batch"${batch}".fwhm"${fwhm}".UKB.txt "${wd}"/FS/"${hemi}"."${moda}"."${batch}".fwhm"${fwhm}".temp.lta > "${wd}"/FS/"${hemi}"."${moda}"."${batch}".fwhm"${fwhm}".temp2.lta |;
cp "${wd}"/FS/"${hemi}"."${moda}"."${batch}".fwhm"${fwhm}".temp2.lta "${wd}"/FS/"${hemi}"."${moda}".batch"${batch}".fwhm"${fwhm}".UKB.txt |;
fi |; done |; fi " 1 8G CortDat_${batch}_${hemi}_${moda}_${fwhm} 10:00:00 " -acct=UQ-IMB-CNSG"
done
done
done
done

for batch in {1..5}  
do
for hemi in lh rh
do
for moda in area thickness
do
for fsav in fsaverage3 fsaverage4 fsaverage5 fsaverage6
do
echo ${batch}.${hemi}.${moda}.${fsav}
${bind}/qsubshcom " firstID=\$(ls "${wd}"/FS/FSresults/batch"${batch}"/ | sort -n | tail -1) |;
echo \${firstID} |;
ls "${wd}"/FS/FSresults/batch"${batch}"/\${firstID}/*"${hemi}"."${moda}".fwhm0."${fsav}"*.asc |;
if [ -f "${wd}"/FS/FSresults/batch"${batch}"/\${firstID}/"${hemi}"."${moda}".fwhm0."${fsav}".asc ] |; 
then |;
echo \"vertexnum\" > "${wd}"/FS/"${hemi}"."${moda}".batch"${batch}".fwhm0."${fsav}".UKB.txt |; 
awk '{print \$1}' "${wd}"/FS/FSresults/batch"${batch}"/\${firstID}/"${hemi}"."${moda}".fwhm0."${fsav}".asc >> "${wd}"/FS/"${hemi}"."${moda}".batch"${batch}".fwhm0."${fsav}".UKB.txt |;
for ID in \$(awk -F\",\" \"NR>0 {print \$1}\" "${medici}"/T1/raw/batch"${batch}"/T1_fetch_batch"${batch}".lis) |;
do |;
ID=\${ID:0:7} |;
echo \${ID} |;
if [ -f "${wd}"/FS/FSresults/batch"${batch}"/\${ID}/"${hemi}"."${moda}".fwhm0."${fsav}".asc ] |;
then |;
echo \${ID} > "${wd}"/FS/"${hemi}"."${moda}"."${batch}".fwhm0."${fsav}".temp.lta |;
awk '{print \$5}' "${wd}"/FS/FSresults/batch"${batch}"/\${ID}/"${hemi}"."${moda}".fwhm0."${fsav}".asc >> "${wd}"/FS/"${hemi}"."${moda}"."${batch}".fwhm0."${fsav}".temp.lta |;
paste "${wd}"/FS/"${hemi}"."${moda}".batch"${batch}".fwhm0."${fsav}".UKB.txt "${wd}"/FS/"${hemi}"."${moda}"."${batch}".fwhm0."${fsav}".temp.lta > "${wd}"/FS/"${hemi}"."${moda}"."${batch}".fwhm0."${fsav}".temp2.lta |;
cp "${wd}"/FS/"${hemi}"."${moda}"."${batch}".fwhm0."${fsav}".temp2.lta "${wd}"/FS/"${hemi}"."${moda}".batch"${batch}".fwhm0."${fsav}".UKB.txt |;
fi |; done |; fi |;
rm -r "${wd}"/FS/*.lta " 1 8G CortDat_${batch}_${hemi}_${moda}_${fsav} 48:00:00 " -acct=UQ-IMB-CNSG"
done
done
done
done

3 Check dimension of files created

for batch in {1..5}
do
for hemi in lh rh
do
for moda in area thickness
do
for fwhm in 5 10 15 20 25
do
echo batch$batch.$hemi.$moda.fwhm$fwhm
awk '{print NF}' ${wd}/FS/$hemi.$moda.batch$batch.fwhm$fwhm.UKB.txt | sort -nu | tail -n 1
cat ${wd}/FS/$hemi.$moda.batch$batch.fwhm$fwhm.UKB.txt | wc -l
done
for fsav in fsaverage3 fsaverage4 fsaverage5 fsaverage6
do
echo batch$batch.$hemi.$moda.$fsav
awk '{print NF}' ${wd}/FS/$hemi.$moda.batch$batch.fwhm0.$fsav.UKB.txt | sort -nu | tail -n 1
cat ${wd}/FS/$hemi.$moda.batch$batch.fwhm0.$fsav.UKB.txt | wc -l
done
done
done 
done

4 Give unique names to vertices and create tab files combining all batches

for fwhm in 5 10 15 20 25
do
${bind}/qsubshcom " awk 'BEGIN{OFS=\"\t\"}\$1=\"rht_\"\$1' "${wd}"/FS/rh.thickness.batch1.fwhm"${fwhm}".UKB.txt > "${wd}"/FS/rh.thickness.fwhm"${fwhm}".UKB.txt |;
awk 'BEGIN{OFS=\"\t\"}\$1=\"lht_\"\$1' "${wd}"/FS/lh.thickness.batch1.fwhm"${fwhm}".UKB.txt > "${wd}"/FS/lh.thickness.fwhm${fwhm}.UKB.txt |;
awk 'BEGIN{OFS=\"\t\"}\$1=\"rha_\"\$1' "${wd}"/FS/rh.area.batch1.fwhm"${fwhm}".UKB.txt > "${wd}"/FS/rh.area.fwhm"${fwhm}".UKB.txt |;
awk 'BEGIN{OFS=\"\t\"}\$1=\"lha_\"\$1' "${wd}"/FS/lh.area.batch1.fwhm"${fwhm}".UKB.txt > "${wd}"/FS/lh.area.fwhm"${fwhm}".UKB.txt |;
for batch in {2..10} |;
do |;
for moda in area thickness |;
do |;
for hemi in rh lh |;
do  |;
echo \${batch}.\${hemi}.\${moda}."${fwhm}" |;
awk '{\$1=\"\"}1' "${wd}"/FS/\${hemi}.\${moda}.batch\${batch}.fwhm"${fwhm}".UKB.txt | awk '{\$1=\$1}1' > "${wd}"/FS/\${hemi}.\${moda}.batch\${batch}.fwhm"${fwhm}".UKB_noColNames.temp |;
paste "${wd}"/FS/\${hemi}.\${moda}.fwhm"${fwhm}".UKB.txt "${wd}"/FS/\${hemi}.\${moda}.batch\${batch}.fwhm"${fwhm}".UKB_noColNames.temp > "${wd}"/FS/\${hemi}.\${moda}.fwhm"${fwhm}".UKB_2.temp |;
cp "${wd}"/FS/\${hemi}.\${moda}.fwhm"${fwhm}".UKB_2.temp "${wd}"/FS/\${hemi}.\${moda}.fwhm"${fwhm}".UKB.txt |;
done |;
done |;
done  |; " 1 10G PutVertexBatchesTogether_${fwhm} 10:00:00 "-acct=UQ-IMB-CNSG" 
done

for fsav in fsaverage3 fsaverage4 fsaverage5 fsaverage6
do
awk 'BEGIN{OFS="\t"}$1="rht_"$1' ${wd}/FS/rh.thickness.batch1.fwhm0.${fsav}.UKB.txt > ${wd}/FS/rh.thickness.fwhm0.${fsav}.UKB.txt
awk 'BEGIN{OFS="\t"}$1="lht_"$1' ${wd}/FS/lh.thickness.batch1.fwhm0.${fsav}.UKB.txt > ${wd}/FS/lh.thickness.fwhm0.${fsav}.UKB.txt
awk 'BEGIN{OFS="\t"}$1="rha_"$1' ${wd}/FS/rh.area.batch1.fwhm0.${fsav}.UKB.txt > ${wd}/FS/rh.area.fwhm0.${fsav}.UKB.txt
awk 'BEGIN{OFS="\t"}$1="lha_"$1' ${wd}/FS/lh.area.batch1.fwhm0.${fsav}.UKB.txt > ${wd}/FS/lh.area.fwhm0.${fsav}.UKB.txt

for batch in {2..10}
do
# Get rid of vertices (row names) of files  
for moda in thickness area
do
for hemi in lh rh
do 
echo ${batch}.${hemi}.${moda}.${fsav}
# Get rid row names (1st column)
awk '{$1=""}1' ${wd}/FS/${hemi}.${moda}.batch${batch}.fwhm0.${fsav}.UKB.txt | awk '{$1=$1}1' > ${wd}/FS/${hemi}.${moda}.batch${batch}.fwhm0.${fsav}.UKB_noColNames.temp
# Paste to the combined file        
paste ${wd}/FS/${hemi}.${moda}.fwhm0.${fsav}.UKB.txt ${wd}/FS/${hemi}.${moda}.batch${batch}.fwhm0.${fsav}.UKB_noColNames.temp > ${wd}/FS/${hemi}.${moda}.fwhm0.${fsav}.UKB_2.temp
cp ${wd}/FS/${hemi}.${moda}.fwhm0.${fsav}.UKB_2.temp ${wd}/FS/${hemi}.${moda}.fwhm0.${fsav}.UKB.txt
done
done 
done
done

# ------------ Tidy up
# rm $wd/FS/*temp

5 Convert all combined files into BOD (binary) format

${bind}/qsubshcom " for hemi in rh lh |;
do |;
for moda in thickness area |;
do |;
for fwhm in 5 10 15 20 25 |;
do |;
"${bind}"/osca --tefile "${wd}"/FS/\${hemi}.\${moda}.fwhm\${fwhm}.UKB.txt  --make-bod --no-fid --out "${wd}"/BodFiles/\${hemi}.\${moda}.fwhm\${fwhm}.UKB --thread-num 5 |;
done |;
for fsav in fsaverage3 fsaverage4 fsaverage5 fsaverage6 |;
do |;
"${bind}"/osca --tefile "${wd}"/FS/\${hemi}.\${moda}.fwhm0.\${fsav}.UKB.txt  --make-bod --no-fid --out "${wd}"/BodFiles/\${hemi}.\${moda}.fwhm0.\${fsav}.UKB --thread-num 5 |;
done |;
done |;
done |; " 5 20G Files_intoBod 10:00:00 "-acct=UQ-IMB-CNSG" 

6 Exclude vertices based on variance

This is a more data driven approach to excluding vertices outside of cortex for which the values are 0 for most/all participants.
If you look at the histogram of SD of vertices, you will see a peak in 0, which corresponds to the vertices we want to remove.

See OSCA options –get-variance and –get-mean and example in RR_4 section.

${bind}/qsubshcom " for hemi in rh lh |;
do |;
for moda in thickness area |;
do |;
for fwhm in 5 10 15 20 25 |;
do |;
"${bind}"/osca --befile "${wd}"/BodFiles/\${hemi}.\${moda}.fwhm\${fwhm}.UKB --sd-min 0.001 --make-bod --out "${wd}"/BodFiles/\${hemi}.\${moda}.fwhm\${fwhm}.UKB.QC |;
done |;
for fsav in fsaverage3 fsaverage4 fsaverage5 fsaverage6 |;
do |;
"${bind}"/osca --befile "${wd}"/BodFiles/\${hemi}.\${moda}.fwhm0.\${fsav}.UKB --sd-min 0.001 --make-bod --out "${wd}"/BodFiles/\${hemi}.\${moda}.fwhm0.\${fsav}.UKB.QC |;
done |;
done |;
done |; " 1 20G Vertex_QC_SD 10:00:00 "-acct=UQ-IMB-CNSG" 

7 Create BRMs

for hemi in rh lh
do
for moda in thickness area
do 
${bind}/qsubshcom "for fwhm in 5 10 15 20 25 |;
do |;
"${bind}"/osca --befile "${wd}"/BodFiles/"${hemi}"."${moda}".fwhm\${fwhm}.UKB.QC  --make-orm --out "${wd}"/BodFiles/"${hemi}"."${moda}".fwhm\${fwhm}.UKB --thread-num 2 |;
done |;
for fsav in fsaverage3 fsaverage4 fsaverage5 fsaverage6 |;
do |;
"${bind}"/osca --befile "${wd}"/BodFiles/"${hemi}"."${moda}".fwhm0.\${fsav}.UKB.QC --make-orm --out "${wd}"/BodFiles/"${hemi}"."${moda}".fwhm0.\${fsav}.UKB --thread-num 2 |;
done  |; " 2 15G BRM_${hemi}_${moda} 10:00:00 "-acct=UQ-IMB-CNSG" 
done 
done 

# Check number of participants - check length of orm.id file (ID list)
for hemi in lh rh
do
for moda in thickness area
do 
for fwhm in 5 10 15 20 25
do
echo $hemi.$moda.fwhm$fwhm
cat $wd/$hemi.$moda.fwhm$fwhm.UKB.orm.id | wc -l
done
done
done


for hemi in lh rh
do
for moda in thickness area
do 
for fsav in fsaverage3 fsaverage4 fsaverage5 fsaverage6
do
echo $hemi.$moda.fwhm0.$fsav
cat $wd/$hemi.$moda.fwhm0.$fsav.UKB.orm.id | wc -l
done
done
done

8 Combine BRMs across modalities and add subcortical BRM data

Here, we create the .flist files (which contain list of BRM to combine) automatically

for fsav in fsaverage3 fsaverage4 fsaverage5 fsaverage6
do
echo ${fsav}
echo "/working/directory/BodFiles/lh.thick.fwhm0.UKB.BRM" > $wd/allVerticesBRM.${fsav}.flist
echo "/working/directory/BodFiles/lh.LogJacs.fwhm0.UKB.BRM" >> $wd/allVerticesBRM.${fsav}.flist
echo "/working/directory/BodFiles/rh.thick.fwhm0.UKB.BRM" >> $wd/allVerticesBRM.${fsav}.flist
echo "/working/directory/BodFiles/rh.LogJacs.fwhm0.UKB.BRM" >> $wd/allVerticesBRM.${fsav}.flist
for hemi in lh rh
do
for moda in thickness area
do
echo ${wd}/${hemi}.${moda}.fwhm0.${fsav}.UKB >> ${wd}/allVerticesBRM.${fsav}.flist
done
done
done

for fwhm in 5 10 15 20 25
do
echo ${fwhm}
echo "/working/directory/BodFiles/lh.thick.fwhm0.UKB.BRM" > $wd/allVerticesBRM.fwhm${fwhm}.flist
echo "/working/directory/BodFiles/lh.LogJacs.fwhm0.UKB.BRM" >> $wd/allVerticesBRM.fwhm${fwhm}.flist
echo "/working/directory/BodFiles/rh.thick.fwhm0.UKB.BRM" >> $wd/allVerticesBRM.fwhm${fwhm}.flist
echo "/working/directory/BodFiles/rh.LogJacs.fwhm0.UKB.BRM" >> $wd/allVerticesBRM.fwhm${fwhm}.flist
for hemi in lh rh
do
for moda in thickness area
do
echo ${wd}/${hemi}.${moda}.fwhm${fwhm}.UKB >> ${wd}/allVerticesBRM.fwhm${fwhm}.flist
done
done
done

9 Run osca to create combined BRMs

for fwhm in 5 10 15 20 25
do
${bind}/osca --merge-orm ${wd}/BodFiles/allVerticesBRM.fwhm${fwhm}.flist --make-orm --out ${wd}/BodFiles/allVerticesBRM.fwhm${fwhm}
done

for fsav in fsaverage3 fsaverage4 fsaverage5 fsaverage6
do
${bind}/osca --merge-orm ${wd}/BodFiles/allVerticesBRM.${fsav}.flist --make-orm --out ${wd}/BodFiles/allVerticesBRM.${fsav}
done

10 Run QC on BRM

setwd("working/directory/BodFiles")
source("RR_4.0_BRM_QC_functions.R")

# Get variance of diagonal and off-diagonal elements
vals=NULL
IDs2ExDiag=NULL
  for (moda in c("allVerticesBRM.fwhm5", "allVerticesBRM.fwhm10", "allVerticesBRM.fwhm15", "allVerticesBRM.fwhm20", "allVerticesBRM.fwhm25", "allVerticesBRM.fsaverage3", "allVerticesBRM.fsaverage4", "allVerticesBRM.fsaverage5", "allVerticesBRM.fsaverage6")){
  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_", moda, ".png"), width = 20, height = 10, units = "cm", res = 400)
  par(mfrow=c(1,2))
  hist(diag(brm), main = paste0( moda, " BRM diagonal"), breaks=500 )
  hist(brm[upper.tri(brm)], main = paste0(moda, " BRM diagonal"), breaks=5000 )
  dev.off()

}
rownames(vals)=c("allVerticesBRM.fwhm5", "allVerticesBRM.fwhm10", "allVerticesBRM.fwhm15", "allVerticesBRM.fwhm20", "allVerticesBRM.fwhm25", "allVerticesBRM.fsaverage3", "allVerticesBRM.fsaverage4", "allVerticesBRM.fsaverage5", "allVerticesBRM.fsaverage6")
colnames(vals)=c("Mean - diag", "Var - diag", "Mean - off diag", "Var - off diag", "BRM cutoff", "min off-diag", "max off-diag")
print(vals)

write.csv(vals, "Summary_BRM_QC.csv")

${bind}/qsubshcom " 
"${bind}"/osca --reml --orm "${wd}"/BodFiles/allVerticesBRM.fwhm5 --orm-cutoff 0.62 --orm-cutoff-2sides --make-orm --out "${wd}"/BodFiles/allVerticesBRM.fwhm5.QC |;
"${bind}"/osca --reml --orm "${wd}"/BodFiles/allVerticesBRM.fwhm10 --orm-cutoff 0.85 --orm-cutoff-2sides --make-orm --out "${wd}"/BodFiles/allVerticesBRM.fwhm10.QC |;
"${bind}"/osca --reml --orm "${wd}"/BodFiles/allVerticesBRM.fwhm15 --orm-cutoff 1.05 --orm-cutoff-2sides --make-orm --out "${wd}"/BodFiles/allVerticesBRM.fwhm15.QC |;
"${bind}"/osca --reml --orm "${wd}"/BodFiles/allVerticesBRM.fwhm20 --orm-cutoff 1.24 --orm-cutoff-2sides --make-orm --out "${wd}"/BodFiles/allVerticesBRM.fwhm20.QC |;
"${bind}"/osca --reml --orm "${wd}"/BodFiles/allVerticesBRM.fwhm25 --orm-cutoff 1.41 --orm-cutoff-2sides --make-orm --out "${wd}"/BodFiles/allVerticesBRM.fwhm25.QC |;
" 1 4G BRM.QC.fwhm 4:00:00 "-acct=UQ-IMB-CNSG" 


${bind}/qsubshcom " 
"${bind}"/osca --reml --orm "${wd}"/BodFiles/allVerticesBRM.fsaverage3 --orm-cutoff 1.06 --orm-cutoff-2sides --make-orm --out "${wd}"/BodFiles/allVerticesBRM.fsaverage3.QC |;
"${bind}"/osca --reml --orm "${wd}"/BodFiles/allVerticesBRM.fsaverage4 --orm-cutoff 0.97 --orm-cutoff-2sides --make-orm --out "${wd}"/BodFiles/allVerticesBRM.fsaverage4.QC |;
"${bind}"/osca --reml --orm "${wd}"/BodFiles/allVerticesBRM.fsaverage5 --orm-cutoff 0.74 --orm-cutoff-2sides --make-orm --out "${wd}"/BodFiles/allVerticesBRM.fsaverage5.QC |;
"${bind}"/osca --reml --orm "${wd}"/BodFiles/allVerticesBRM.fsaverage6 --orm-cutoff 0.42 --orm-cutoff-2sides --make-orm --out "${wd}"/BodFiles/allVerticesBRM.fsaverage6.QC |;
" 1 4G BRM.QC.fsavg 4:00:00 "-acct=UQ-IMB-CNSG" 

source("RR_4.0_BRM_QC_functions.R")

# Get variance of diagonal and off-diagonal elements
vals=NULL
IDs2ExDiag=NULL
  for (moda in paste0(c("allVerticesBRM.fwhm5", "allVerticesBRM.fwhm10", "allVerticesBRM.fwhm15", "allVerticesBRM.fwhm20", "allVerticesBRM.fwhm25", "allVerticesBRM.fsaverage3", "allVerticesBRM.fsaverage4", "allVerticesBRM.fsaverage5", "allVerticesBRM.fsaverage6"), ".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_", moda, ".png"), width = 20, height = 10, units = "cm", res = 400)
  par(mfrow=c(1,2))
  hist(diag(brm), main = paste0( moda, " BRM diagonal"), breaks=500 )
  hist(brm[upper.tri(brm)], main = paste0(moda, " BRM diagonal"), breaks=5000 )
  dev.off()

}
rownames(vals)=paste0(c("allVerticesBRM.fwhm5", "allVerticesBRM.fwhm10", "allVerticesBRM.fwhm15", "allVerticesBRM.fwhm20", "allVerticesBRM.fwhm25", "allVerticesBRM.fsaverage3", "allVerticesBRM.fsaverage4", "allVerticesBRM.fsaverage5", "allVerticesBRM.fsaverage6"), ".QC") 
colnames(vals)=c("Mean - diag", "Var - diag", "Mean - off diag", "Var - off diag", "BRM cutoff", "min off-diag", "max off-diag")
print(vals)

write.csv(vals, "Summary_BRM_afterQC.csv")