Use docker in nextflow pipelines

Nextflow可以很方便地搭建数据分析流程,而且对docker,云和集群计算环境也有很好的兼容。这里仅记载我使用nextflow搭建整合了docker功能的分析流程中,常用的命令。

docker相关常用命令

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
## build image
docker build -f ./Dockerfile .
## add a tag
docker build -t user/repo:tag .
## push local docker image to the remote repo
docker tag <md5_of_docker_image> user/repo:tag
docker push user/repo:tag
## one line to build an image with predefined name
docker build -t user/repo:tag -f Dockerfile .

## run the docker
docker run -i -t user/repo:tag

## manage the images
### list all the images
docker images
### remove a specific image, might need to add -f to force the operation
docker image rm <md5_of_image>

## export
docker save hifibio/pipeline:scRNAseq_nf_v1.06 | gzip -c > scRNAseq_nf_v1.06.tar.gz

## import
gunzip -c scRNAseq_nf_v1.06.tar.gz | docker import
## or import after gunzip
docker import scRNAseq_nf_v1.06.tar

对于docker镜像导入和导出,可参考这里save and export docker image

nextflow流程启动

1
2
3
##这里默认 nextflow.nf 在当前目录,如果不在,需通过 -c/-C 指定配置文件(-C为覆盖默认参数)
## -with-tower 为使用tower的后台监控功能,开启此功能需在 https://tower.nf/ 注册账号,并将自己的token写到配置文件的 tower section
nextflow run main.nf -with-docker -N user@example.com -with-tower -name name_of_this_run

email通知

如果使用-with-tower可直接打开tower内置的邮件通知,这样流程全部任务结束后,会收到邮件通知。如果需要在不用tower的情况下收到流程运行相关的通知,可以使用如下配置:

1
2
3
4
5
6
7
8
9
10
mail {
from = 'user@163.com'
smtp.host = 'smtp.163.com'
smtp.ssl.enable = true
smtp.port = 465
smtp.user = 'user'
smtp.password = 'passw0rd'
smtp.auth = true
debug = true
}

docker file permission

由于docker流程是以root身份运行,输出文件所有者是root,这样可能会给后续文件修改和清理带来麻烦,要避免此类问题,可在配置中指定docker.fixOwnershiptrue
注意对process的参数修饰,支持通过process名或通配符的方式,这体现了nextflow的灵活性,这里也支持通过指定label的方式完成,具体可参见官方文档

1
2
3
4
5
6
7
8
9
10
11
process {
container = 'aafe3c024859'
docker.enabled = true
docker.fixOwnership = true

withName: 'runFastQC.*' { cpus = 6 }
withName: 'clusterBarcodes.*' { cpus = 62 }
withName: mapTranscripts { cpus = 62 }
withName: createMatrix { cpus = 4 }
withName: createConsensus { cpus = 62 }
}

如何剪辑flv格式视频

参加公司年会节目,我们部门准备的是一个舞蹈,因为选的舞蹈视频是来自一个在线视频网站,这种视频格式是flv格式,年会上表演需要用的音乐伴奏需要剪辑其中的一部分,而手里可直接用的软件iMovie剪辑不支持flv格式。

经过网上搜索,采用以下方法解决:

  • flv视频转码为avi格式
  • 用iMovie剪辑并提取音频

flv转码avi采用ffmpeg完成。官网提供Mac下编译好的程序,使用也很简单: ffmpeg -i input.flv output.avi

你好2021

好久没有写东西了,上次听《捕蛇者说》laixintao提到写作,他形容自己就是喜欢写文章,以至于如果很久没写,会感到不舒服,而这种不舒服直到写完一篇文章发出来,才会缓解,才会重新感到浑身畅快。刚才看了下,我这的文章主要是是读书期间积累的,工作后5年积累的不足1/3。今天终于新的电脑上配置好环境,以后要多写。

