Supercomuter: 3.84 TB RAM 512 CPS ############### LIKELIHOOD: G-H/G G-D/G G-L/G H-D/H Min. :0.4340 Min. :0.6057 Min. :0.4686 Min. :0.5838 1st Qu.:0.6021 1st Qu.:0.7083 1st Qu.:0.5761 1st Qu.:0.6744 Median :0.6362 Median :0.7329 Median :0.6069 Median :0.7000 Mean :0.6357 Mean :0.7348 Mean :0.6063 Mean :0.7018 3rd Qu.:0.6697 3rd Qu.:0.7610 3rd Qu.:0.6355 3rd Qu.:0.7285 Max. :0.7864 Max. :0.8643 Max. :0.7358 Max. :0.8198 H-L/H D-L/D G-H/H G-D/D Min. :0.4943 Min. :0.4857 Min. :0.5925 Min. :0.6211 1st Qu.:0.5894 1st Qu.:0.5819 1st Qu.:0.6844 1st Qu.:0.7150 Median :0.6166 Median :0.6115 Median :0.7121 Median :0.7424 Mean :0.6169 Mean :0.6113 Mean :0.7116 Mean :0.7424 3rd Qu.:0.6453 3rd Qu.:0.6374 3rd Qu.:0.7397 3rd Qu.:0.7668 Max. :0.7733 Max. :0.7358 Max. :0.8391 Max. :0.8622 G-L/L H-D/D H-L/L D-L/L Min. :0.5815 Min. :0.4519 Min. :0.4737 Min. :0.5475 1st Qu.:0.6768 1st Qu.:0.6008 1st Qu.:0.6094 1st Qu.:0.6745 Median :0.7067 Median :0.6318 Median :0.6432 Median :0.7052 Mean :0.7081 Mean :0.6335 Mean :0.6439 Mean :0.7068 3rd Qu.:0.7373 3rd Qu.:0.6650 3rd Qu.:0.6782 3rd Qu.:0.7374 Max. :0.8452 Max. :0.8130 Max. :0.8378 Max. :0.8650 > diag(var(LD)) G-H/G G-D/G G-L/G H-D/H H-L/H D-L/D 0.002364580 0.001545414 0.001866014 0.001644833 0.001789776 0.001776777 G-H/H G-D/D G-L/L H-D/D H-L/L D-L/L 0.001682689 0.001543011 0.002016111 0.002430299 0.002896852 0.002068685 0|0:0.000:-0.19,-0.45,-2.24 0|1:0.900:-0.19,-0.46,-2.51 0|0:0.000:-0.20,-0.44,-1.89 1|0:0.900:-0.21,-0.43,-1.85 1|1:1.400:0,0,0 1|0:1.200:0,0,0 0|0:0.000:-0.48,-0.48,-0.48 1|0:1.000:-0.48,-0.48,-0.48 10^(-0.48,-0.48,-0.48)=(0.33,0.33,0.33) 1|0:1.000:-0.41,-0.44,-0.62 10^(-0.41,-0.44,-0.62)=(0.39,0.36,0.24) 1|0:0.950:-0.10,-0.68,-5.00 1|0:0.950:-0.10,-0.67,-3.38 10^(-0.10,-0.67,-3.38)=(0.79,0.21,0.00) Type <- c("G","H","D","L") TypeL <- c("Genotype","Haplotype","Dosage","Likelihood") [1] "" [1] "" [1] "Genotype-Haplotype_1040000_1050000" [1] "#IBD Genotype:" [1] 177 [1] "#IBD Haplotype:" [1] 182 [1] "#IBD Genotype-Haplotype:" [1] 181 [1] "#IBD PairsGenotype-Haplotype:" [1] 131 [1] "" [1] "" [1] "Genotype-Dosage_1040000_1050000" [1] "#IBD Genotype:" [1] 177 [1] "#IBD Dosage:" [1] 190 [1] "#IBD Genotype-Dosage:" [1] 185 [1] "#IBD PairsGenotype-Dosage:" [1] 137 [1] "" [1] "" [1] "Genotype-Likelihood_1040000_1050000" [1] "#IBD Genotype:" [1] 177 [1] "#IBD Likelihood:" [1] 181 [1] "#IBD Genotype-Likelihood:" [1] 195 [1] "#IBD PairsGenotype-Likelihood:" [1] 125 [1] "" [1] "" [1] "Haplotype-Dosage_1040000_1050000" [1] "#IBD Haplotype:" [1] 182 [1] "#IBD Dosage:" [1] 190 [1] "#IBD Haplotype-Dosage:" [1] 193 [1] "#IBD PairsHaplotype-Dosage:" [1] 130 [1] "" [1] "" [1] "Haplotype-Likelihood_1040000_1050000" [1] "#IBD Haplotype:" [1] 182 [1] "#IBD Likelihood:" [1] 181 [1] "#IBD Haplotype-Likelihood:" [1] 207 [1] "#IBD PairsHaplotype-Likelihood:" [1] 112 [1] "" [1] "" [1] "Dosage-Likelihood_1040000_1050000" [1] "#IBD Dosage:" [1] 190 [1] "#IBD Likelihood:" [1] 181 [1] "#IBD Dosage-Likelihood:" [1] 207 [1] "#IBD PairsDosage-Likelihood:" [1] 128 [1] "" [1] "" [1] "Genotype-Haplotype_1000000_1010000" [1] "#IBD Genotype:" [1] 205 [1] "#IBD Haplotype:" [1] 206 [1] "#IBD Genotype-Haplotype:" [1] 218 [1] "#IBD Pairs Genotype-Haplotype:" [1] 143 [1] "" [1] "" [1] "Genotype-Dosage_1000000_1010000" [1] "#IBD Genotype:" [1] 205 [1] "#IBD Dosage:" [1] 207 [1] "#IBD Genotype-Dosage:" [1] 205 [1] "#IBD Pairs Genotype-Dosage:" [1] 155 [1] "" [1] "" [1] "Genotype-Likelihood_1000000_1010000" [1] "#IBD Genotype:" [1] 205 [1] "#IBD Likelihood:" [1] 210 [1] "#IBD Genotype-Likelihood:" [1] 226 [1] "#IBD Pairs Genotype-Likelihood:" [1] 144 [1] "" [1] "" [1] "Haplotype-Dosage_1000000_1010000" [1] "#IBD Haplotype:" [1] 206 [1] "#IBD Dosage:" [1] 207 [1] "#IBD Haplotype-Dosage:" [1] 222 [1] "#IBD Pairs Haplotype-Dosage:" [1] 137 [1] "" [1] "" [1] "Haplotype-Likelihood_1000000_1010000" [1] "#IBD Haplotype:" [1] 206 [1] "#IBD Likelihood:" [1] 210 [1] "#IBD Haplotype-Likelihood:" [1] 235 [1] "#IBD Pairs Haplotype-Likelihood:" [1] 136 [1] "" [1] "" [1] "Dosage-Likelihood_1000000_1010000" [1] "#IBD Dosage:" [1] 207 [1] "#IBD Likelihood:" [1] 210 [1] "#IBD Dosage-Likelihood:" [1] 238 [1] "#IBD Pairs Dosage-Likelihood:" [1] 136 [1] "" [1] "" [1] "Genotype-Haplotype_1500000_1510000" [1] "#IBD Genotype:" [1] 274 [1] "#IBD Haplotype:" [1] 249 [1] "#IBD Genotype-Haplotype:" [1] 285 [1] "#IBD Pairs Genotype-Haplotype:" [1] 176 [1] "" [1] "" [1] "Genotype-Dosage_1500000_1510000" [1] "#IBD Genotype:" [1] 274 [1] "#IBD Dosage:" [1] 250 [1] "#IBD Genotype-Dosage:" [1] 272 [1] "#IBD Pairs Genotype-Dosage:" [1] 180 [1] "" [1] "" [1] "Genotype-Likelihood_1500000_1510000" [1] "#IBD Genotype:" [1] 274 [1] "#IBD Likelihood:" [1] 227 [1] "#IBD Genotype-Likelihood:" [1] 288 [1] "#IBD Pairs Genotype-Likelihood:" [1] 167 [1] "" [1] "" [1] "Haplotype-Dosage_1500000_1510000" [1] "#IBD Haplotype:" [1] 249 [1] "#IBD Dosage:" [1] 250 [1] "#IBD Haplotype-Dosage:" [1] 274 [1] "#IBD Pairs Haplotype-Dosage:" [1] 167 [1] "" [1] "" [1] "Haplotype-Likelihood_1500000_1510000" [1] "#IBD Haplotype:" [1] 249 [1] "#IBD Likelihood:" [1] 227 [1] "#IBD Haplotype-Likelihood:" [1] 273 [1] "#IBD Pairs Haplotype-Likelihood:" [1] 153 [1] "" [1] "" [1] "Dosage-Likelihood_1500000_1510000" [1] "#IBD Dosage:" [1] 250 [1] "#IBD Likelihood:" [1] 227 [1] "#IBD Dosage-Likelihood:" [1] 272 [1] "#IBD Pairs Dosage-Likelihood:" [1] 158 [1] "" [1] "" [1] "Genotype-Haplotype_2000000_2010000" [1] "#IBD Genotype:" [1] 162 [1] "#IBD Haplotype:" [1] 152 [1] "#IBD Genotype-Haplotype:" [1] 164 [1] "#IBD Pairs Genotype-Haplotype:" [1] 112 [1] "" [1] "" [1] "Genotype-Dosage_2000000_2010000" [1] "#IBD Genotype:" [1] 162 [1] "#IBD Dosage:" [1] 178 [1] "#IBD Genotype-Dosage:" [1] 187 [1] "#IBD Pairs Genotype-Dosage:" [1] 117 [1] "" [1] "" [1] "Genotype-Likelihood_2000000_2010000" [1] "#IBD Genotype:" [1] 162 [1] "#IBD Likelihood:" [1] 161 [1] "#IBD Genotype-Likelihood:" [1] 185 [1] "#IBD Pairs Genotype-Likelihood:" [1] 106 [1] "" [1] "" [1] "Haplotype-Dosage_2000000_2010000" [1] "#IBD Haplotype:" [1] 152 [1] "#IBD Dosage:" [1] 178 [1] "#IBD Haplotype-Dosage:" [1] 182 [1] "#IBD Pairs Haplotype-Dosage:" [1] 109 [1] "" [1] "" [1] "Haplotype-Likelihood_2000000_2010000" [1] "#IBD Haplotype:" [1] 152 [1] "#IBD Likelihood:" [1] 161 [1] "#IBD Haplotype-Likelihood:" [1] 172 [1] "#IBD Pairs Haplotype-Likelihood:" [1] 103 [1] "" [1] "" [1] "Dosage-Likelihood_2000000_2010000" [1] "#IBD Dosage:" [1] 178 [1] "#IBD Likelihood:" [1] 161 [1] "#IBD Dosage-Likelihood:" [1] 193 [1] "#IBD Pairs Dosage-Likelihood:" [1] 109 [1] "" [1] "" [1] "Genotype-Haplotype_2500000_2510000" [1] "#IBD Genotype:" [1] 242 [1] "#IBD Haplotype:" [1] 224 [1] "#IBD Genotype-Haplotype:" [1] 270 [1] "#IBD Pairs Genotype-Haplotype:" [1] 146 [1] "" [1] "" [1] "Genotype-Dosage_2500000_2510000" [1] "#IBD Genotype:" [1] 242 [1] "#IBD Dosage:" [1] 236 [1] "#IBD Genotype-Dosage:" [1] 262 [1] "#IBD Pairs Genotype-Dosage:" [1] 169 [1] "" [1] "" [1] "Genotype-Likelihood_2500000_2510000" [1] "#IBD Genotype:" [1] 242 [1] "#IBD Likelihood:" [1] 198 [1] "#IBD Genotype-Likelihood:" [1] 265 [1] "#IBD Pairs Genotype-Likelihood:" [1] 137 [1] "" [1] "" [1] "Haplotype-Dosage_2500000_2510000" [1] "#IBD Haplotype:" [1] 224 [1] "#IBD Dosage:" [1] 236 [1] "#IBD Haplotype-Dosage:" [1] 275 [1] "#IBD Pairs Haplotype-Dosage:" [1] 142 [1] "" [1] "" [1] "Haplotype-Likelihood_2500000_2510000" [1] "#IBD Haplotype:" [1] 224 [1] "#IBD Likelihood:" [1] 198 [1] "#IBD Haplotype-Likelihood:" [1] 263 [1] "#IBD Pairs Haplotype-Likelihood:" [1] 123 [1] "" [1] "" [1] "Dosage-Likelihood_2500000_2510000" [1] "#IBD Dosage:" [1] 236 [1] "#IBD Likelihood:" [1] 198 [1] "#IBD Dosage-Likelihood:" [1] 274 [1] "#IBD Pairs Dosage-Likelihood:" [1] 134 For above 4 cases: Ratio Likelihood-Genotype/Genotype: 0.706 0.702 0.609 0.654 0.566 Thus take: 1000000_1010000 ############### ssh -X k337870@mach.edvz.uni-linz.ac.at scp compa*R k337870@mach.edvz.uni-linz.ac.at:/BIGtmp/scratch/k3307/k337870/SFSsimulation/scripts/ scp ./compareMethodsPlink*R k337870@mach.edvz.uni-linz.ac.at:/BIGtmp/scratch/k3307/k337870/SFSsimulation/scripts/ scp * k337870@mach.edvz.uni-linz.ac.at:/BIGtmp/scratch/k3307/k337870/SFSsimulation/scripts/ scp -P 5792 * hochreit@140.78.90.2:/seppdata/sepp/linkage/SFSsimulation/simAlong ################################### /system/apps/biosoft/R-2.14.0/bin/R source("compareMethodsD1.R") ################################### simAlong 1MB 6samples simBlong 1MB 2samples simClong 2MB 6samples screen -r 9958.pts-60.tiger seppdata simDlong 2MB 2samples 1 MB -> 140 - 250 minor alleles and 8000-9500 bases 2 MB -> 280 - 500 minor alleles and 17000-19000 bases hapFABIA: [1] "Start" [1] "Tue Mar 19 08:08:02 2013" [1] "End Fabia" [1] "Tue Mar 19 08:10:29 2013" [1] "End HapFabia" [1] "Tue Mar 19 08:11:04 2013" 3 min 2 s GERMLINE: [1] "Start GERMLINE" [1] "Wed Mar 20 12:08:51 2013" [1] "End GERMLINE" [1] "Wed Mar 20 12:09:27 2013" 36 s BEAGLE: [1] "Start BEAGLE" [1] "Wed Mar 20 12:15:24 2013" [1] "End BEAGLE" [1] "Wed Mar 20 22:44:10 2013" 10 h 28 m 46 s PLINK: [1] "Start PLINK" [1] "Wed Mar 20 16:03:10 2013" [1] "End PLINK" [1] "Sat Mar 23 11:17:14 2013" 2 days 19h 14m 4s 67h 14m 4s ##################### ##################### ##################### ##################### ##################### ##################### ##################### ##################### ### in each data... directory: sfs_code 1 1 -n 5000 -N 2000 -L 1 300000 -I -r 0.001 -W 2 0 0 0 0.184 0.00040244 -Td 0 1.98 -Td 0.264196 0.12859 -Tg 0.340566 44.75089 -Td 0.340566 0.554513 -Tg 0.39085 282.7192 -TE 0.404845 -o out.txt -e err.txt -F freq.txt convertSFS_CODE out.txt --structure F struc.txt doConvert.R ###(calls convert.R for each data... directory) combine.R ### combines the simulations ##################### ## DNA length: 44,094,874 ## SNVs: 1,148,822 ## individuals: 5,000 ## 85,951,308 minor alleles ## 440,948,740,000 bases ## avSNV in 1000 individuals: 418,000 ## DNA length: 44,094,874 ## av. SNV dist in 1000 individuals: 44,094,874/41,800 = 105.49 ## private # S1 <- which(colSums(GalleleI)==1) # length(S1) # [1] 689548 # 60% are private! 1000Genomes: 21.4% are private cd /seppdata/sepp/linkage/SFSsimulation/simE /system/apps/biosoft/R-2.14.0/bin/R source("read1000.R") /system/apps/biosoft/R-2.14.0/bin/R source("../compareMethods.R") simA noImplanted <- 1 ## how many haplotypes implanted <- 10 ## in how many samples lengthI <- 20 ## IBD length in kb simB noImplanted <- 1 ## how many haplotypes implanted <- 10 ## in how many samples lengthI <- 10 ## IBD length in kb simC noImplanted <- 1 ## how many haplotypes implanted <- 6 ## in how many samples lengthI <- 20 ## IBD length in kb simD noImplanted <- 20 ## how many haplotypes implanted <- 10 ## in how many samples lengthI <- 10 ## IBD length in kb simE noImplanted <- 5 ## how many haplotypes implanted <- 10 ## in how many samples lengthI <- 20 ## IBD length in kb PLINK! | v1.07 | 10/Aug/2009 Germline: germline -haploid -min_m 0.001 -err_hom 0 -err_het 0 -bits 128 -bits 128 260 and 390 as filter! too many hits (30 Mio) DASH and Fabia: minsamples 5 (Dash and Fabia) Plink: pruned because of all SNVs cd ./simA /system/apps/biosoft/R-2.14.0/bin/R source("../testSignificance.R") q() n cd .. LONG: ##### BEAGLE RUN: java -Djava.io.tmpdir=/tmp -jar beagle.jar unphased=databeagle fastibd=true missing=N out=beagleres.txt Filter A-D: res1: 1.0e-50 ; res2: 1.0e-80 E F: res1: 1.0e-40 ; res2: 1.0e-60 PLINK RUN: Filter A-D: plink --file dataplink --indep 50 5 2 --out dataplink plink --file dataplink --extract dataplink.prune.in --make-bed --out dataplinkpruned plink --bfile dataplinkpruned --noweb --genome --allow-no-sex --out plink plink --bfile dataplinkpruned --noweb --read-genome plink.genome --segment --all-pairs --segment-length 600 --segment-snp 6000 --out plink Filter E und F: plink --bfile dataplinkpruned --noweb --read-genome plink.genome --segment --all-pairs --segment-length 300 --segment-snp 3000 --out plink GERMLINE RUN: file "germline_input": 1 dataplink.map dataplink.ped germline_output Filter A-D: germline -haploid -min_m 0.5 -err_hom 0 -err_het 0 -bits 1000 < germline_input res1: ## Filter 7000 res2: ## Filter 5000 Filter E und F: germline -haploid -min_m 0.3 -err_hom 0 -err_het 0 -bits 700 < germline_input res1: ## Filter 4000 res2: ## Filter 3000 FABIA: p=10 #number bicluster per iteration cyc=50 #number cycles alpha=0.03 #sparseness parameter iter=1 # variable HapFabia: thresCount <- 1e-5 # p-value of random histogram hit, default 1e-5 minSNPsFactor <- 3/4 # percentage of IBD overlap; 1/2 for large to 3/4 for for sps <- 0.9 #top L used for identifying accumulations psZ <- 0.8 #top Z to identify chromosomes/individuals I <- 200000 # variable # number SNVs per interval IBDlength <- 5000 # variable # Length of IBD interval of interest in kb pMAF <- 0.035 # variable # average minor allele frequency Filter E und F: IBDlength <- 1000 SHORT: ##### BEAGLE RUN: java -Djava.io.tmpdir=/tmp -jar beagle.jar unphased=databeagle fastibd=true missing=N out=beagleres.txt Filter res1: 1.0e-10 ; res2: 1.0e-13 --------------- Shuffeled data: filter: > 5kb --------------- PLINK RUN: plink --file dataplink --indep 50 5 2 --out dataplink plink --file dataplink --extract dataplink.prune.in --make-bed --out dataplinkpruned plink --bfile dataplinkpruned --noweb --genome --allow-no-sex --out plink plink --bfile dataplinkpruned --noweb --read-genome plink.genome --segment --all-pairs --segment-length 5 --segment-snp 50 --out plink --------------- Simulated data had better results with (less private SNVs): plink --bfile dataplinkpruned --noweb --read-genome plink.genome --segment --all-pairs --segment-length 0.1 --segment-snp 20 --out plink ---------------- --------------- Shuffeled data: filter: > 5kb plink --bfile dataplinkpruned --noweb --read-genome plink.genome --segment --all-pairs --segment-length 0.1 --segment-snp 20 --out plink --------------- GERMLINE RUN: file "germline_input": 1 dataplink.map dataplink.ped germline_output germline -haploid -min_m 0.001 -err_hom 0 -err_het 0 -bits 170 < germline_input ## Filter 180 ## Filter 200 --------------- Simulated data gave better results with (less private SNVs): germline -haploid -min_m 0.001 -err_hom 0 -err_het 0 -bits 30 < germline_input ## Filter 150 ## Filter 200 For up to 6 mismatches AMis, BMis, CMis, DMis germline -haploid -min_m 0.001 -err_hom 6 -err_het 6 -bits 30 < germline_input ------------- --------------- Shuffeled data above setting gave 100Mio hits and FDR approx 1 Therefore: bits 128 and stringenter Filters filter: > 5kb germline -haploid -min_m 0.001 -err_hom 0 -err_het 0 -bits 128 < germline_input ## Filter 260 ## Filter 390 --------------- FABIA: p=10 #number bicluster per iteration cyc=50 #number cycles alpha=0.03 #sparseness parameter iter=1 # variable HapFabia: thresCount <- 1e-5 # p-value of random histogram hit, default 1e-5 minSNPsFactor <- 3/4 # percentage of IBD overlap; 1/2 for large to 3/4 for for sps <- 0.9 #top L used for identifying accumulations psZ <- 0.8 #top Z to identify chromosomes/individuals I <- 10000 # variable # number SNVs per interval IBDlength <- 50 # variable # Length of IBD interval of interest in kb pMAF <- 0.035 # variable # average minor allele frequency --------------- Simulated data had MAF of 0.04, therefore: pMAF <- 0.04 iter=4 # for 10 or 20 IBD segments --------------- --------------- Shuffeled data: filter: > 5kb and more than 4 chromosomes iter=2 for all settings --------------- DASH: cut -f 1,2,4 germline_result_filtered | dash_cc dataplink.fam dash_output -min 2 --------------- Shuffeled data: cut -f 1,2,4 germline_result_filtered | dash_cc dataplink.fam dash_output -min 5 filter: > 5kb --------------- 1000Genomes: FABIA: p=10 #number bicluster per iteration cyc=50 #number cycles alpha=0.03 #sparseness parameter iter=40 # variable HapFabia: thresCount <- 1e-5 # p-value of random histogram hit, default 1e-5 minSNPsFactor <- 3/4 # percentage of IBD overlap; 1/2 for large to 3/4 for for sps <- 0.9 #top L used for identifying accumulations psZ <- 0.8 #top Z to identify chromosomes/individuals I <- 10000 # variable # number SNVs per interval IBDlength <- 50 # variable # Length of IBD interval of interest in kb pMAF <- 0.035 # variable # average minor allele frequency