為什麼會有這個需求呢？很簡單，因為國內的一些基因檢測公司支援匯入23andme的晶片資料做解讀，而我正想看看一下他們的技術功底到底如何？

23andme和wegene都是用的一款特製的晶片，可以捕獲基因組上面的一些特定位點而已，既然我已經測了全基因組，那麼分分鐘就可以從我的基因組分析結果裡面提取出23andme的晶片位點，偽裝成23andme的晶片資料！

我從谷歌裡面找到了一個公共的資料，點選閱讀原文檢視這個公共資料的下載連結！

這很容易明白23andme的晶片資料是什麼格式的，基因組座標的轉換對我等生物資訊學工程師來說比吃飯還簡單！（當然，我其實拿到了新版的資料，但是由於隱私問題，不便傳播）

轉換很簡單：

第一步，把晶片設計的rsID全部拿出來

第二步，根據rsID從我的VCF檔案中挑取位點，並賦予純合雜合基因型

第三步，去dbSNP資料庫檔案裡面對映我VCF檔案沒有記錄的點為野生型

(perl -alne '{print if /^rs/}' dm_23andme_v3_110219.txt  |cut -f 1 >23andme.rsID.listcat ../variation/autochr.highQuali.dbsnp.vcf  23andme.rsID.list |perl -alne '{if($F[2]=~/^rs/){if(/1/1/){$gt=$F[4].$F[4]}else{$gt=$F[3].$F[4]};$h{$F[2]}="$F[0]t$F[1]t$gt" }  print "$_t$h{$_}" if /^rs/}' >my_23andme.1.txtzcat ~/annotation/variation/human/dbSNP/All_20160601.vcf.gz |perl -alne 'BEGIN{ open FH,"my_23andme.1.txt";while(<FH>){chomp;@F=split;if(/^rs/){ $pos{$.}=$_;if($F[3]){$h{$F[0]}=$_}else{$tmp{$F[0]}=1}  }} }{if(exists $tmp{$F[2]}){ $tmp{$F[2]}="$F[0]t$F[1]t$F[2]$F[2]"  }}END{foreach(sort{$a<=>$b} keys %pos){ if(exists $h{$pos{$_}} ){$value=$h{$pos{$_}}}else{$value=$tmp{$pos{$_}} } ;print "$pos{$_}t$value" }}'
 

wegene的晶片資料在格式上是一模一樣，因為他們用的都是illumina公司出品的定製化晶片。

本來是想上傳這個公共資料去這個網站上面做一次免費報告生成，但是他們要求很多，搞了好幾次還沒成功，最後還是嫌棄我晶片版本太低了，所以我又用了下面的程式碼把舊基因組版本晶片資料轉換成新的。

zcat ~/annotation/variation/human/dbSNP/All_20160601.vcf.gz |perl -alne 'BEGIN{ open FH,"dm_23andme_v3_110219.txt";while(<FH>){chomp;@F=split;if(/^rs/){$pos{$.}=$_;$h{$F[0]}=$F[3]} } }{if(exists $h{$F[2]}){ $h{$F[2]}="$F[0]t$F[1]t$h{$F[2]}"  }}END{print "$pos{$_}t$h{$pos{$_}}" foreach sort{$a<=>$b} keys %pos}' >dm_23andme_v3_hg19.txt
 

這個難度有點高,程式設計功底不夠就不用看了，想看看具體是怎麼回事，點選閱讀原文檢視!

參考連結：

https://www.wegene.com/demo/

https://www.mygene.com/demo

http://online.cambridgecoding.com/notebooks/cca_admin/genetic-ancestry-analysis-python

文：Jimmy

圖文編輯：吃瓜群眾