nievergeltlab
3/22/2017 - 10:08 PM

2 file IBD, if strand mismatches and allele mismatches

2 file IBD, if strand mismatches and allele mismatches

bimfile=pts_meg2_mix_am-qc.bim
grep -P "A\tT" $bimfile > ambiguous_snps.txt
grep -P "T\tA" $bimfile >> ambiguous_snps.txt
grep -P "C\tG" $bimfile >> ambiguous_snps.txt
grep -P "G\tC" $bimfile >> ambiguous_snps.txt


#filer out ambiguous
$plink_location --bfile pts_meg2_mix_am-qc --exclude ambiguous_snps.txt --make-bed --out pts_meg2_mix_am-qc-ambig

$plink_location --bfile pts_gtpc_mix_am-qc --bmerge pts_meg2_mix_am-qc-ambig.bed pts_meg2_mix_am-qc-ambig.bim pts_meg2_mix_am-qc-ambig.fam --make-bed --out gtp_meg2

$plink_location --bfile pts_meg2_mix_am-qc-ambig --flip gtp_meg2-merge.missnp --make-bed --out pts_meg2_mix_am-qc-ambig-flip

$plink_location --bfile pts_gtpc_mix_am-qc --bmerge pts_meg2_mix_am-qc-ambig-flip.bed pts_meg2_mix_am-qc-ambig-flip.bim pts_meg2_mix_am-qc-ambig-flip.fam --make-bed --out gtp_meg2-flip

#and multi allelic
 

$plink_location --bfile pts_meg2_mix_am-qc-ambig-flip --exclude gtp_meg2-flip-merge.missnp --make-bed --out pts_meg2_mix_am-qc-ambig-flip-diallele

$plink_location --bfile pts_gtpc_mix_am-qc --bmerge pts_meg2_mix_am-qc-ambig-flip-diallele.bed pts_meg2_mix_am-qc-ambig-flip-diallele.bim pts_meg2_mix_am-qc-ambig-flip-diallele.fam --allow-no-sex --make-bed --out gtp_meg2-flip-diallele

#Filter to just the overlap
LC_ALL=C join <(awk '{print $2}' pts_meg2_mix_am-qc.bim | LC_ALL=C sort -k1b,1) <(awk '{print $2}' pts_gtpc_mix_am-qc.bim | LC_ALL=C sort -k1b,1) > overlap.snps

$plink_location --bfile gtp_meg2-flip-diallele --extract overlap.snps --maf 0.05 --allow-no-sex --make-bed --out gtp_meg2_use



mkdir temporary_files
$plink_location --bfile gtp_meg2_use --maf 0.05 --geno 0.02 --mind 0.15 --indep-pairwise 50 5 0.2 --out temporary_files/gtp_meg2_use_ibd


 $plink_location --bfile gtp_meg2_use --mind 0.15 --allow-no-sex --extract  temporary_files/gtp_meg2_use_ibd.prune.in --genome --min 0.5 --out  temporary_files/gtp_meg2_use_ibd

 awk '{if(NR ==1) print "FID","IID"; else print $3,$4}' temporary_files/gtp_meg2_use_ibd.genome > temporary_files/gtp_meg2_use_ibd.remove