移步 github 實(shí)現(xiàn)目錄跳轉(zhuǎn)试幽，獲得更好的閱讀體驗

目錄

• 1. 基因型填充
  • 1.1. 問題描述
  • 1.2. 技術(shù)來源的基因型缺失
  • 1.3. 缺失的判斷：缺失率
  • 1.4. 基因型缺失的影響
  • 1.5. 基因型填充的原理
  • 1.6. 實(shí)現(xiàn)工具
• 2. 實(shí)操：用IMPUTE2實(shí)現(xiàn)基因型填充
  • 2.1. 兩種應(yīng)用場景
  • 2.2. Best Practices
  • 2.3. 1000 Genomes Imputation Cookbook
    • 2.3.1. Before Imputation
    • 2.3.2. Pre-Phasing
      • 2.3.2.1. using IMPUTE2
      • 2.3.2.2. using SHAPEIT (recommended)
    • 2.3.3. Imputation

1. 基因型填充

1.1. 問題描述

基因型缺失：樣本中沒有被測序數(shù)據(jù)覆蓋到的區(qū)域状飞，基因型就屬于未知的，我們將之稱為缺失位點(diǎn)

基因型數(shù)據(jù)的缺失又分為遺傳性缺失和檢測性缺失：

• 遺傳性缺失：個體遺傳信息的變異（例如谆膳，這個位點(diǎn)DNA片段真實(shí)缺失）導(dǎo)致的基因型缺失

• 檢測性缺失：由于檢測技術(shù)的局限懂衩、錯誤等導(dǎo)致的信息丟失华坦。各類基因型檢測技術(shù)都會產(chǎn)生檢測性的基因型缺失愿吹。

1.2. 技術(shù)來源的基因型缺失

1. 全基因組重測序技術(shù)

   全基因組重測序理論上應(yīng)該覆蓋整個基因組，因此未覆蓋的區(qū)域都可以被定義為缺失惜姐。那么群體研究中的低深度測序（一般平均深度低于10X）犁跪，不可避免會產(chǎn)生大量隨機(jī)缺失。

2. 簡化基因組測序

   簡化基因組測序是通過酶切歹袁，并富集限制性內(nèi)切酶周邊的片段并進(jìn)行測序的策略坷衍。針對簡化基因組，我們稱的缺失一般指的是沒有被檢測到的酶切片段相關(guān)的位點(diǎn)条舔。簡化基因組的缺失枫耳，主要與酶切效率有關(guān)。酶切效率越高孟抗，缺失率越低迁杨。

3. 外顯子測序以及目標(biāo)區(qū)域捕獲測序

   同簡化基因組測序類似，基于探針雜交的DNA捕獲以及測序技術(shù)凄硼，同樣會產(chǎn)生大量的缺失铅协。這種缺失主要是由于探針雜交捕獲的效率所致。

4. SNP芯片

   SNP芯片利用芯片雜交后的熒光信號摊沉，來判斷某個位點(diǎn)的基因型狐史。SNP芯片同樣也會產(chǎn)生大量缺失。但在實(shí)際的研究中说墨，SNP 芯片主要面臨的問題是芯片型號不同预皇，甚至來源不同的廠商，那么芯片中包含的SNP位點(diǎn)也不同婉刀。當(dāng)來源不同的數(shù)據(jù)一起分析的時候，將面臨數(shù)據(jù)不一致的問題序仙。簡單說來突颊，就是你有的我沒有，我有的你沒有潘悼。如下圖律秃，Affymetrix和illuminate兩大SNP 芯片廠商生產(chǎn)的人類芯片就使用的是不同的SNP集，當(dāng)放在一起分析的時候?qū)⒚媾RSNP不一致的問題治唤。

注意棒动！

基因型缺失是一個相對性的概念。以上缺失的概念都是針對同種技術(shù)的比較宾添。不同的技術(shù)比較船惨，也可以定義為缺失。例如，同樣一份樣本扁眯，我們使用全部以上4種技術(shù)檢測钮呀。如果以全基因組高深度測序（>30X）為參照標(biāo)準(zhǔn)，后續(xù)的3種技術(shù)都有大量位點(diǎn)沒有檢測到怜浅，處于基因型缺失的狀態(tài)铐然。

