【直播】我的基因组59:把我的数据伪装成23andme或wegene的芯片数据

为什么会有这个需求呢?很简单,因为国内的一些基因检测公司支持导入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

图文编辑:吃瓜群众

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-02-23

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏Crossin的编程教室

【每周一坑】房贷计算器 |【解答】生成九宫格图片

因为是“刚需”,所以网上早有无数的版本。有人已经用过,有人以后可能会用。有没有想过,类似这种小工具,其实你自己也可以实现。

962
来自专栏FreeBuf

GSM Hacking Part ②:使用SDR捕获GSM网络数据并解密

本文作者:雪碧0xroot@漏洞盒子安全团队 0×00 在文章第一部分 GSM Hacking Part ① :使用SDR扫描嗅探GSM网络搭建了嗅探GSM流量...

2048
来自专栏腾讯大讲堂的专栏

朋友圈遇到乱码后怎么办

最近,在朋友圈里,经常会看到这样的乱码,举个例子: මම ඔබට කියන්න අවශ්ය, මම ඔයාට ආදරෙයි ไม่ว่าจะเกิดอะไรขึ...

1639
来自专栏不二小段

【一起学Python】爬取网易云歌词

说在前面:这是公众号第一篇来自小伙伴的投稿。我之前挖过一个坑,说想抓取歌词以后做文本分析,后面不了了之了。刚好Ricky作为爬虫的初学者,需要小项目练手,他就把...

33910
来自专栏玉树芝兰

如何用Sikuli自动录入成绩?

手里明明有一份学生成绩Excel表格,却还得一一手动把它们输入到教务系统?类似这样的简单重复枯燥操作,其实你都可以一键让电脑自动替你完成。

562
来自专栏Youngxj

EMLOG博客添加OwO表情教程

1746
来自专栏IMWeb前端团队

前端优化的技巧

最近有几个搭档开通了自己博客,但却诉苦说因为的买的虚拟空间,所以自己博客网站翻开速度很慢。关于这种景象,依照一般的状况来看,一个网 最近有几个搭档开通了自己博客...

1700
来自专栏数据小魔方

think-cell chart系列20——使用建议及附加功能

今天是think-cell chart系列收尾篇——使用建议及附加功能。 由于think-cell chart图表插件是office平台的第三方插件,而且图表...

3324
来自专栏工科狗和生物喵

【闲来无事,py写game】一个问答游戏Trivia -来自《Python游戏编程入门》

正文之前 嗯,没错,我只是为了规范化,就写这么多了!要洗澡了,明早有事! ? 正文 Trivia是一款书籍阅读类软件,支持Android 2.3.3。具体的...

3577
来自专栏Jerry的SAP技术分享

SAP UI 搜索分页技术

搜索分页技术往往和另一个术语Lazy Loading(懒加载)联系起来。今天由Jerry首先介绍S/4HANA,CRM Fiori和S4CRM应用里的UI搜索分...

1843

扫描关注云+社区