imputationDensity_workflow.sh

#!/bin/bash

# run as: bash imputationDensity_Workflow.sh -f <plink_filename> -s <species> -d <low_density_array> -n <sample_size> -l <ld_size> -o <outdir> -c <config>
# <plink_filename>: Plink name without .ped/.map extension
# <species> species identifier to be used in Plink
# <low_density_array>: list of SNP names from the desired low density SNP array 
# <sample_size>: n. of individuals that will be sampled initialliy from the original file 
# <ld_size>: n. of individuals that will be assigned the low density SNP array 
# <outdir> root output directory
# <config> path to config file with parameters
# Use -gt 1 to consume two arguments per pass in
# the loop (e.g. each
# argument has a corresponding value to go with it).
# Use -gt 0 to consume one or more arguments per pass in the loop (e.g.
# some arguments don't have a corresponding value to go with it such
# as in the --default example).

Help()
{
   # Display Help
   echo "This script runs the pipeline to measure the imputation accuracy from low to high density SNP arrays"
   echo "Paths to required software packages (e.g. R, Plink, Java, Beagle) are to be set in the file pathNames.txt"
   echo
   echo "Syntax: imputationDensity_Workflow.sh [-h|f|s|d|n|l|o|c]"
   echo "options:"
   echo "h     print this help"
   echo "f     Plink filename, binary (path to) [required]"
   echo "s     species [required]"
   echo "d     low density array (list of SNP names) [required]"
   echo "n     sample size (to be sampled randomly from the dataset) [required]"
   echo "l     n. of samples assigned to the low density SNP array [required]"
   echo "o     output directory [required]"
   echo "c     path to config file [if not provided, default is used]"
   echo
}

# Get the options
while getopts ":h" option; do
   case $option in
      h) # display Help
         Help
         exit;;
   esac
done

while [[ $# -gt 1 ]]
do
key="$1"

case $key in
    -f|--input_file)
    INPUTFILE="$2"
    shift # past argument
    ;;
    -s|--species)
    SPECIES="$2"
    shift
    ;;
    -d|--low_density_array)
    LOWDENSITY="$2"
    shift # past argument
    ;;
    -n|--sample_size)
    SAMPLESIZE="$2"
    shift # past argument
    ;;
    -l|--ld_size)
    LDSIZE="$2"
    shift # past argument
    ;;
    -o|--outdir)
    OUTDIR="$2"
    shift # past argument
    ;;
    -c|--config)
    CONFIG="$2"
    shift
    ;;
    *)
            # unknown option
    ;;
esac
shift # past argument or value
done

echo "########################################################"
echo "## WRAPPER BASH SCRIPT TO RUN THE IMPUTATION EXPERIMENTS"
echo "########################################################"

echo INPUT FILE  = "${INPUTFILE}"
echo SPECIES  = "${SPECIES}"
echo LOW DENSITY FILE     = "${LOWDENSITY}"
echo SAMPLE SIZE     = "${SAMPLESIZE}"
echo LD SAMPLE SIZE     = "${LDSIZE}"
echo OUT FOLDER     = "${OUTDIR}"
configFile="${CONFIG:-pathNames.txt}"
echo CONFIG FILE = "${configFile}"

if [[ -n $1 ]]; then
    echo "Last line of file specified as non-opt/last argument:"
    tail -1 $1
fi

## CHECK IF INPUT FILE EXISTS
if [ -f "${INPUTFILE}.bed" ]; then
    echo "$INPUTFILE exists."
else
	echo "Input file does not exist - TERMINATING"; 
	exit 1; 
fi 

cwd=`pwd`
echo "current directory is $cwd"
source $configFile

### Parameters from the config file
echo "#######################"
echo "## EXPERIMENT         #"
echo "#######################"

echo "Experiment type: ${PREFIX}"

echo "################################################"
echo "## MAIN PATHS TO SOFTWARE - from pathNames.txt #"
echo "################################################"

echo "Main path to software is $MAINPATH"
echo "Path to Rscript is $RPATH"
echo "Path to Plink is $PLINKPATH"
echo "Path to Bagle is ${BEAGLEPATH}"

echo "#################################################"
echo "## GENOTYPE FILTERING THRESHOLDS               ##"
echo "#################################################"

echo "MAF threshold (to remove unimputable monomorphic loci): $MAF"
echo "MIND threshold (to remove samples with excess missing data): $MIND"
echo "GENO threshold (to remove loci with excess missing data): $GENO"


