ABCA6.sh

#!/bin/bash
#$ -S /bin/bash
#$ -cwd

# Check if CNVs are LD with SNPs.

for file in Samples*txt
do
	gene=`echo $file|sed 's/.txt//'|sed 's/Samples//'`
	# Extract chr and (one) position from file
	chr=`tail -n+2 Samples${gene}.txt | cut -f2 | head -n1|cut -d'r' -f2`
	pos=`tail -n+2 Samples${gene}.txt | cut -f4 | head -n1`
	echo "chr=$chr pos=$pos gene=$gene"

	# Extract relevant part of chromosome as plink file
	[ -f chr${chr}.map ] || plink19 --bfile /eva/brutus/J8/data/deepExom/burden5/clean/lrt24qibin7R2BeagleV2 --chr $chr --recode --out chr$chr

	# affected lists particids  with CNV
	tail -n+2 Samples${gene}.txt | cut -f17 > affected

	# Make a genotype column by going throug each particid in plink file 
	# and  checking if it is among affected
	rm -f chr${chr}_${gene}_column.txt
	for i in `cat chr${chr}.ped|cut -d' ' -f1`
	do
		echo $i
		if [ `grep -c -w "$i" affected` -eq 0 ]; then
			echo  "N N" >> chr${chr}_${gene}_column.txt
		else
			echo "Y N" >> chr${chr}_${gene}_column.txt
		fi
	done

	# Create new set of plink files with CNV column added 
	paste -d' ' chr${chr}.ped chr${chr}_${gene}_column.txt > chr${chr}_${gene}.ped
	cp chr${chr}.map chr${chr}_${gene}.map
	echo -en "${chr}\t${gene}\t0\t${pos}\n" >> chr${chr}_${gene}.map

	# Check LD in new plink file	
	plink19 --file chr${chr}_${gene} --r2 --ld-snp ${gene} --ld-window 99999 --ld-window-kb 10000  --ld-window-r2 0 --out ${gene}


cat << RscriptEOF > plot.$gene.R
library(ggplot2)
ld <- read.table("$gene.ld", header=T)
pdf("$gene.pdf")
ggplot(ld)+geom_point(aes(x=BP_B - BP_A, y = R2, col=(SNP_B == "$gene")))+guides(colour=FALSE)+ggtitle("$gene LD plot") + xlab("Distance to CNV")
dev.off()
RscriptEOF
	Rscript plot.$gene.R
done