https://wiki.anunna.wur.nl/index.php?action=history&feed=atom&title=Extract_noncall_snps_from_soyExtract noncall snps from soy - Revision history2026-04-18T01:49:33ZRevision history for this page on the wikiMediaWiki 1.43.1https://wiki.anunna.wur.nl/index.php?title=Extract_noncall_snps_from_soy&diff=1293&oldid=prevMegen002: Created page with "== Extract the SNPs that are not called in Soy == From the <code>.ped</code> file, only the soy samples were extracted, and call rate calculated like this: <source lang='bas..."2014-03-26T12:28:43Z<p>Created page with "== Extract the SNPs that are not called in Soy == From the <code>.ped</code> file, only the soy samples were extracted, and call rate calculated like this: <source lang='bas..."</p>
<p><b>New page</b></p><div>== Extract the SNPs that are not called in Soy ==<br />
<br />
From the <code>.ped</code> file, only the soy samples were extracted, and call rate calculated like this:<br />
<br />
<source lang='bash'><br />
./plink --file soy --out soy --missing-genotype N --freq --noweb<br />
</source><br />
<br />
The output data looks like this:<br />
CHR SNP A1 A2 MAF NCHROBS<br />
1 MARC0044150 0 0 NA 0<br />
1 ASGA0000014 A C 0.5 26<br />
1 ASGA0000021 G T 0.375 16<br />
1 ALGA0000009 G A 0.1333 30<br />
1 ALGA0000014 T C 0.5 28<br />
1 H3GA0000032 A C 0.5 38<br />
1 ASGA0000005 C A 0.4 20<br />
1 INRA0000015 A G 0.4643 28<br />
1 M1GA0000060 T C 0.06667 30<br />
<br />
Subsequently, only SNPs with less than 4 calls (i.e. 2 individuals) were extracted. This amounts to only a single Soy sample being called for that SNP.<br />
<source lang='bash'><br />
cat soy.frq | sed 's/ \+/\t/g' | sort -k6 -nr | awk '$6<4' | cut -f3 >noncallsnps2.txt<br />
</source><br />
Note that the 3rd column is selected because the first column is shifted due to leading whitespace.<br />
<br />
== Make reduced SNP set in Plink ==<br />
From the original SNP set (bones2.ped), that contains all samples (including the soy samples) for 49802 SNPs that map to autosomes in Ssc10.2, a reduced SNP set is made that only contains a subset of 8552 SNPs for the same individuals (including soy).<br />
<source lang='bash'><br />
./plink --file bones2 --out bones2_reduced --missing-genotype N --noweb --extract noncallsnps2.txt --recode<br />
</source><br />
<br />
== Calculate heterozygosities from <code>.ped</code> file ==<br />
<br />
Python3 script to calculate heterozygosity from <code>.ped</code> file:<br />
<source lang='python'><br />
import argparse<br />
import sys<br />
import os<br />
<br />
parser = argparse.ArgumentParser( description='calculate heterozygosities for individuals in a ped file')<br />
parser.add_argument("-p", "--pedfile", help="name of pedfile", nargs=1)<br />
<br />
args = parser.parse_args()<br />
ped=args.pedfile[0]<br />
<br />
def pedfile_to_list_of_lists(pedfile):<br />
flist=[]<br />
with open(pedfile) as pedf:<br />
for l in pedf.readlines():<br />
l=l.rstrip().split()<br />
flist.append(l[:])<br />
return flist<br />
<br />
def calc_het(geno):<br />
het=0<br />
hom=0<br />
for i in range(0,len(geno),2):<br />
pair=geno[i:i+2]<br />
if pair[0] == pair[1] and pair[0] != 'N':<br />
hom+=1<br />
elif pair[0] != pair[1] and pair[0] != 'N':<br />
het+=1<br />
obshet=het/(het+hom)<br />
return [obshet,het,hom]<br />
<br />
if __name__=="__main__":<br />
allgenos=pedfile_to_list_of_lists(ped)<br />
for genos in allgenos:<br />
het=calc_het(genos[6:])<br />
print(genos[0],genos[1],het[0],het[1],het[2],sep='\t')<br />
<br />
</source><br />
<br />
The script is invoked like this:<br />
<br />
<source lang='bash'><br />
python3 het_stats.py -p bones2_reduced.ped<br />
</source><br />
<br />
The output will look something like this:<br />
BABA HAPmap3_box3_G9_BABA_29 0.003549362 27 7580<br />
BABA HAPmap3_box3_H9_BABA_31 0.007805386 60 7627<br />
BK22 WG0093132-DNAH08_SSBK22M01 0.247806248 2118 6429<br />
BK23 WG0093131-DNAA02_SSBK23F09 0.243493552 2077 6453<br />
DU22 WG0093131-DNAC01_SSDU22F01 0.287118287 2454 6093<br />
DU23 WG0093132-DNAD11_SSDU23M06 0.290866511 2484 6056<br />
GL100 HAPmap5_Box7_B03_SSRB417_U01 0.137879911 186 1163<br />
GL101 HAPmap5_Box7_C03_SSRB420_U01 0.010806217 89 8147<br />
GL102 HAPmap5_Box7_D03_SSRB421_U01 0.014296464 114 7860<br />
<br />
== See also ==<br />
[[Bioinformatics_protocols_ABG_Chairgroup | Bioinformatics tips, tricks, workflows at ABGC]]<br />
<br />
== External links ==<br />
[http://pngu.mgh.harvard.edu/~purcell/plink/index.shtml Plink website and manual]</div>Megen002