google深度学习用于糖尿病视障碍检测

Published:

就在11月底,google research博客,发出一篇文章,介绍其使用深度学习检测糖尿病诱发的视网膜病变的成果,该研究使用由9,963张视网膜照片组成的验证数据进行测试,结果机器给出的结果,其F-score值(综合反映敏感性和准确性的度量指数,最大为1,越大越好)为0.95,高于8为医学专家的中位数值0.91。

简单搜索了下,微信公众号已经有一篇文章介绍这个成果的,但一看就是从某个英文网站翻译过来,而且整体的逻辑也不清晰,所以我觉得有必要写一篇新文章,这里的文字并非是google research博客的原文翻译,而是在保证准确转述的基础上加了我自己的理解。

目前全球有4.15亿糖尿病人,而Diabetic retinopathy(简称DR,查了维基百科的解释,这里译为『糖尿病视障碍』或『糖尿病眼盲』)是增长最快的导致失明的发病诱因。这种病发现早的话,是可以治愈的,而发现晚的话,就无可奈何了。

目前对DR病人的诊断,是由医学专家对病人的眼睛黑区照片检查,依据损伤的存在和面积判断疾病的发生和严重度,进行这种人工检查需要专门的训练,而很多这类疾病高发地区恰恰稀缺这样的医学专家。

google和印度、美国的医生合作开发了这个算法:

  • 首先解决的是训练数据,在医生的帮助下收集整理了128,000张图片,这些图片都经过由54名医学专家组中的3-7人手工检查并做了标记
  • 接着,使用这批数据进行深度神经网路算法的训练
  • 然后使用~12,000张的两批新数据集进行算法效力的评估,而评估的金标准,大部分是来自一个7-8人的专家团队的意见,这个团队是由前述54人的专家委员会内部推举产生

结果就是下面这张图,深度学习算法的表现和专家判读结果相当,甚至优于专家组的平均值。
最后googler也提醒我们,不要高兴太早,原因是:

  • 尽管以常规的衡量方式看,深度学习算法的表现已经很好,他们认为还要收集更多更准确的参考结果,用于提升算法表现
  • 解读2D的病人眼部图片是诊断糖尿病眼疾的多个步骤中的一个,有些情况下还需结合3D成像技术,OCT来细致地检查病人视网膜的不同层。好消息是,利用机器学习解决3D成像结果的识别,已经由google的另一个部门DeepMind在做了,相信结合两个方法可以帮医生诊断更多眼部疾病

千万富翁的财富之旅

Published:

昨天花了一上午的时间,看完了刘欣的书《刘大猫的财富之旅》。读完他的故事,再对比自己看,给我很大启发。他很小的时候,就对财富有很强的渴望。在路上看到豪车,就会在心里想『很多人连一百万可支配的收入都没有,他/她却可以花一百万买一辆车,他是怎么做到的?』我自己则没有这种想法,甚至对创造财富也没有那么强烈的愿望,所以到现在为止仍然为买房苦恼纠结也很容易理解了。

他小学时候开始接触电脑,自己捣鼓了一个传奇私服,这段经历也很奇特的。发现他有个特点,做事很明确的目标,而且为此坚韧不拔的付出努力,这方面是我还不具备的。他从而接触电脑到后来自己架设传奇服务器,再到自己建网站的经历,期间学习很多技术和知识,但他坦诚自己并非对技术很感兴趣,而是出于赚钱才这么做的(可能起先架设传奇私服不是,后面建网站就确实无疑了)。我则有区别,我自己也买域名,弄过自己的博客,但就是出于兴趣的折腾。
赚到人生第一个一百万,对于自己的财富之旅无疑是一个里程碑意义的时刻。我很佩服他拿到这笔钱之后的智慧,给自己定了一个原则:

这笔钱不能用于互联网创业

他的解释是,互联网的生意要么是不需投入大笔资金,而是大量时间和精力的,就像他之前做的生意一样,要么就是需要投入很大资金去砸市场,比如滴滴出行,这一百万根本不够,而且可能自己当时也驾驭不了。基于这样的观察,他的投资决策是

投资房产
买域名
投资自己

而这三个决策后来逐一证明都是正确的,房产的收益具体没说,估计也在几十倍的量级,域名收益十倍,投资自己就是买了很多书来看,报各种学习班,另外买了一部车,这样很大程度上增加自己活动半径,在长三角地区谈生意办事都很方便。这让我产生一个感觉,有些东西,你不拥有过,就不会产生特定的想法,而没有想法,断不会有行动。

总的来说,从他身上,我学到的是,一个人要想成事儿,要有清晰的目标和准确的判断力。

Clinvar学习笔记

Published:

ClinVar是NCBI开发的一个存储人类变异和表型关系的数据库,一般会有证据支持信息(来自临床诊断还是仅仅源自文献报道)。

数据库建立

2013年发布第一个版本,初始数据由OMIM,dbSNP,GeneViewer等数据库,以及一部分临检实验室提交的数据。

数据库结构

