Changes between Version 8 and Version 9 of GwasQcPipeline
- Timestamp:
- Feb 10, 2011 5:35:13 PM (14 years ago)
Legend:
- Unmodified
- Added
- Removed
- Modified
-
GwasQcPipeline
v8 v9 36 36 Get missing data information: 37 37 38 38 {{{p-link --bfile GvNL_raw_final --missing --out GvNL_raw_final}}} 39 39 40 40 Get heterozygocity information: … … 45 45 * In R: 46 46 47 {{{het <- read.table("GvNL_raw_final.het", head=TRUE) mis <- read.table("GvNL_raw_final.imiss", head=TRUE) het$HET_RATE = (het$"N.NM." - het$"O.HOM.")/het$"N.NM." plot (mis$F_MISS,het$HET_RATE, xlab="Missing SNPs", ylab="Heterozygocity Rate")}}} 47 {{{ 48 het <- read.table("GvNL_raw_final.het", head=TRUE) 49 mis <- read.table("GvNL_raw_final.imiss", head=TRUE); 50 het$HET_RATE = (het$"N.NM." - het$"O.HOM.")/het$"N.NM." plot (mis$F_MISS,het$HET_RATE, xlab="Missing SNPs", ylab="Heterozygocity Rate")}}} 48 51 49 52 No samples with callrate < 95% (Filtered in lab) … … 53 56 In R: 54 57 55 {{{het_fail = subset(het, (het$HET_RATE < mean(het$HET_RATE)-3*sd(het$HET_RATE)) | (het$HET_RATE > mean(het$HET_RATE)+3*sd(het$HET_RATE))) het_fail$HET_DST = (het_fail$HET_RATE-mean(het$HET_RATE))/sd(het$HET_RATE) write.table(het_fail, "fail-het-qc.txt", row.names=FALSE)}}} 58 {{{ 59 het_fail = subset(het, (het$HET_RATE < mean(het$HET_RATE)-3*sd(het$HET_RATE)) | (het$HET_RATE > mean(het$HET_RATE)+3*sd(het$HET_RATE))); 60 het_fail$HET_DST = (het_fail$HET_RATE-mean(het$HET_RATE))/sd(het$HET_RATE); 61 write.table(het_fail, "fail-het-qc.txt", row.names=FALSE)}}} 56 62 57 63 Result: