分类目录归档:Default Category

cosine similarity

在SNV分析中,我们在算signature和样本mutation spectrum之间的相似性时,会用到cosine similarity。cosine similarity (distance)的公式,其实就是两个向量的夹角的cosine值,计算公式如下

它与欧式距离的差别如下图,cosθ就是similarity,而d则是欧氏距离Euclidean distance。

有些时候,距离也算作一种相似性,因为距离越远,说明两个样本越不相似。Euclidean distance和cosine similarity要根据情况来选择,最重要的是,是否要考虑weight or magnitude,参考下图。

在文本挖掘分析的时候,计算两个文本的相似性,我们可以统计每个词出现的次数,然后计算相似性(距离),因为文章有短有长,如果考虑单词出现的次数,那么字数多的文章一定与字数少的文章不一样(欧氏距离),所以如果我们不考虑这个量(magnitude)的时候,用cosine计算更加合适,结果也与欧氏距离不一样。

参考:

https://cmry.github.io/notes/euclidean-v-cosine#:~:text=While%20cosine%20looks%20at%20the,to%20actually%20measure%20the%20distance.

https://www.machinelearningplus.com/nlp/cosine-similarity/

基于DNA或RNA的NGS数据进行HLA分型

写这个原因呢,最近又要对样本的HLA分子进行分型,然后看到某公司的微信公众号讲的HLA的分型软件,全文讲了那么多,要么巨难用,要么下载不到,反正不如我自己正在用的这两个。另外一方面,没必要太纠结非常高的精度,除非你用得到。4-digital resolution,我觉得已经够了。

HLA分子