1.3. 缺失的判斷：缺失率

分為樣本水平的缺失率和位點(diǎn)水平的缺失率

例如下圖，0恶座、1搀暑、2 分別代表三種檢測到的基因型，圖中缺失位點(diǎn)使用“跨琳？”表示自点。那么樣本1的缺失率=20%（總體10個位點(diǎn)，有兩個位點(diǎn)缺失）湾宙，而位點(diǎn)2的缺失率=60%（總體5個位點(diǎn)樟氢，有3個位點(diǎn)缺失）

1.4. 基因型缺失的影響

基因型缺失最直接的影響就是這個位置的信息缺失，從而影響下游分析（包括遺傳圖譜構(gòu)建侠鳄，QTL定位埠啃，選擇壓力分析，GWAS分析等）的信息完整性和準(zhǔn)確性伟恶。

例如碴开，（b）中紅色的點(diǎn)是（a）中缺失的位點(diǎn)。而與性狀關(guān)聯(lián)的SNP位點(diǎn)博秫，恰恰位于虛線所在的區(qū)域內(nèi)潦牛。這些顯著位點(diǎn)在（a）中是缺失的，所以（a）沒有檢測到關(guān)聯(lián)信號挡育，從丟失了非常關(guān)鍵的信息

基因型缺失對GWAS分析巴碗、選擇壓力分析影響都比較大

1.5. 基因型填充的原理

原理：

基于家系樣本的遺傳特性。具有已知親緣關(guān)系的個體之間具有共享的單體型（haplotype）即寒，這些由有限個遺傳標(biāo)記所構(gòu)成的單體型隨祖先一起遺傳橡淆，反映連鎖不平衡。

在具有相同單體型的家系中母赵，遺傳標(biāo)記少的樣本可以參照遺傳標(biāo)記多的樣本進(jìn)行基因型填充逸爵。對于沒有親緣關(guān)系的樣本，以上理論也基本適用凹嘲，主要的差別在于無血緣關(guān)系的樣本之間共享的單體型比家系樣本之間的要短很多师倔。對無親緣關(guān)系樣本進(jìn)行基因型填充需要一個高密度遺傳標(biāo)記構(gòu)成的單體型圖譜作為參照。

通過對比待填充樣本和參考模板周蹭，找到兩者之間共有的單體型趋艘，然后就可以將匹配上的參考模板中的位點(diǎn)復(fù)制到目標(biāo)數(shù)據(jù)集中疲恢。

常見imputation的基本邏輯包括兩步：

• 從目標(biāo)位點(diǎn)/區(qū)域非缺失的位點(diǎn)中，總結(jié)這個區(qū)域的基因型規(guī)律致稀，并分類冈闭。其實(shí)就是分析各個區(qū)域的單體型組成；

• 根據(jù)某樣本缺失位點(diǎn)的上下其他非缺失位點(diǎn)抖单，判斷這個區(qū)域?qū)儆谀姆N單倍型萎攒。然后根據(jù)所屬單倍型的基因型補(bǔ)充該樣本的缺失位點(diǎn)；

根據(jù)缺失樣本有限的基因型信息（僅有3個位點(diǎn)）矛绘，就可以判斷這個樣本與參考單倍型集中的哪種單倍型最為相似（圖中分別對應(yīng)紫色耍休、綠色、黃色三種單倍型）货矮。然后羊精，將對應(yīng)的最相似的單倍型賦予給該樣本，從而讓該樣本獲得完整的基因型囚玫，圖b

1.6. 實(shí)現(xiàn)工具

(1) 計算密集型喧锦，比如IMPUTE、 IMPUTE2抓督、MACH燃少、 和fastPHASE/BIMBAM

這種類型的方法在填充的過程中充分考慮到全部可以觀察到的基因型信息，使得對缺失值的估算更加精確铃在；但以上大部分軟件都是針對人類的開發(fā)的阵具。人類種群的遺傳特性是個體雜合率較高、近交率低定铜、系譜關(guān)系來源隨機(jī)阳液。很多植物，尤其作物的遺傳特性則和人類相反揣炕。