echo "#######################################"
echo "## STEP -1"
echo "## create unique folders for each run"
echo "#######################################"
tmstmp=$(date +%s)
currDate=$(date +%d-%m-%Y)
folderName=${PREFIX}_$( basename $INPUTFILE).${SAMPLESIZE}_$( basename $LDSIZE)_${tmstmp}.${currDate}

## CHECK IF OUTPUT FOLDER EXISTS
if [ -d "$OUTDIR" ]; then
    echo "$OUTDIR exists."
else
        echo "Output folder does not exist - TERMINATING";
        exit 1;
fi

echo "Folder name is: $folderName"
cd $OUTDIR

if [ -d "$folderName" ]; then
>&2 echo "!! ERROR: folder $folderName already exists !!"
exit
fi

mkdir $folderName
cd $folderName

echo "#######################################"
echo "## STEP 0"
echo "## sample individuals from the ped file"
echo "#######################################"
$PLINKPATH --$SPECIES --allow-extra-chr --bfile ${MAINPATH}/${INPUTFILE} --recode transpose --out transposed
$RPATH --vanilla ${MAINPATH}/heterogeneousImputation/scripts/sampleRows.R ${INPUTFILE}.ped $SAMPLESIZE $MAINPATH
$PLINKPATH --$SPECIES --allow-extra-chr --bfile ${MAINPATH}/${INPUTFILE} --keep keepIDs.txt --maf $MAF --bp-space 1 --recode --out subset
$PLINKPATH --$SPECIES --allow-extra-chr --file subset --freq --out subset

echo "#######################################"
echo "## STEP 1"
echo "## sample individuals from the ped file"
echo "#######################################"
## STEP 1
## Prepare file with HD and LD (missing values) SNP genotypes
# randomly sample individuals to assign to the LD array
if [[ "$OSTYPE" == "darwin"* ]]; then
  cut -f1-2 -d' ' subset.ped | gshuf -n $LDSIZE > keep.ids
else
  cut -f1-2 -d' ' subset.ped | shuf -n $LDSIZE > keep.ids
fi

# create low-density and high-density subsets
$PLINKPATH --$SPECIES --allow-extra-chr --file subset --keep keep.ids --extract ${LOWDENSITY} --recode --out subsetLD
$PLINKPATH --$SPECIES --allow-extra-chr --file subset --remove keep.ids --recode --out subsetHD
# put HD and LD subsets together into a combined file
$PLINKPATH --$SPECIES --allow-extra-chr --file subsetHD --merge subsetLD --allow-no-sex --recode --out combined
rm subsetHD* subsetLD.ped

echo "#######################################"
echo "## STEP 1.5"
echo "## recode the ped file into a .raw file"
echo "#######################################"
$PLINKPATH --$SPECIES --allow-extra-chr --file subset --recode A --out originalRaw
$PLINKPATH --$SPECIES --allow-extra-chr --file combined --recode A --out combinedRaw
rm combinedRaw.nosex combinedRaw.log combined.bed combined.bim originalRaw.nosex originalRaw.log

echo "#######################################"
echo "## STEP 2"
echo "## imputation of missing genotypes"
echo "#######################################"
## STEP 2
## Imputation of missing genotypes
date +%s > anfangZeit
## Imputation of missing genotypes
$PLINKPATH --$SPECIES --allow-extra-chr --file combined --recode vcf --out combined

if [ $use_singularity = 1 ]; then
	echo "beagle run from the singularity container"
	singularity run $BEAGLESING beagle gt=combined.vcf out=imputed
else
	echo "beagle run from the native java runtime environment"
	java -Xss5m -Xmx4g -jar $BEAGLEPATH gt=combined.vcf out=imputed
fi

$PLINKPATH --$SPECIES --allow-extra-chr --vcf imputed.vcf.gz --recode --out imputed
$PLINKPATH --$SPECIES --allow-extra-chr --file imputed --recode A --out imputed
rm *.vcf *.gz

echo "#######################################"
echo "## STEP 3"
echo "## Caclulate MAF"
echo "#######################################"
## STEP 3
## MAF calculation
$PLINKPATH --$SPECIES --allow-extra-chr --file imputed --freq --out freq 
rm freq.log freq.nosex

echo "#######################################"
echo "## STEP 4"
echo "## parsing results"
echo "#######################################"
## STEP 4
## parsing results
$RPATH --vanilla ${MAINPATH}/heterogeneousImputation/scripts/parseResults_density.R originalRaw.raw combinedRaw.raw imputed.raw ${LDSIZE} $( basename $INPUTFILE) ${MAINPATH}
rm originalRaw.raw imputed.raw imputed.map imputed.ped combinedRaw.raw freq.frq combined.* subset.* *.hh subsetLD.* transposed.*