|
上一次我們經(jīng)過去掉缺失，去掉錯誤的性別信息，得到的文件為： HapMap_3_r3_6.bed HapMap_3_r3_6.fam HapMap_3_r3_6.log
HapMap_3_r3_6.bim HapMap_3_r3_6.hh
這里，我們根據(jù)最小等位基因頻率（MAF）去篩選。 「為什么要根據(jù)MAF去篩選？」 ?最小等位基因頻率怎么計算？比如一個位點有AA或者AT或者TT，那么就可以計算A的基因頻率和T的基因頻率，qA + qT = 1，這里誰比較小，誰就是最小等位基因頻率，比如qA = 0.3, qT = 0.7， 那么這個位點的MAF為0.3. 之所以用這個過濾標(biāo)準(zhǔn)，是因為MAF如果非常小，比如低于0.02，那么意味著大部分位點都是相同的基因型，這些位點貢獻的信息非常少，增加假陽性。更有甚者MAF為0，那就是所有位點只有一種基因型，這些位點沒有貢獻信息，放在計算中增加計算量，沒有意義，所以要根據(jù)MAF進行過濾。 ? 1. 去掉性染色體上的位點「思路：」 「提取常染色體上的位點名稱：」 因為這里是人的數(shù)據(jù)，所以染色體只需要去1~22的常染色體，提取它的家系ID和個體ID，后面用于提取。 awk '{ if($1 >=1 && $1 <= 22) print $2}' HapMap_3_r3_6.bim > snp_1_22.txt
wc -l snp_1_22.txt
1399306 snp_1_22.txt
常染色體上共有1399306位點。 2. 提取常染色體上的位點這里，用到了位點提取參數(shù)--extract plink --bfile HapMap_3_r3_6 --extract snp_1_22.txt --make-bed --out HapMap_3_r3_7
「查看日志：」 PLINK v1.90b6.5 64-bit (13 Sep 2018) www.cog-genomics.org/plink/1.9/
(C) 2005-2018 Shaun Purcell, Christopher Chang GNU General Public License v3
Logging to HapMap_3_r3_7.log.
Options in effect:
--bfile HapMap_3_r3_6
--extract snp_1_22.txt
--make-bed
--out HapMap_3_r3_7
515185 MB RAM detected; reserving 257592 MB for main workspace.
1431211 variants loaded from .bim file.
163 people (79 males, 84 females) loaded from .fam.
112 phenotype values loaded from .fam.
--extract: 1399306 variants remaining.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 112 founders and 51 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.998153.
1399306 variants and 163 people pass filters and QC.
Among remaining phenotypes, 56 are cases and 56 are controls. (51 phenotypes
are missing.)
--make-bed to HapMap_3_r3_7.bed + HapMap_3_r3_7.bim + HapMap_3_r3_7.fam ...
done.
可以看到，共有163個基因型，共有1399306個SNP 3. 計算每個SNP位點的基因頻率首先，通過參數(shù)--freq，計算每個SNP的MAF頻率，通過直方圖查看整體分布。可視化會更加直接。 plink --bfile HapMap_3_r3_7 --freq --out MAF_check
結(jié)果文件：MAF_check.frq預(yù)覽： 4. 對基因頻率作圖「R代碼：」 maf_freq <- read.table("MAF_check.frq", header =TRUE, as.is=T)
pdf("MAF_distribution.pdf")
hist(maf_freq[,5],main = "MAF distribution", xlab = "MAF")
dev.off()
「圖形如下：」 可以看出，很多基因頻率為0，說明沒有分型，這些位點需要刪掉。 4. 去掉MAF小于0.05的位點plink --bfile HapMap_3_r3_7 --maf 0.05 --make-bed --out HapMap_3_r3_8
日志： PLINK v1.90b6.5 64-bit (13 Sep 2018) www.cog-genomics.org/plink/1.9/
(C) 2005-2018 Shaun Purcell, Christopher Chang GNU General Public License v3
Logging to HapMap_3_r3_8.log.
Options in effect:
--bfile HapMap_3_r3_7
--maf 0.05
--make-bed
--out HapMap_3_r3_8
515185 MB RAM detected; reserving 257592 MB for main workspace.
1399306 variants loaded from .bim file.
163 people (79 males, 84 females) loaded from .fam.
112 phenotype values loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 112 founders and 51 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.998153.
325518 variants removed due to minor allele threshold(s)
(--maf/--max-maf/--mac/--max-mac).
1073788 variants and 163 people pass filters and QC.
Among remaining phenotypes, 56 are cases and 56 are controls. (51 phenotypes
are missing.)
--make-bed to HapMap_3_r3_8.bed + HapMap_3_r3_8.bim + HapMap_3_r3_8.fam ...
done.
可以看到，325518個位點被刪掉了，剩余1073788 個位點。 「結(jié)果文件：」 HapMap_3_r3_8.bed HapMap_3_r3_8.bim HapMap_3_r3_8.fam HapMap_3_r3_8.log
后面我們用這個文件，進行后續(xù)的質(zhì)控。 相關(guān)系列：
筆記GWAS操作流程1：下載數(shù)據(jù)
筆記GWAS操作流程2-1：缺失質(zhì)控
筆記GWAS操作流程2-2：性別質(zhì)控
|
