大家好,我是鄧飛。 GWAS分析完成后,進行單倍型圖分析的核心目的是驗證顯著性位點的可靠性并深入理解其遺傳背景,具體原因包括以下幾點: 排除假陽性結果 GWAS發現的顯著性位點可能因多重檢驗或群體結構干擾產生假陽性。單倍型圖通過分析顯著性位點上下游區域的連鎖不平衡(LD)結構,判斷這些位點是否處于高度連鎖的區塊(Block)中。若基因位于Block內且被高LD區域覆蓋,則說明該基因的關聯性更穩定,假陽性風險顯著降低。 識別連鎖區域與單倍型結構 單倍型圖可直觀展示SNP之間的LD關系,劃分Block(如通過D’或r2閾值)。Block內的SNP組合形成少數單倍型,這些單倍型代表共同遺傳的DNA片段。若目標基因或標記位于Block內,其功能或調控作用更可能因連鎖效應而真實存在。 輔助基因定位與功能驗證 當GWAS位點位于非編碼區時,單倍型圖可幫助鎖定候選基因的調控區域。例如,若Block覆蓋某基因的啟動子區,可推測該基因通過調控表達水平影響表型,從而指導后續實驗驗證方向。 優化分子標記選擇 在分子育種或醫學遺傳學中,單倍型圖可篩選出Block內代表性SNP(如標簽SNP),減少冗余檢測成本,同時保留遺傳信息完整性。 1,比較流行的圖,Block + GWAS pvalue + LD result- 最下方的熱圖是兩兩SNP之間的LD值,越高越紅,比較紅的區域構成一個Block(用黑線連起來)
- 如果提供gff文件,可以顯示基因的上游、下游、外顯子、內含子區域
- 位點之間,也可以根據LD值進行可視化,以最顯著的位點為四方形,其它位點與其LD值的大小呈現不同的顏色
參考代碼:
LDBlockShow -InVCF Test.vcf.gz -OutPut re2 -Region chr11:24100000:24200000 -InGWAS gwas.pvalue -OutPng -SeleVar 1
結果:
re2.blocks.gz re2.png re2.site.gz re2.svg re2.TriangleV.gz

ShowLDSVG -InPreFix re2 -OutPut temp -InGWAS gwas.pvalue -Cutline 7 -ShowNum -PointSize 3
命令:
LDBlockShow -InVCF Test.vcf.gz -OutPut re3 -Region chr11:24100000:24200000 -InGWAS gwas.pvalue -OutPng -SeleVar 1 -InGFF In.gff
 也可以增加SNP的名稱:
$ cat Spe.snp chr11 24142660 chr11 24142669 SpeA chr11 24142760 SpeB
命令:
LDBlockShow -InVCF Test.vcf.gz -OutPut re3 -Region chr11:24100000:24200000 -InGWAS gwas.pvalue -OutPng -SeleVar 1 -InGFF In.gff -SpeSNPName Spe.snp
2,比較流行的圖,Haploview + 箱線圖 + 顯著性分析 + 單倍型頻率統計進一步分析呢?比如將單倍型block看一下不同樣本的分布,統計一下不同單倍型和表型數據的關系,做一個顯著性分析不是更有說服力嗎?特別是找到的單倍型應用于育種時,通過單倍型的顯著性分析更有說服力:優勢單倍型,優良基因,呼之欲出,一圖勝千言,這樣更有說服力。3,對單倍型和表型數據進行顯著性檢驗,如果是兩個單倍型就用T檢驗,如果是3個及以上,就用方差分析4,對結果可視化就是上面圖的結果了,有時候還有小提琴圖的樣式:發現,10號染色體的,這個區間,有幾個位點位于一個block里面。## 導入基因型vcf數據 library(geneHapR)
vcf = import_vcf("df2.vcf")
# 單倍型分型 hapResult <- vcf2hap(vcf)
write.csv(hapResult,"df2-1-hapresult.csv")
上面結果中,共有11個單倍型,每個個體都會給出具體的單倍型分型。表型數據準備,兩列數據,第一列是ID,第二列是分析的性狀,整理成txt文件,tab分割。使用函數import_AccINFO,讀取txt文件,自動將第一列ID變為行名。結果文件中,可以看出,共有兩個單倍型:H001和H002,其中H001有10個個體,H002有8個個體,觀測性狀為y1,兩個單倍型之間沒有顯著性。data("geneHapR_test")# plot the figs directlyhapVsPheno(hap = hapResult, pheno = pheno, phenoName = "GrainWeight.2021", minAcc = 3)
|