VAF (variant allele frequency) vs CCF (cancer cell fraction)

图片来源于网络MesKit

VAF - variant allele frequency

Variant allele fraction or frequency (VAF): the fraction of mutated reads for a given variant, which is a readout for the proportion of DNA mutated in the sequenced tissue.

测序时特定位点突变的reads数比上总的reads数,可以从VCF中获得。

CCF - cancer cell fraction

Cancer cell fraction (CCF): the fraction of cancer cells from the sequenced sample carrying a set of SNVs.

携带突变的癌细胞比例,可以通过pyclone(https://github.com/Roth-Lab/pyclone-vi)或sciclone(https://github.com/genome/sciclone)计算。

$$ CCF = VAF *\frac{1}{p}[pCN_t + CN_n(1-p)] $$

VAF: corresponds to the variant allele frequency at the mutated base

p: the tumor purity肿瘤纯度

CNt: the tumor locus specific copy number所在位置的拷贝数

CNn: the normal locus specific copy number (CNn was assumed to be 2 for autosomal chromosomes)正常样本的拷贝数

时间依赖的ROC

三种时间依赖的ROC

诊断模型的ROC是大家最熟悉的,一组二分类的真实标签,一组风险分值,不同的cutoff下有不同的灵敏度和特异性,就能画出ROC曲线。

生存分析一般建立Cox模型,根据Cox模型也会有一组风险分值,生存结局也是一个二分类的标签,但病例多了时间的信息。

三种不同的定义来估计删失事件的时间依赖的敏感性和特异性,即(1)cumulative/dynamic累积/动态(C/D),(2)incident/dynamic事件/动态(I/D)和(3) incident/static事件/静态 (I/S)。不同的定义下,灵敏度和特异性的计算不一样。其中C/D比I/D和I/S更具有临床相关性,在临床中普遍使用。

t: 目标时间,t* 固定的随访时间,c 阈值,A-F为研究中的个体病例,其中ABC的指标高于c,DEF个体的指标小于c,实心圆表示发生事件的个体,空心圆表示Censored个体。

图a,C/D、I/D、I/S中的相对于基线时间点的实验和对照个体说明,图b,I/S(纵向)的说明

Cumulative sensitivity and dynamic specificity (C/D)

$$ \begin{array}{l} S{e}^C\left( c, t\right)= P\left({X}_i> c\Big|{T}_i\le t\right)\\ {} S{p}^D\left( c, t\right)= P\left({X}_i\le c\Big|{T}_i> t\right)\\ {} AU{C}^{C, D}(t)= P\left({X}_i>{X}_j\Big|{T}_i\le t,{T}_j> t\right), i\ne j.\end{array} $$

cumulative/dynamic(C/D)中cumulative是指Cumulative sensitivity,dynamic是指dynamic specificity。C/D是用的最广泛的ROC模型。

灵敏度:生存时间小于t的人群之中,Xi大于阈值c的人群(A和B个体)所占总体(A、B和E)的比例

特异度:生存时间大于t的人群之中,Xi小于等于阈值c的人群(D和F)所占总体(C,D和F)的比例

用二分类模型来类比,时间ROC,在选定特定时间点(比如1年、3年或者5年)时,不同的阈值c可以将队列人群分成高于c的一组(高风险,认为发生了死亡/事件)和小于等于c(低风险,认为没有发生死亡/事件),但这个时间点有真实的死亡或者事件的标签,于是可以计算在特定时间点特定c时的灵敏度和特异性,进行画ROC计算AUC。所以文献中经常看到1yr, 3yr, 5yr ROC的图,有了ROC的图,其实也可以画DCA进行DCA分析。

通过snakemake向slurm提交任务

集群任务调度,最初读研究生的时候,接触的是实验室用的condor(https://research.cs.wisc.edu/htcondor/),目前还在维护。在脚本里面设定好请求的计算资源,交给master节点即可,调度系统会自动分配和管理任务。这也是我对集群管理中的任务调度的初始了解。后来不管是工作还是在实验室,遇到最多的是SGE(sun grid engine)。再后来SGE被Oracle抛弃,又接触了Slurm和Torque。不管何种系统,通过集群调度来实现分析的分布式计算,而不用关系具体的任务分配,极大的提高了集群的利用率和分析效率。

关于流程管理,如果只有一两个样本的情况下,我都是直接把pipeline放到sh脚本上直接跑,懒得用流程管理。但如果样本很多,成百上千个样本的时候,最好用流程管理,来追踪大规模任务的分析状态,是否报错,是否成功结束,以免手工check造成遗漏。我有时候会在命令之后加&& touch “Done"或者判断$0的状态,来确保程序正确执行,但这样做确实很繁琐。后来snakemake开始流行,又是和python结合,易读性很好,就尝试开始用snakemake。一个Snakefile文件,可以搞定N多样本,还能监控分析的进度。

正好最近需要在一个slurm调度管理的集群上进行分析,而流程本身就封装在Snakefile文件中,如果体验了snakemake和slurm合体,体验非常好,颠覆了我对向集群投递任务的繁琐印象。要是放在以前,我要分析1000个样本,我可能要生成1000个script,需要专门写一个生成job script的script,然后再投递。我也知道可以通过传参进行批量投递,但体验非常不好。所以还是感慨技术的进步,也推荐snakemake和slurm一起用。本文主要提一下如何微调下资源请求的命令。

假设我的Snakemake文件有1000个job,我需要向集群提交任务,如果和slurm或者SGE配合使用,需要用到–cluster和–jobs选项。

1
2
3
snakemake --cluster "sbatch -N 1 --ntasks=1  --cpus-per-task=10" --jobs 10
# --cluster 后面跟的是提交命令,如果是SGE的话,就是qsub
# --jobs是最大同时投递的任务数目

很多文档中提到要在Snakefile中的每个rule中设置resources,比如resources: mem_mb = 40000,但我发现这样提交的任务,是获取不了具体的请求资源。这个时候用scontrol show job jobid查看请求的资源,可以看到可以用10个cpu,但请求的内存依然是默认的。

查了很多教程,都是通过配置yaml或者json文件,我觉得这样很繁琐,而且每个rule的请求资源是不一样的,通过配置文件,相当于多了一个文件,多了很多工作。看到有人说–cluster中可以用wildcards,就打开了我的思路。比如在rule设置了相对应的内存和时间,可以通过wildcards来调用,不同的rule的任务,提交时就对应不同的资源请求。实现起来是这样的

1
snakemake --cluster "sbatch -N 1--mem={resources.mem_mb} --ntasks=1  --cpus-per-task=10 " --jobs 10

这个时候,可以进一步把snakefile里面的threads信息利用起来

1
snakemake --cluster "sbatch -N 1--mem={resources.mem_mb} --ntasks=1  --cpus-per-task={threads} " --jobs 10

这个时候scontrol show job jobid显示请求资源是我们想要的,但squeque的job名是snakemake.job,分不清现在正在运行的任务对应哪个rule,可以这么修改。

1
snakemake --cluster "sbatch -N 1--mem={resources.mem_mb} --ntasks=1  --cpus-per-task={threads} --job-name={rule}" --jobs 10

那么,进一步让输出日志和错误日志可读,可以这么修改

1
snakemake --cluster "sbatch -N 1--mem={resources.mem_mb} --ntasks=1  --cpus-per-task={threads} --job-name={rule} --output={rule}.%j.out --error={rule}.%j.err" --jobs 10

如果想进一步修改job名的可读性和日志文件的可读性,可以把rule里面的wildcards传进来。比如我的每个rule中都匹配了sample,那在传给 job-name和output, errror的时候,把{wildcards.sample}传过来就行。如果我不想把任务投递给一个node节点,可以通过exclude(比如–exclude=gpu1)来指定。

真空采血管的颜色分类

是的,没看错,我开始研究采血管的颜色了。

不同颜色的头盖和标签,表示不同的添加剂种类和试验用途。举个例子,提取血清和血浆,要分别用促凝管和抗凝管。又比如促凝管有两种,一个是橘色的,一种是含有分离胶的黄色管,抗凝管的成分有添加柠檬酸钠的浅蓝色管,也有EDTA紫色管。具体的添加剂和使用范围如下图。

删除tmp文件

服务器的tmp文件夹满了,导致不能正常登陆和运行软件,需要清理tmp文件夹下面的临时文件。Linux的机制是在重启的时候清理tmp文件夹,或者一段时间才删除(比如下面我们看到的服务器默认是10天),find命令和tmpwacth命令可以在不重启和不影响用户的情况下快速删除目标文件。

Find命令

1
2
# ctime即change time文件状态改变时间为超过一天的
find /tmp -ctime +1 exec rm -rf {} \;

tmpwatch

1
2
3
4
5
yum install tmpwatch -y
# 超过一天
tmpwatch 1d /tmp
# 超过32小时
tmpwatch -afv 32 /tmp/

定时清理

CentOS 6

看到说在CentOS下的/etc/cron.daily/tmpwatch可以设置自动清理时间,但我没找到,-_-||,例如下面设置的是/var/tmp下30d,也就是30天,/tmp下10天,但服务器/usr/sbin下面并没有tmpwatch,如果大家的服务器上有这个,可以通过这个设置定期清理一定时间的文件。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
#! /bin/sh 
flags=-umc 
/usr/sbin/tmpwatch "$flags" -x /tmp/.X11-unix -x /tmp/.XIM-unix \ 
        -x /tmp/.font-unix -x /tmp/.ICE-unix -x /tmp/.Test-unix \ 
        -X ‘/tmp/hsperfdata_*’ 10d /tmp 
/usr/sbin/tmpwatch "$flags" 30d /var/tmp 
for d in /var/{cache/man,catman}/{cat?,X11R6/cat?,local/cat?}; do 
    if [ -d "$d" ]; then 
        /usr/sbin/tmpwatch "$flags" -f 30d "$d" 
    fi 
done