Skip to content

Instantly share code, notes, and snippets.

@nievergeltlab
Created March 22, 2017 22:08
Show Gist options
  • Save nievergeltlab/892ea41687e40a67bff1804b85e5ed78 to your computer and use it in GitHub Desktop.
Save nievergeltlab/892ea41687e40a67bff1804b85e5ed78 to your computer and use it in GitHub Desktop.
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment