Created
March 22, 2017 22:08
-
-
Save nievergeltlab/892ea41687e40a67bff1804b85e5ed78 to your computer and use it in GitHub Desktop.
2 file IBD, if strand mismatches and allele mismatches
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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