先回顾下百度百科对HLA的介绍(https://baike.baidu.com/item/HLA/9504270?fr=aladdin):

HLA(human leukocyte antigen ,人类白细胞抗原)是人类的主要组织相容性复合体(MHC)的表达产物,该系统是所知人体最复杂的多态系统。

HLA是具有高度多态性的同种异体抗原,其化学本质为一类糖蛋白,由一条α重链(被糖基化的)和一条β轻链非共价结合而成。其肽链的氨基端向外(约占整个分子的3/4),羧基端穿入细胞质,中间疏水部分在胞膜中。HLA按其分布和功能分为Ⅰ类抗原和Ⅱ类抗原。

HLA-I类分子:内源性抗原的递呈分子, HLA-Ⅱ类分子:外源性抗原的递呈分子

我想介绍的是seq2hla和optitype这两个软件

软件的大致原理,都是和HLA数据库中的数据进行比对,然后分型。我拿过7个样本的RNA数据用这两个软件进行比较,如下表,虽然个别分型不一致(后两个resolution不一致),但大部分结果还是相当一致的。

   
Sample   
   
   
   
A1   
   
A2   
   
B1   
   
B2   
   
C1   
   
C2   
1    
seq2HLA   
   
A*02:01   
   
A*24:18   
   
B*35:102   
   
B*35:102   
   
C*16:02   
   
C*04:01   
   
OptiType   
   
A*02:01   
   
A*24:02   
   
B*35:14   
   
B*35:14   
   
C*16:02   
   
C*04:01   
2    
seq2HLA   
   
A*29:01   
   
A*02:07   
   
B*46:01   
   
B*15:01   
   
C*04:01   
   
C*01:02   
   
OptiType   
   
A*29:01   
   
A*02:07   
   
B*46:01   
   
B*15:01   
   
C*04:01   
   
C*01:02   
3    
seq2HLA   
   
A*02:07   
   
A*11:01   
   
B*52:01   
   
B*15:01   
   
C*12:02   
   
C*03:03   
   
OptiType   
   
A*02:07   
   
A*11:01   
   
B*52:01   
   
B*15:01   
   
C*12:02   
   
C*03:03   
4    
seq2HLA   
   
A*02:07   
   
A*11:01   
   
B*46:01   
   
B*40:01   
   
C*01:02   
   
C*07:02   
   
OptiType   
   
A*02:07   
   
A*11:01   
   
B*46:01   
   
B*40:01   
   
C*01:02   
   
C*07:02   
5    
seq2HLA   
   
A*33:03   
   
A*11:01   
   
B*58:01   
   
B*15:02   
   
C*03:02   
   
C*08:01   
   
OptiType   
   
A*33:03   
   
A*11:01   
   
B*58:01   
   
B*15:02   
   
C*03:02   
   
C*08:01   
6    
seq2HLA   
   
A*02:07   
   
A*29:01   
   
B*46:01   
   
B*07:02   
   
C*01:02   
   
C*15:05   
   
OptiType   
   
A*02:07   
   
A*29:01   
   
B*46:01   
   
B*07:05   
   
C*01:02   
   
C*15:05   
7    
seq2HLA   
   
A*02:03   
   
A*11:01   
   
B*56:01   
   
B*56:01   
   
C*01:02   
   
C*01:02   
   
OptiType   
   
A*02:03   
   
A*11:01   
   
B*56:01   
   
B*56:01   
   
C*04:01   
   
C*01:02   

seq2hla

首先介绍的是seq2hla,这个其实一个命令就可以进行HLA的分型(很实用吧),但只基于RNA Seq的数据,我见有人用WES的数据来做,但不确定靠不靠谱。

官网: https://bitbucket.org/sebastian_boegel/seq2hla/wiki/Home

利用conda安装

conda create -n seq2hla
conda install -n seq2hla -c bioconda seq2hla

运行

seq2HLA -1  readfile1  -2  readfile2  -r  runname  [-p int] [-3 int]
# 提供R1和R2,现在已经支持gzip格式,-p是线程,-3是trim read末端的碱基数(因为read末端的碱基质量不高)

seq2hla的结果

下面的结果基于RNA-seq的数据得到的,在介绍optitype的时候,我用的是WES-seq的数据,可以看到虽然软件和数据类型不同,但两个软件的结果是一致的,综合前面表格中两种软件基于同一样本RNA-seq的结果,和两种软件基于同一样本不同数据类型的结果来看,都能说明两个软件都比较靠谱。

#Locus  Allele 1        Confidence      Allele 2        Confidence
A       A*24:02 0.001426754     A*11:01 0.003463472
B       B*38:02'        0.0002728549    B*46:01 0.0002107604
C       C*01:02 0.01282529      C*07:02 0.01516768
#Locus  Allele 1        Confidence      Allele 2        Confidence
DQA     no      NA      no      NA
DQB     DQB1*05:02      NA      DQB1*05:02      NA
DRB     DRB1*11:52      0.297521        DRB1*04:05'     0.42202

optitype

另外一个软件是optitype,这个软件呢,稍微多了几步,但绝对在让人接受的范围内,而且这个软件支持DNA数据和RNA数据。

官网:https://github.com/FRED-2/OptiType

利用conda安装

这个时候真的要安利下conad,有了conda之后,各种软件包的安装和管理,变得非常轻松。拿optitype来说,我第一次安装的时候,自己下载,编译,修改部分文件,还要自己安装razers,而有了conda之后,只需一个命令

conda create -n optitype
conda install -n optitype -c bioconda optitype

optitype需要先利用razers3先过滤一下reads,把和HLA相关的序列提出来,然后在运行optitype进行分型,从而加快分析速度。

第一步,Razers3

razers3 -i 95 -m 1 -dr 0 --thread-count 10 -o fished_1.bam /path/to/hla_reference_dna.fasta sample_1.fastq

我是用conda装的,但在conda的env目录下没找到hla的fasta文件,不过不用着急,官网上有,可以从https://github.com/FRED-2/OptiType/tree/master/data下载对应的HLA的DNA或RNA的数据。

如果是paired的read,则需要分别跑一次,得到fished_1.bam和fished_2.bam文件。

另外官网提示Razers3会把所有的fastq读到内容,所以计算服务器的内存最好大点。

第二步,将bam文件转成fastq

samtools bam2fq fished_1.bam > fished_1.fastq
samtools bam2fq fished_2.bam > fished_2.fastq

第三步,运行optitype

OptiTypePipeline.py -i fished_1.fastq fished_2.fastq --dna -v -o /path/to/outdir -p sample.optitype.dna

–dna需要和前面Razers3用的HLA的序列类型一致,如果前面用的RNA的数据,这里也是–rna,当然用RNA还是DNA取决于测序数据。

optitype的结果

这个结果和seq2hla基于RNA-seq的数据得到的HLA亚型是一致的。

        A1      A2      B1      B2      C1      C2      Reads   Objective
0       A*11:01 A*24:02 B*38:02 B*46:01 C*01:02 C*07:02 200.0   191.00000000000003

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

Agilent全外芯片的目标区域下载

安捷伦的全外芯片捕获的目标区域的bed文件是可以下载的,下载网址

https://earray.chem.agilent.com/suredesign

这个网站我没保存过,毕竟用的频率不高,但每次想用的时候还要搜一下这个网站(Σ( ° △ °|||)︴)

注册登录之后,找到Find Design -> SureSelect DNA -> Agilent Catalog Designs,进一步筛选之后,可以下载对应的文件。

download-agilent-bed-target-region

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

立个走完麦理浩径的flag

麦理浩径一共十段,长达100公里,前几段风景比较优美,我只走过二段的一半,我想把十段都走完,先把flag立下(3年实现)。

2019-06-23 走了二段的一部分,景色真的非常好,山径、海岸、溪流、沙滩。又乘坐快艇去了西贡,享受海鲜大餐。

2020-06-20 走完了一段和剩余的二段的一部分。

一段和部分二段图片

起点
万宜水库
部分二段

二段图片

继续阅读

分类模型的性能评估

最常用的就是灵敏度和特异性,不过还有其他的,比如阴性预测值(negative predictive value, NPV)。


通常,先画一个ROC曲线,计算曲线下面积。ROC上的每个点是特定阈值下,分类的sensitivity和specificity,没多点连起来组成ROC,曲线下面积就是AUC。面积越大越好,如果AUC是1,说明模型能够完全区分要预测的类别。

如果不是1,就要考虑阈值取哪里比较好,这里就涉及到Youden index。Youden index 其实就是为了找到使得sensitivity和specificity之和最大max(sensitivities+specificities)的阈值。

另外就是考虑其他指标来评估分类模型的性能:specificity, sensitivity, accuracy, npv, ppv, precision, recall, tpr, fpr, tnr, fnr, fdr。这些指标可谓琳琅满目,不过这之间有重复的,如下,都是基于tn(真阴), tp(真阳), fn(假阴), fp(假阳)的个数进行计算。

 

预测

P

N

实际

P

TP

FN

N

FP

TN

因为经常用到,就罗列了一下。

具体描述 公式 别名
tn True negative count真阴数
tp True positive count真阳数
fn False negative count假阴数
fp False positive count假阳数
specificity Specificity特异度 tn / (tn + fp) tnr
sensitivity Sensitivity灵敏度 tp / (tp + fn) recall, tpr
accuracy Accuracy正确率 (tp + tn) / N
npv Negative Predictive Value阴性预测值 tn / (tn + fn)
ppv Positive Predictive Value阳性预测值 tp / (tp + fp) precision
precision Precision精准率 tp / (tp + fp) ppv
recall Recall正确率 tp / (tp + fn) sensitivity, tpr
tpr True Positive Rate真阳性率 tp / (tp + fn) sensitivity, recall
fpr False Positive Rate假阳性率 fp / (tn + fp) 1-specificity
tnr True Negative Rate真阴性率 tn / (tn + fp) specificity
fnr False Negative Rate假阴性率 fn / (tp + fn) 1-sensitivity
fdr False Discovery Rate伪发现率 fp / (tp + fp) 1-ppv