(2) 計算高效型帘皿，比如PLINK、TUNA畸陡、WHAP和BEAGLE

此種算法僅僅關(guān)注與特定位點(diǎn)相鄰的一小部分標(biāo)記的基因型矮烹，因此在計算上更加快捷

2. 實(shí)操：用IMPUTE2實(shí)現(xiàn)基因型填充

2.1. 兩種應(yīng)用場景

Impute2的基因填充 (genotype imputation) 分為兩種應(yīng)用情景：

1. ONE REFERENCE PANEL

   ./impute2 \
    -m ./Example/example.chr22.map \
    -h ./Example/example.chr22.1kG.haps \
    -l ./Example/example.chr22.1kG.legend \
    -g ./Example/example.chr22.study.gens \
    -strand_g ./Example/example.chr22.study.strand \
    -int 20.4e6 20.5e6 \
    -Ne 20000 \
    -o ./Example/example.chr22.one.phased.impute2
   

   參數(shù)說明：

   • -m <file>: 目標(biāo)區(qū)域重組率圖譜文件(Fine-scale recombination map for the region to be analyzed)，記錄的是基因組中各個位點(diǎn)的重組率和彼此間物理距離的關(guān)系

   這個文件應(yīng)該包含三列：

   (1) physical position: in base pairs
   (2) recombination rate: between current position and next position in map (in cM/Mb)
   (3) genetic map position: in cM
   
   例如：
   position COMBINED_rate(cM/Mb) Genetic_Map(cM)
   35326 0.251801 0.000000
   35411 0.482009 0.000021
   40483 0.598191 0.002466
   

   ---

   • -h <file 1> <file 2>: 已知的單體型信息文件罩锐，每行表示一個SNP位點(diǎn)，每列表示一個單體型 (one row per SNP and one column per haplotype)

   所有的allele必須表示成0或1的形式

   一旦用-h參數(shù)指定一個單體型文件卤唉，就需要用-l參數(shù)指定一個對應(yīng)的Legend文件

   Impute2允許同時指定兩個單體型文件： -h <file 1> <file 2>

   ---

   • -l <file 1> <file 2>：與單體型文件對應(yīng)的Legend文件涩惑，保存的是對每個SNP位點(diǎn)的描述信息

   這個文件包含四列：

   rsID, physical position (in base pairs), allele 0, and allele 1
   
   最后兩列的 allele 0 和 allele 1 是對堿基組成的說明
   

   ---

   • -g <file>: 包含目標(biāo)研究群體的genotypes的文件，即Genotype File Format桑驱，對它進(jìn)行后續(xù)的基因型填充 (impute) 和分型 (phase)

   該文件每行表示一個SNP竭恬，前五列分別為：

   (1) SNP ID：這一列一般表示為染色體號
   (2) RS ID of the SNP
   (3) base-pair position of the SNP
   (4) the allele coded A
   (5) the allele coded B
   

   緊接著的3列是群體中的一個個體的三種可能的基因型：AA跛蛋，AB或BB

   再接著3列是第二個個體的，以此類推痊硕，示例文件如下：

   SNP1 rs1 1000 A C 1 0 0 1 0 0
   SNP2 rs2 2000 G T 1 0 0 0 1 0
   SNP3 rs3 3000 C T 1 0 0 0 1 0
   SNP4 rs4 4000 C T 0 1 0 0 1 0
   SNP5 rs5 5000 A G 0 1 0 0 0 1
   

   ---

   • -strand_g <file>: 指定SNP所在的鏈的方向

   該文件每行表示一個SNP赊级，包含兩列（列之間用單個空格隔開）：（1）SNP所在的堿基位置；（2）鏈的方向岔绸，+或-理逊；

   • -int <lower> <upper>: 用于基因型推斷的基因組間隔的長度，可以以長格式表示盒揉，如 -int 5420000 10420000晋被，也可以以指數(shù)形式表示，如 -int 5.42e6 10.42e6

   • -Ne <int>: 這個參數(shù)的說明看不懂刚盈，把原文貼在下面：

   "Effective size" of the population (commonly denoted as Ne in the population genetics literature) from which your dataset was
    sampled. This parameter scales the recombination rates that IMPUTE2 uses to guide its model of linkage disequilibrium patterns.
   When most imputation runs were conducted with reference panels from HapMap Phase 2, we suggested values of 11418 for imputation 
   from HapMap CEU, 17469 for YRI, and 14269 for CHB+JPT. 
   

   • -o <file>: 輸出文件名羡洛，文件格式與-g參數(shù)指定的文件相同，即都是Genotype File Format