很多人的生活在2020年都被深刻的改变了。我在这一年,搬了家,换了工作,从一段失败的婚姻中挣扎着走出来。社会的变化和个人生活上的大事纠缠着,驶过2020的年轮。细想起来,2020也有不少值得庆祝的事情。

  • 年底入职一家心仪的公司,之前一直想到从事创新药的公司试一试,这下终于如愿。新的同事很好,新的工作有挑战,但迎接挑战的过程也有激动和满足
  • 从2019年底开始做的理财,经过2020年市场动荡洗礼,算下来收益也还不错,接下来继续学习,交了一些市场的学费后,以后要牢牢记住敬畏市场,敬畏风险。
  • 本来心心念念的,盼望年会能中一个switch奖品的,没想到中了一个更大的一等奖iPhone12 Pro Max。

下面记录一下2020读的书,新买的电子产品,用过的数字服务。
2020年读的书不多,发现自己喜欢上了微信读书,买了一年的无限卡,用的最多的书听书功能,而且喜欢的公众号文章,也可以用微信读书聆听。年底的时候,把手里不怎么用的iPad补足差价,换了一台boox,督促自己以后多用它读书。

这一年在微信读书听完了10本书,最近读的一本是陆铭的大国大城, 核心观点就是中国应该顺应城市化和大城市规模增加的规律,不应该用政策和户籍制度对人口的流动施加限制。传言也是陆教授这本书和他的不断呼吁,才使得近两年政府的人口落户的限制有所转向和中大城市落户的放松。

陪伴了我6年的Macbook Air因为电池老化,把它回收了之后,新买了一台Mac mini,另外配了一台4K显示器。当时在Mac mini和iMac之间有过犹豫,后来想到iMac的大屏显示器是不错,但不能做外界显示器用,而Mac mini + 4k显示器,虽然价格和iMac接近,但4显示器还外接windows PC作扩展用,这样更方便。Mac mini的HDMI接口不太好用,经常从待机唤醒的时候,屏幕不能点亮,试了更换显示器,没用。网上查过之后,发现别人也遇到同样问题,解决方法是使用type C转HDMI的转接头连接。更换好,再也没有出现过上述屏幕点不亮问题。

记一次午餐会

今天下班走在路上,一阵风吹过过,一股凉意袭来,不由自主地打了一个哆嗦。不觉间秋天已经悄悄地爬上树梢,拨弄下片片红叶。每当这时,总会想起大学好友Ben,想到他面对此情此景,幽幽地吐出三个字:“秋—深了”。

中午公司午餐会,照例请到了Gary来给大家打气,他是创始人里演讲能力最强的,也难怪每次公司活动基本都会请他为大家讲话,这次一样不负众望。不过还好,说好讲20分钟,没有拖堂,不然大家肚子都要咕咕叫了。Gary这次先从自己喜得老五宝贝女儿,克服很多困难,说到创业往往九死一生,勉励大家鼓足干劲,迎接挑战。回想起来,Gary每次跟大家讲话,主旨大体都是号召大家不畏困难,拼搏奋斗,但是每次讲的内容,用的例子都不一样,几乎每次都做到联系时下热门,举例论证。要不怎么说,人家能拿到投资,能做老板呢。

十二点半左右,饥肠辘辘的同事们终于结束了“会”,下面进入激动人心的“午餐”环节。顷刻间茶水间挤满了人,看到怀里抱着三明治、沙拉、香蕉的可爱同事们从一边穿梭到另一半的时候,手里拿着一块披萨的我,可能显得有些另类。不过我也早已看惯了这些,因为好吃的有些,所以不抢先那一些,等下回来再拿可能就没了。其实看到他们那样,我也想拿,只是总觉得不好意思那么做。心想或许那就是“吃相不雅”吧。当然,想要吃相雅观的话,就只能吃别人挑剩下的,或者干脆没得吃,比如Gavin就是。

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

就在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在做了,相信结合两个方法可以帮医生诊断更多眼部疾病

千万富翁的财富之旅

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

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

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

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

投资房产
买域名
投资自己

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

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

Clinvar学习笔记

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学习笔记

最近在学习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
2
3
4
5
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新增用户配置

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

1
2
3
4
5
6
#增加一个名为"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
2
3
4
5
6
7
8
9
10
11
12
13
#将"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

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:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    #!/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.