|
一般自然群體，基因型個體的雜合度過高或者過低，都不正常，我們需要根據(jù)雜合度進(jìn)行過濾。偏差可能表明樣品受到污染，近親繁殖。我們建議刪除樣品雜合率平均值中偏離±3 SD的個體。 ?我的理解：非自然群體中，比如自交系，雜交種F1，這些群體不需要過濾雜合度。 ? 「參數(shù)過濾和手動過濾」plink有個特點(diǎn)，所有的過濾標(biāo)準(zhǔn)，都可以生成過濾前的文件，然后可以手動過濾，也可以用參數(shù)進(jìn)行過濾。 - 比如：
--missing生成結(jié)果，也可以用--geno和--mind過濾。 - 比如：
--hardy生成結(jié)果，可以使用--hwe過濾 - 比如：
--freq生成結(jié)果，可以用--maf過濾
但是雜合度--het，沒有過濾的函數(shù)，只能通過編程去提取ID，然后用--remove去實(shí)現(xiàn)。
1. 計(jì)算雜合度plink --bfile HapMap_3_r3_9 --het --out R_check
結(jié)果文件： .het (method-of-moments F coefficient estimates)
Produced by --het.
A text file with a header line, and one line per sample with the following six fields:
FIDFamily ID # 家系ID
IIDWithin-family ID # 個體ID
O(HOM)Observed number of homozygotes # 實(shí)際純合個數(shù)
E(HOM)Expected number of homozygotes # 期望純合個數(shù)
N(NM)Number of non-missing autosomal genotypes # 總個數(shù)
FMethod-of-moments F coefficient estimate #
所以，雜合度 = (N-O)/N 2. 雜合度可視化R代碼： het <- read.table("R_check.het", head=TRUE)
pdf("heterozygosity.pdf")
het$HET_RATE = (het$"N.NM." - het$"O.HOM.")/het$"N.NM."
hist(het$HET_RATE, xlab="Heterozygosity Rate", ylab="Frequency", main= "Heterozygosity Rate")
dev.off()
可視化： 3. 計(jì)算雜合度三倍標(biāo)準(zhǔn)差以外的個體首先，查看哪些個體在3倍標(biāo)準(zhǔn)差以外： het <- read.table("R_check.het", head=TRUE)
het$HET_RATE = (het$"N.NM." - het$"O.HOM.")/het$"N.NM."
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)
結(jié)果可以看出，這兩個個體雜合度在3倍標(biāo)準(zhǔn)差以外： 4. 去掉這些個體cat fail-het-qc.txt
"FID" "IID" "O.HOM." "E.HOM." "N.NM." "F" "HET_RATE" "HET_DST"
1330 "NA12342" 703770 691000 1066256 0.03402 0.33996151018142 -4.75272486494316
1459 "NA12874" 709454 695300 1072858 0.03758 0.338725162137021 -5.18288854902555
先對數(shù)據(jù)進(jìn)行清洗，去掉引號，然后提取家系和個體ID sed 's/"http://g' fail-het-qc.txt |awk '{print $2}' > het_fail_ind.txt
sed 's/"http://g' fail-het-qc.txt |awk '{print $1,$2}' > het_fail_ind.txt
使用remove去掉這兩個個體 plink --bfile HapMap_3_r3_9 --remove het_fail_ind.txt --make-bed --out HapMap_3_r3_10
5. 結(jié)果文件HapMap_3_r3_10.bed HapMap_3_r3_10.bim HapMap_3_r3_10.fam HapMap_3_r3_10.log
GWAS系列： 筆記 | GWAS 操作流程1：下載數(shù)據(jù)
筆記 | GWAS 操作流程2-1：缺失質(zhì)控
筆記 | GWAS 操作流程2-2：性別質(zhì)控
筆記GWAS 操作流程2-3：MAF過濾
筆記 | GWAS 操作流程2-4：哈溫平衡檢驗(yàn)
PS吐槽：
微信編輯器中的文章引用功能有bug，上面的超鏈接有亂碼&；|，這都是什么鬼，我不想手動刪除了，let it be！ 我的公眾號：
|