2. TWO REFERENCE PANELS

   在這種應(yīng)用情景中藕漱，用到了兩個refrence panel欲侮，分別記作 panel 0 和 panel 1

   例如，panel 0 可以是1000 Genomes Project的haplotype肋联，包含了基因組中幾乎全部的常見SNPs威蕉；panel 1 可以是HapMap Phase 3的haplotype，僅包含了基因組中的部分的常見SNPs牺蹄；panel 3 是用商用SNPs芯片得到的一系列的case和control的樣本的genotype

   之所以使用兩套reference panels忘伞，是想通過兩個reference panels互相填充對方中缺失的部分來得到SNP密度更高的一套merged reference panel

2.2. Best Practices

1. 基因型填充前 (pre-imputation) 進(jìn)行g(shù)enotypes質(zhì)控

   過濾低質(zhì)量的變異位點(diǎn)和樣本

   質(zhì)控方法可以參照：

   Anderson, C.A. et al. Data quality control in genetic case-controlassociation studies. Nat. Protoc. 5, 1564–1573 (2010)

2. 保證分析中使用的基因組坐標(biāo)系統(tǒng)一致

   NCBI build number (e.g., "b36" or "b37") 對應(yīng)于 UCSC version (e.g., "hg18" or "hg19")

3. 選擇reference panel

   之前的GWAS研究中，研究人員一般都是選擇與對應(yīng)人群遺傳距離最相近的reference panel沙兰，而Impute2推薦使用worldwide reference panel氓奈，程序能夠從中選出最合適的haplotype用于基因型填充

   這樣做的好處是：

   • 不需要費(fèi)時費(fèi)力去挑選haplotypes來構(gòu)造reference panel；

   Good results can be obtained in any study population by tuning a single software parameter (-k_hap) with a simple rule of thumb
   

   • 該策略適用于各種人群的研究鼎天；

   Our group and others have used this approach to successfully impute populations ranging from homogeneous isolates to recent and
   complex admixtures
   

   • 填充效果往往優(yōu)于研究人員自己挑選構(gòu)造的小reference panel

   This is because individuals from "diverged" populations may still share genomic segments of recent common ancestry, and IMPUTE2
    can use this haplotype sharing to improve accuracy. At the same time, the software can ignore haplotypes that are not helpful.
   

   • 對于大的reference panel舀奶，Impute2也能進(jìn)行高效地處理，不需要擔(dān)心會帶來的計算負(fù)擔(dān)

4. 填充后 (post-imputation) 質(zhì)控

   質(zhì)控方法可以參照：

   Verma, S.S. et al. Imputation and quality control steps for combining multiple genome-wide datasets. Front. Genet. 5, 370 (2014).

2.3. 1000 Genomes Imputation Cookbook

2.3.1. Before Imputation

1. 對Genotype數(shù)據(jù)進(jìn)行質(zhì)控

   包括樣本水平的質(zhì)控和marker水平的質(zhì)控

   • 樣本水平：
     • call rate
     • heterozygosity
     • relatedness between genotyped individuals
     • correspondence between sex chromosome genotypes and reported gender
   • marker水平：
     • call rates
     • deviations from Hardy-Weinberg Equilibrium
     • excluding low frequency SNPs斋射，for older genotyping platforms,

   質(zhì)控代碼參考自：http://sites.google.com/site/mikeweale/software/gwascode

2. 將Genotype數(shù)據(jù)轉(zhuǎn)換為Build 37

   目前的1000 Genome Project的數(shù)據(jù)使用的是NCBI genome build 37 (hg19)的坐標(biāo)系統(tǒng)育勺，因此在基因型填充之前需要保證你的Genotype文件也是hg19的坐標(biāo)系統(tǒng)，且位點(diǎn)是落在正鏈上

   若坐標(biāo)系統(tǒng)不一致罗岖，可以使用LiftOver進(jìn)行坐標(biāo)轉(zhuǎn)換涧至，但是轉(zhuǎn)換過程中可能有少量的SNP轉(zhuǎn)換失敗

