Analyze dosages data from Ricopili

analyze_dosages.pbs

#PBS -lnodes=1
#PBS -lwalltime=0:05:00

#!/bin/bash

while getopts a:b:o:p:d: option
do
  case "${option}"
    in
      a) phenotype=${OPTARG};;
      b) covariate=${OPTARG};;
      o) outputfile=${OPTARG};;
      p) probdir=${OPTARG};;
      d) outdir=${OPTARG};;
    esac
done

 
 #Write the start and stop points of the file
 jstart=$((($PBS_ARRAYID-1)*$nodeuse +1))
 jstop=$(($PBS_ARRAYID*$nodeuse))

 for j in $(seq -w $jstart 1 $jstop)
 do
  file_use=$(awk -v lineno=$j '{if(NR==lineno) print}' $outputfile)
  # zcat "$probdir"/"$file_use" | awk -v s=$nsub '{ printf $1 "," $2 "," $3; for(i=1; i<=s; i++) printf "," $(i*2+2)*2+$(i*2+3); printf "\n" }' | gzip > "$outdir"/"$file_use".doscnt.gz & #Convert file to a single entry dosage format
  $plink_location --fam "$probdir"/"$file_use".fam --dosage "$probdir"/"$file_use".gz format=2 --pheno $phenotype --covar $covariate --out "$outdir"/$file
  
 done
wait

analyze_dosages.sh

#Bash version, no job script neded 
cut -d " " -f1,2,4 /home/maihofer/vets/qc/pca/pts_vets_mix_am-qc-eur_pca.menv.mds_cov > ea_mds.cov

PHENO=pts_vets_mix_am-qc.pheno
cov=ea_mds.cov

#Specify where probability format genotypes are stored
probdir=/home/maihofer/vets/qc/imputation/dasuqc1_pts_vets_mix_am-qc.hg19.ch.fl/qc1

#Specify where output will be stored 
outdir=/home/maihofer/vets/qc/infoscores

#List dosages files
ls $probdir | grep .gz$ | sed 's/.gz//g' > files_to_analyze.txt

for file in $(head -n1 files_to_analyze.txt)
do
 $plink_location --fam "$probdir"/"$file".fam --dosage "$probdir"/"$file".gz format=2 --pheno $PHENO --covar $cov --out "$outdir"/$file
done

analyze_dosages_torque.sh

#Job script version
cut -d " " -f1,2,4 /home/maihofer/vets/qc/pca/pts_vets_mix_am-qc-eur_pca.menv.mds_cov > ea_mds.cov

PHENO=pts_vets_mix_am-qc.pheno
cov=ea_mds.cov

#Specify where probability format genotypes are stored
probdir=/home/maihofer/vets/qc/imputation/dasuqc1_pts_vets_mix_am-qc.hg19.ch.fl/qc1

#Specify where output will be stored 
outdir=/home/maihofer/vets/qc/infoscores

#List dosages files
ls $probdir | grep .gz$ | sed 's/.gz//g' > files_to_analyze.txt

qsub analyze_dosages.pbs  -t 1-$totjobs -d $workingdir -F "-a $PHENO -b $cov -d $outdir -p $probdir -n $nodeuse -o files_to_analyze.txt"