同一个submitter不同批次提交的数据,都会被记录(前缀是SCV,不同批次间做了区分)。不同submitter提交的针对同一个位点的变异,聚合为一个RCV记录,如果他们的结论有冲突,会有一个专家组进行review后给出结论,同一个位点不一致的结论每个月会发布一个release。

数据收集方法

  • 临床检测:符合CLIA或ISO 1589标准的实验室报告出来
  • 科学研究:科研项目鉴定到的
  • 仅文献报道:被第三方引用并未声明出处,和第二点的区别是,科学研究来源的有据可查,而这种没有原文可追溯

需要注意的是,GWAS的结果,除了经人工鉴别并给出临床解释的之外,其他的都没有收录。

数据库使用

数据获取方式

  • web 可使用基因名,rs number,疾病名,HGVS表达式(人类基因组变异学会定义的用于精确描述变异信息的格式)
  • ftp XML格式存储的是全部ClinVar的信息,VCF存储的是SNP信息,大于50nt的SNV不包括在内
  • api 支持esearch efetch esummary elink共4种方法,这里是文档

sqlite学习笔记

Published:

最近在学习sqlite,sqlite官网上面的文档无疑是最权威的,但对于新手来说上手有些困难,而是比较适合已经入门的初级使用者查阅。搜索之后发现下面这两个网站,还不错,我主要是跟着第一个学习的:

简介

获取帮助 :.help
显示当前设置:.show 比如显示table时是否打印表头,各列之间以什么字符分隔,默认‘|’

数据库和列表打开与查看

打开一个数据库:.open <database_name> 通过这种方法打开的数据库,默认的名字是main,如果需要同时打开多个数据库,可以为每个数据库取一个别名,ATTACH DATABASE '/home/manager/Documents/sqlite3/ath.db' as 'ath';
注意,maintemp是保留名,在每个数据库连接中都有,用于作为primary database,存放临时列表和其他数据对象。

显示当前已打开的全部数据库使用.database,同样地,如果需要查看当前session全部的table,可使用.table(s)。如果有多个数据库,注意看到拥有alias name的数据库中的列表,会以databasename.tablename的格式显示,但是在引用他们的名字时,使用全名(databasename.tablename)和列表名均可。

查看一个表共有哪些列,.schema OrganismsPRAGMA table_info(Organisms);
注意,sqlite中commands全部是小写字母组成,且以点号’.’开头,而statements是大写字母组成,以分号结尾’;’ 你可能已经看出来了,sqlite是区分大小写的,这点和MySQL不同(关于sqlite、MySQL和NoSQL三者的比较,这里的一篇文章很值得一看)。注释分为两种,--为行内注释,有限范围从字符往后到行尾,/*是C语言风格的注释,为块注释,位于/**/之间的内容会被sqlite引擎略过。

常用的操作:简单统计、查询

类似于linux的head,有的时候需要从列表中选出几条记录看一下,可以这样SELECT * FROM Organisms LIMIT 3; 返回如下:

1
abbr        name
----------  -----------------
aaa         Acidovorax avenae
aac         Alicyclobacillus
aad         Alicyclobacillus

此外,还有其他常见操作

1
2
3
4
5
6
7
8
9
SELECT COUNT(abbr) FROM Organisms; --统计表中全部记录的数目
SELECT COUNT(abbr) FROM Organisms WHERE abbr LIKE 'a__'; --统计abbr以「a」开头,共3个字符长的记录数
COUNT(abbr)
-----------
232
SELECT COUNT(abbr) FROM Organisms WHERE abbr LIKE 'a%'; --统计abbr以「a」开头的所有记录数
COUNT(abbr)
-----------
258

SGE新增用户配置

Published:

日常的服务器管理中,经常需要添加或删除用户,对于linux用户的添加比较简单,使用useradd就可以:

1
#增加一个名为"zhangshan"的用户,并把他加入"dev"群组,同事创建其home目录(注: 这里所有的操作,都是以root用户身份进行的)
useradd -m -g dev zhangshan
groupadd qa #新建一个组
#更改密码
passwd zhangshan
	<p@ssw0rd>

添加完之后,新的用户就可以正常使用终端ssh登录集群(服务器)了,如果集群配置了SGE之类的作业调度系统,通常还需要对这些用户开通相应的集群队列权限,还是以用户”zhangshan”为例,将他加入”arusers”这个userset(SGE 称为Access List),这样他就可以使用”arusers”可调度的集群队列(all.q)了。

1
#将"zhangshan" 添加到"arusers"用户组
[root@ngnode1 cuijie]# qconf -au zhangshan arusers
added "zhangshan" to access list "arusers"
#显示"arusers"用户列表,看到"zhangshan"已经加入这个列表
[root@ngnode1 cuijie]# qconf -su arusers
name    arusers
type    ACL
fshare  0
oticket 0
entries chengchaoze,cuijie,zhouyi,leijie,zhangshan
#把"zhangshan"从"arusers"删除
[root@ngnode1 cuijie]# qconf -du zhangshan arusers
deleted user "zhangshan" from access list "arusers"

更多选项和参数,可参考qconf -help