3. 將Genotype文件轉(zhuǎn)換為IMPUTE格式

   在格式轉(zhuǎn)換之前，需要先按照坐標(biāo)進(jìn)行排序

   GTOOL可以將PLINK PED轉(zhuǎn)換為IMPUTE格式

2.3.2. Pre-Phasing

對于大規(guī)模的reference panels桑包，基因型填充建議分兩步進(jìn)行：

• pre-phasing：推斷每個樣本的單體型
• imputation：對分型得到的單體型 (phased haplotypes) 中缺失的allele進(jìn)行基因型填充

IMPUTE2 或 SHAPEIT 都可以執(zhí)行pre-phasing操作南蓬，Drs. Bryan Howie 和 Jonathan Marchini推薦使用SHAPEIT進(jìn)行pre-phasing，因為該工具采用的phasing方法更準(zhǔn)確

2.3.2.1. using IMPUTE2

IMPUTE2的pre-phasing推薦采用滑動窗口法 (Sliding Window Analyses) 進(jìn)行：

將一個染色體分割成若干Mb的塊 (blocks)，對這個塊中的genotypes進(jìn)行phasing赘方。分塊需要通過設(shè)置-int <start> <end>參數(shù)實(shí)現(xiàn)

用這種phasing方法可能遇到的兩種棘手的情況：

• 若每個塊的大小一致烧颖，可能不同的塊之間SNP的分布密度存在很大差異；
• 若某個塊正好跨過著絲點(diǎn) （centromere）,而著絲點(diǎn)附近區(qū)域的SNP密度極低窄陡，往往意味著是 large gap in 1000 Genome SNPs炕淮；

若SNP密度太低是很難進(jìn)行phasing的，此時可以采取的解決策略是：

• 若某一個塊的SNP密度太低（例如少于200個）跳夭，可以將它與鄰居的塊合并成一個更大的塊一起phasing涂圆；
• 避免構(gòu)造跨著絲點(diǎn)的塊；

$ impute2 \
 -prephase_g \
 -m ./Example/example.chr22.map \
 -g ./Example/example.chr22.study.gens \
 -int 20.4e6 20.5e6 \
 -Ne 20000 \
 -o ./Example/example.chr22.prephasing.impute2 \
 -allow_large_regions

2.3.2.2. using SHAPEIT (recommended)

對整條染色體進(jìn)行pre-phasing建議使用SHAPEIT

SHAPEIT接受PLINK PED和IMPUTE的格式輸入

$ shapeit \
 --input-bed gwas.bed gwas.bim gwas.fam \
 --input-map genetic_map.txt \
 --thread 8 \
 --effective-size 11418 \
 --output-max gwas.phased.haps gwas.phased.sample \
 --output-log gwas.phasing.log

2.3.3. Imputation

$ impute2 \
 -use_prephased_g \
 -m ./Example/example.chr22.map \
 -h ./Example/example.chr22.1kG.haps \
 -l ./Example/example.chr22.1kG.legend \
 -known_haps_g ./Example/example.chr22.prephasing.impute2_haps \
 -strand_g ./Example/example.chr22.study.strand \
 -int 20.4e6 20.5e6 \
 -Ne 20000 \
 -o ./Example/example.chr22.one.phased.impute2 \
 -phase

參考資料：

(1) 【簡書】群體遺傳學(xué)習(xí)筆記-基因型缺失數(shù)據(jù)的填充

(2) Impute2官方文檔

(3) Genotype File Format

(4) IMPUTE2: 1000 Genomes Imputation Cookbook

(5) Weale M (2010) Quality Control for Genome-Wide Association Studies. Methods Mol. Biol. 628:341–372

(6) van Leeuwen EM, et al. Population-specific genotype imputations using minimac or IMPUTE2[J]. Nature Protocols, 2015, 10(9):1285-1296.