setup and run PIPS to identify Putative Pathogenicity Islands

Published:

Recently I need to do PAI search using PIPS, after going throught all the little traps, I think it worth to write the whole process down.

  1. follow the README instructions to install hammer and tRNA-scan. Maybe you want to change the default prefix of INSTALL derectory
  2. Theoretically, it should be good enough to run, but you’ll definitely be stuck at some point. OK, here are the tricks to try below
  3. replace gbk2embl.pl with a new one, such as this:
#!/usr/bin/env perl
use strict;
use warnings;
use Bio::SeqIO;

if (@ARGV != 2) {    die "USAGE: gb2embl.pl <in.gbk> <out.embl>\n"; }

my $seqio = Bio::SeqIO->new('-format' => 'genbank', '-file' => "$ARGV[0]");
my $seqout = new Bio::SeqIO('-format' => 'embl', '-file' => ">$ARGV[1]");
while( my $seq = $seqio->next_seq) {
  $seqout->write_seq($seq)
}

For the last step, modify mergfiles.pl. change for e.g. %hash2->{$chave} to $hash2{$chave}, note that a couple of similar syntax mistake in this script.

Fast genome and metagenome distance estimation using MinHash

Published:

mash 是一个可快速计算两个(或多个)基因组之间相似度的软件,可输入基因组序列或测序结果(fastq文件),使用计算机科学中一个经典的方法MinHash,一般用来判断两个网页内容是否雷同,以决定是否从搜索结果中去掉某些类似的结果。

使用示例:

mash dist -p 4 -l 、/path/to/ref.fa query_fa.txt > query_mashed.txt

用于计算所有query序列和ref序列的相似度,而相似度和物种间的亲缘关系是相关的,其中 是一个记录了所有query序列的文件

parsing ABI file

Published:

ABI测序仪下机的结果包含.ab1和.seq文件,通常可以导入BioEdit这样的图形化软件检查峰形和运行blast,而如果想要实现自动化的检测杂合位点并做比对,最终得到变异位点,就需要花些时间研究下。

主要的问题是,ab1文件是二进制的,关于文件格式的文档稀少,唯一靠谱的是ABI的官方文档 仔细看过这个文档你会发现,看完之后作为新手的你尽管对于这个格式有所了解,知道它分为Header、Directory以及Data,以及有很多的tag。但知道这些对于想要从ab1提取trace value并进行heterogeneous base calling是远远不够的。所好的是,已经有人有过类似需求并作了尝试。

  • abifpy 实现的一个解析ab1文件的模块,现已并入BioPython。优点是兼容Python2和Python3,部分meta data已经给你提取出来了,而如果想要更多,可以使用get_data(‘youkey‘)获取。

  • 如果只需要取出每个通道的trace value,有了官方文档和abifpy就够了,如果还想学习更多,可以读一读SeqTrace源码,从中可以学到怎样组织一个项目的代码,以及怎么写测试。

  • 前面两个都研究一番之后,还是不能满足你的好奇心?可以把io lib这个库研究下,学习使用它的API.

hello again

Published:

翻了一下以前的文章发现,最近的上一篇文章是2013年11月的,那时候还在南京的学校,从找工作开始,然后面试,来北京,上班之后就更没有留出时间静下来写东西了。后来换了Mac,还特意把blog的资料作了备份,想着不久之后,重搭环境另起炉灶的,没想一下子拖到不久前才在新电脑上配好环境。

上班的这一年多,改变了很多。工作上从项目执行到流程负责人,再到信息分析主管而后技术副总监,见证了公司规模的扩张和变迁,也目睹了身边同事的来来去去,“送走”了三波运营岗同事,现在接触的大都是第四波的人。公司在成长,自己的工作也在经历变化,以前只需要自己干活,现在要考量人和事,把相应的事交给合适的人去做,然后跟踪进度,培养新人。

除了工作,生活上也有一些变化。刚到北京住的没有独立卫生间的小院,后来的院子有单独卫生间,尽管很少,还是方便多了。随公司搬到天津之后,住的更好。经过两段感情经历的洗礼,也成长了很多,对自己也有了更多地认识。写东西是一个特别好的总结的方式,通过它可以和自己对话,从而自省;和别人交流,从而结交同好。所以,以后我会坚持下去的。

retrieving and sending infomation over the web

Published:

I found this when hanging around online. To be honest, I am dying to learn web crawling long time ago, can not wait to play with this.

modules involved:

This post is very easy to follow with basic familarity of Python, the difficulty of web crawling is lots of, if not all, websites follow HTML specification exactly, there are often errors here and there if you try to parse useful bits of infomation out there, and here comes BeautifulSoup which handle most of these for us.

The second section discuss sending email with Python. One thing need to keep in mind is check your internet connection first if your code refuse to work :-)- I have complete the spec of No.5, and below are some commands I find useful:

  • execfile('script.py') #load script into console
  • python -i script.py #execute script then come to Python console, ‘i’ short for ‘interactive’
  • sometimes you need to insert a time stamp into you file, and time.ctime() is probably the simplest solution
    source code is here