0%

在安装R包的时候遇到报错,C++14 standard requested but CXX14 is not defined

查了很多办法,刚开始是根据https://github.com/stan-dev/rstan/issues/892修改.R下面的Makevars,

但是包另外一个错g++: error: unrecognized command line option ‘-std=c++14’

于是继续查到c++1y这个问题,但依然没有解决问题。

复盘了一下,感觉是gcc的问题,所以升级了最新的gcc

1
2
3
4
# 系统是CentOS
sudo yum install centos-release-scl
sudo yum install devtoolset-10
scl enable devtoolset-10 bash

但是装包的时候新版的gcc依然不能别识别,所以修改Makevars,最终用了如下的配置,重点是指定了新版的g++和c++的路径,这样问题就解决了

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
27
28
29
30
31
MAKEFLAGS = -j18

## C++ flags
CXX=g++
CXX11=g++
CXX14=/opt/rh/devtoolset-9/root/usr/bin/g++
CXX17=g++

CXXFLAGS=-O3 -march=native -Wno-ignored-attributes
CXX11FLAGS=-O3 -march=native -Wno-ignored-attributes
CXX14FLAGS=-O3 -march=native -Wno-ignored-attributes
CXX17FLAGS=-O3 -march=native -Wno-ignored-attributes

CXXPICFLAGS=-fPIC
CXX11PICFLAGS=-fPIC
CXX14PICFLAGS=-fPIC
CXX17PICFLAGS=-fPIC

CXX11STD=-std=c++11
CXX14STD=-std=c++14
CXX17STD=-std=c++17

## C flags
CC=/opt/rh/devtoolset-10/root/usr/bin/gcc
FLAGS=-O3 -march=native

## Fortran flags
FC=gfortran
F77=gfortran
FFLAGS=-O3 -march=native
FCFLAGS=-O3 -march=native

总结:
1,upgrade gcc

2, specify the absolute gcc and g++ path

如果可以的话,建议把整个系统的gcc都替换成新版的

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

下载数据的时候,偶然碰到了SRA-Explorer,感觉挺好用的,地址:https://sra-explorer.info/

这个页面本身非常小,见https://github.com/ewels/sra-explorer,利用的是SRA API。

检索好之后,选择你想下载的样本,点击Add ** to collection,然后点击右上角saved datasets,页面下方就可以原始的fastq的链接,用curl下载fastq的命令,用aspera下载fastq的命令,还有下载SRA的命令,以及样本的metadata。非常好用。

和SRA的Run Selector类似,https://www.ncbi.nlm.nih.gov/Traces/study/?

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

我们常用ROC曲线来衡量模型的预测能力,但很少关注模型的校准度calibration。

Calibration curve的横坐标是我们用模型预测的probability,比如我预测的是可能是肿瘤患者的概率,risk probability,纵坐标是真实的事件event的概率或者事件的proportion。

我们画图的时候,通常是点,需要进行平滑smooting,常用的是 loess方法。

如果我们预测的概率和真实的事件概率完全匹配,那么calibration plot则是一条斜率为1的直线(perfect calibration)。如果我们画出的先在perfect calibration上面,则表示underestimated,就是我们预测的风险或者概率低,真实情况却比较高;反之如果在perfect calibration在下,则是overestimated ,我们预测的风险比真实的风险要高。

根据吻合的程度,又分为poor calibration或者strong calibration,也可以做Hosmer–Lemeshow (HL) goodness-of-fit test。

参考:https://bmcmedicine.biomedcentral.com/articles/10.1186/s12916-019-1466-7

1
2
3
4
5
6
7
8
9
10
# rms包可以画calibration plot
library(rms)
val.prob(predicted.probability, TRUE.Label)

如果模型用rms建立的话,还可以用
rms::calibrate的方法

#我自己也参考别人写了一个函数,结果如下图,欢迎使用
library(loonR)
riskCalibrationPlot(TRUE.Label, predicted.probability)

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

特征选择是机器学习中的一个重要步骤,通过特征选择挑选出对预测起重要作用的变量,既可以减少数据的维度,也可以减少计算的消耗,同时也有助于我们对自己的数据的理解。有许多方法都可以应用到特征选择,比如大家常用的LASSO。我在用R做数据分析的时候,看到过这个帖子,进而了解了很多算法,所以对这个帖子进行了翻译,方便自己复习,也方便大家学习。

原文参考: https://www.machinelearningplus.com/machine-learning/feature-selection/

  1. Boruta
  2. Variable Importance from Machine Learning Algorithms
  3. Lasso Regression
  4. Step wise Forward and Backward Selection
  5. Relative Importance from Linear Regression
  6. Recursive Feature Elimination (RFE)
  7. Genetic Algorithm
  8. Simulated Annealing
  9. Information Value and Weights of Evidence
  10. DALEX Package

Introduction

真实的数据中,有些变量可能仅是噪声,并没有多少重要意义。

这类变量占用内存空间、消耗计算资源,我们最好去除这类变量,特别是在很大的数据集中。

有时候,我们有一个具有业务意义的变量,但不确定它是否确实有助于预测Y。还有一个事实:在一个机器学习算法中有用的特征(例如决策树) 可能其他算法中(例如回归模型)不被选用或者低估。

同时,有些变量单独预测Y的性能不好,但与其他预测变量/特征组合的情况下却非常显著。比如说有些变量与预测指标的相关性很低,但在其他变量参与的情况下,它可以帮助解释某些其他变量无法解释的模式/现象。

在这些情况下,很难决定包含还是去掉这些变量/特征。

这里讨论的策略可以解决这些问题,同时可以帮助理解对于一个模型而言,变量的重要性与否importance,以及对模型有多少贡献。

重要的一点是,我们最希望使用的变量是既具有业务意义同时也有重要性方面的指标。

我们这里导入Glaucoma 数据集,此数据集的目标是通过63个不同的生理测量指标来预测青光眼的与否。

1
2
3
4
5
6
# Load Packages and prepare dataset
library(TH.data)
library(caret)
data("GlaucomaM", package = "TH.data")
trainData <- GlaucomaM
head(trainData)

Glaucoma数据集

1. Boruta

Boruta 是一个基于随机森林对特征进行排名和筛选的算法。

Boruta 的优势是它可以明确的决定变量是否重要,并且帮助选择那些统计显著的变量。此外,可以通过调整p值和maxRuns来调整算法的严格性。

maxRuns是算法运行的次数,该参数值越大,会选择出更多的变量。默认是100.

在决定特征是否重要的过程中,一些特征可能会被标为Tentative。有时增加maxRuns或许可以解决特征的暂定行Tentativeness。

以TH.data的Glaucoma数据集为例子。

1
2
# install.packages('Boruta')
library(Boruta)

Boruta的方程用的公式同其他预测模型类似,响应变量response在左边,预测变量在右边。

如果预测变量处输入的是点,则意味着所有变量都会纳入评价中。

doTrace参数控制打印到终端的数目。该值越高,打印的log信息越多。为了节省空间,这里设置0,你可以设置1和2试试。

结果输出在boruta_output.

1
2
# Perform Boruta search
boruta_output <- Boruta(Class ~ ., data=na.omit(trainData), doTrace=0)

看看里面包含哪些内容

1
2
3
4
5
6
7
8
9
10
11
names(boruta_output)
1. ‘finalDecision’
2. ‘ImpHistory’
3. ‘pValue’
4. ‘maxRuns’
5. ‘light’
6. ‘mcAdj’
7. ‘timeTaken’
8. ‘roughfixed’
9. ‘call’
10. ‘impSource’
1
2
3
4
5
6
7
# 得到显著的或者潜在的变量
boruta_signif <- getSelectedAttributes(boruta_output, withTentative = TRUE)
print(boruta_signif)
[1] "as" "ean" "abrg" "abrs" "abrn" "abri" "hic" "mhcg" "mhcn" "mhci"
[11] "phcg" "phcn" "phci" "hvc" "vbss" "vbsn" "vbsi" "vasg" "vass" "vasi"
[21] "vbrg" "vbrs" "vbrn" "vbri" "varg" "vart" "vars" "varn" "vari" "mdn"
[31] "tmg" "tmt" "tms" "tmn" "tmi" "rnf" "mdic" "emd"

如果不确定潜在的变量是否保留,可以对boruta_output进行TentativeRoughFix

1
2
3
4
5
6
7
# Do a tentative rough fix
roughFixMod <- TentativeRoughFix(boruta_output)
boruta_signif <- getSelectedAttributes(roughFixMod)
print(boruta_signif)
[1] "abrg" "abrs" "abrn" "abri" "hic" "mhcg" "mhcn" "mhci" "phcg" "phcn"
[11] "phci" "hvc" "vbsn" "vbsi" "vasg" "vbrg" "vbrs" "vbrn" "vbri" "varg"
[21] "vart" "vars" "varn" "vari" "tmg" "tms" "tmi" "rnf" "mdic" "emd"

这样Boruta就代表我们来决定是否保留Tentative变量。查看这些变量的重要性分值。

1
2
3
4
# Variable Importance Scores
imps <- attStats(roughFixMod)
imps2 = imps[imps$decision != 'Rejected', c('meanImp', 'decision')]
head(imps2[order(-imps2$meanImp), ]) # descending sort
meanImp decision
varg 10.279747 Confirmed
vari 10.245936 Confirmed
tmi 9.067300 Confirmed
vars 8.690654 Confirmed
hic 8.324252 Confirmed
varn 7.327045 Confirmed

用画图的形式来展示变量的重要性

1
2
# Plot variable importance
plot(boruta_output, cex.axis=.7, las=2, xlab="", main="Variable Importance")

Variable Importance Boruta

这幅图展示了每个变量的重要性。

绿色表示的是筛选出的变量confirmed,红色则是需要提出的变量,蓝色不是真实的特征,待变的是ShadowMaxShadowMin,用来决定变量重要性与否。

2. Variable Importance from Machine Learning Algorithms

另外一种特征选择的方法是将各种ML算法最常用的变量视为最重要的变量。

根据机器学习算法学习X与Y之间关系的方式,不同的机器学习算法选出不同的变量(大多数是重叠的),但赋予的权重并不相同。

例如,在基于树的算法(如“ rpart”)中被证明有用的变量在基于回归的模型中可能没那么有用。 因此,不同算法没必要用相同的变量。

那么对于一个特定的机器学算法,如何找到变量的重要性?

  1. 利用caret 包中的train()一个特定的模型

  2. 使用 varImp() 来决定变量的重要性。

可能尝试多种算法,以了解各种算法之间的重要变量。

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
27
28
29
30
31
# 训练rpart模型,计算变量重要性
library(caret)
set.seed(100)
rPartMod <- train(Class ~ ., data=trainData, method="rpart")
rpartImp <- varImp(rPartMod)
print(rpartImp)
rpart variable importance

only 20 most important variables shown (out of 62)

Overall
varg 100.00
vari 93.19
vars 85.20
varn 76.86
tmi 72.31
vbss 0.00
eai 0.00
tmg 0.00
tmt 0.00
vbst 0.00
vasg 0.00
at 0.00
abrg 0.00
vbsg 0.00
eag 0.00
phcs 0.00
abrs 0.00
mdic 0.00
abrt 0.00
ean 0.00

rpart仅使用了63个功能中的5个,如果仔细观察,这5个变量位于boruta选择的前6个中。

让我们再做一件事:正则随机森林(Regularized Random Forest ,RRF)算法中的变量重要性。

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
27
28
29
30
31
# Train an RRF model and compute variable importance.
set.seed(100)
rrfMod <- train(Class ~ ., data=trainData, method="RRF")
rrfImp <- varImp(rrfMod, scale=F)
rrfImp
RRF variable importance

only 20 most important variables shown (out of 62)

Overall
varg 24.0013
vari 18.5349
vars 6.0483
tmi 3.8699
hic 3.3926
mhci 3.1856
mhcg 3.0383
mv 2.1570
hvc 2.1357
phci 1.8830
vasg 1.8570
tms 1.5705
phcn 1.4475
phct 1.4473
vass 1.3097
tmt 1.2485
phcg 1.1992
mdn 1.1737
tmg 1.0988
abrs 0.9537
plot(rrfImp, top = 20, main='Variable Importance')

正则随机森林中的变量重要性。

最前面的重要变量也于Boruta选择的吻合

其它的在train()中可以用的算法也可以用varImp计算变量的重要性:包括

ada, AdaBag, AdaBoost.M1, adaboost, bagEarth, bagEarthGCV, bagFDA, bagFDAGCV, bartMachine, blasso, BstLm, bstSm, C5.0, C5.0Cost, C5.0Rules, C5.0Tree, cforest, chaid, ctree, ctree2, cubist, deepboost, earth, enet, evtree, extraTrees, fda, gamboost, gbm_h2o, gbm, gcvEarth, glmnet_h2o, glmnet, glmStepAIC, J48, JRip, lars, lars2, lasso, LMT, LogitBoost, M5, M5Rules, msaenet, nodeHarvest, OneR, ordinalNet, ORFlog, ORFpls, ORFridge, ORFsvm, pam, parRF, PART, penalized, PenalizedLDA, qrf, ranger, Rborist, relaxo, rf, rFerns, rfRules, rotationForest, rotationForestCp, rpart, rpart1SE, rpart2, rpartCost, rpartScore, rqlasso, rqnc, RRF, RRFglobal, sdwd, smda, sparseLDA, spikeslab, wsrf, xgbLinear, xgbTree.

3. Lasso Regression

最小绝对收缩和选择算子(LASSO,Least Absolute Shrinkage and Selection Operator)回归是一种用L1范数惩罚的正则化方法。

从根本上来说,要增加权重(系数值)会带来成本。 它被称为L1正则化,增加的成本,与权重系数的绝对值成正比。

结果,在收缩系数的过程中,最终将某些不需要的特征的系数全部减小到零。 也就是说,它删除了不重要的变量。

LASSO回归也可以被视为有效的变量选择技术。

1
2
3
4
5
6
7
8
9
10
11
12
library(glmnet)
trainData <- read.csv('https://raw.githubusercontent.com/selva86/datasets/master/GlaucomaM.csv')

x <- as.matrix(trainData[,-63]) # all X vars
y <- as.double(as.matrix(ifelse(trainData[, 63]=='normal', 0, 1))) # Only Class

# Fit the LASSO model (Lasso: Alpha = 1)
set.seed(100)
cv.lasso <- cv.glmnet(x, y, family='binomial', alpha=1, parallel=TRUE, standardize=TRUE, type.measure='auc')

# Results
plot(cv.lasso)

LASSO变量重要性

X轴是log之后的lambda,当是2的时候,lambda的真实值是100。

图的最上面显示了模型使用了多少个变量,相对应的红点Y值则是在这些变量使用的情况下,模型可以达到多少的AUC。

可以看到两条垂直虚线,左边的第一个指向具有最小均方误差的lambda。 右边的一个表示在1个标准偏差内偏差最大的变量的数量。

最优的lambda值存储在cv.lasso $ lambda.min中。

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
27
28
29
30
31
# plot(cv.lasso$glmnet.fit, xvar="lambda", label=TRUE)
cat('Min Lambda: ', cv.lasso$lambda.min, '\n 1Sd Lambda: ', cv.lasso$lambda.1se)
df_coef <- round(as.matrix(coef(cv.lasso, s=cv.lasso$lambda.min)), 2)

# See all contributing variables
df_coef[df_coef[, 1] != 0, ]
Min Lambda: 0.01166507
1Sd Lambda: 0.2513163

Min Lambda: 0.01166507
1Sd Lambda: 0.2513163

(Intercept) 3.65
at -0.17
as -2.05
eat -0.53
mhci 6.22
phcs -0.83
phci 6.03
hvc -4.15
vass -23.72
vbrn -0.26
vars -25.86
mdt -2.34
mds 0.5
mdn 0.83
mdi 0.3
tmg 0.01
tms 3.02
tmi 2.65
mv 4.94

上面的输出显示了LASSO认为重要的变量。绝对值越高表示该变量越重要。

Read more »

Pacbio的工具更新实在是太快了,https://github.com/PacificBiosciences/pbbioconda

原来Isoseq的数据处理之后还叫high quality reads,现在叫HiFi reads,而且Pacbio的工具包有这么个倾向,把Isoseq分析过程中用的工具都囊括在内,比如后续的比对和collapse的工作。在tofu停了之后,cDNA_Cupcake也顶上来了。下面是几个名词解释,CLR, CCS, subreads, scrapes等,在处理原始的raw data时会碰到。

The contiguous sequence generated by the polymerase during sequencing is referred to as a “polymerase read” or a Continuous Long Read (CLR).

This CLR read may include sequence from adapters and multiple copies of inserts, because it traverses the circular template many times. The CLRs are processed to remove adapter sequences and to retain only the insert sequence, called “subreads”.

All other sequences sequenced from the CLR are called “scraps”.

Multiple copies of subreads generated from the single SMRTBell can then be collapsed to a single, high-quality sequence, called the “read of insert” ROI or Circular Consensus Sequence (CCS).

The scraps file includes all the sequence not in the HQ (high quality) region, and the adapter sequences. The old bax format saved everything with indexes into the usable data. The subreads.bam contains all the usable data.

关于scrapes和subreads的bam文件的区别,只用subreads即可
The scraps file includes all the sequence not in the HQ (high quality) region, and the adapter sequences. The old bax format saved everything with indexes into the usable data. The subreads.bam contains all the usable data.

.subreads.bam A file contain usable reads without adapter;The Sequel System outputs one subreads.bam file per collection/cell, which contains unaligned base calls from high-quality regions.

.scraps.bam A file contain sequence data outside of the High Quality region, rejected subreads,excised adapter and possible barcode sequences, as well as spike-in control sequences. (The basecaller marks regions of single molecule sequence activity as high-quality.)

Ref:
http://www.zxzyl.com/archives/1044

https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-020-03751-8

http://seqanswers.com/forums/showpost.php?p=202171&postcount=9

帅旸谈一谈在变异解读过程中用到的几个不太熟悉的预测指标:

z score

z score:这个指标指的是某个基因对missense的耐受程度,具体是指该基因所期望的missense数比上观察

到的missense数,如果z score>3.09,则认为该基因对missense不耐受,根据公式我们可以看出如果比值越大,则基因对missense越不耐受。利用z score可以在我们使用ACMG指南PP2的时候使用。

REVEL score

REVEL score:ClinGen SVI建议使用REVEL用来预测missense致病性。与其他常用missense致病性预测软件不同,REVEL整合了包括SIFT、PolyPhen、GERP++在内的13个软件的预测结果,对罕见变异的预测结果更加出色。当REVEL score>0.75,<0.15时分别使用ACMG指南PP3和BP4。

GERP++

GERP++ rejected substitutions” (RS) score:GERP++从基因进化速率角度预测位点保守性,具体是指该基因位点所期望的碱基替换次数减去观察到的碱基替换次数,可见分数越大,该位点保守性较强,当GERP++ RS score>6.8时,认为该位点保守。当分析一个不影响剪切的同义突变时,如果RS score<6.8,则可以使用ACMG指南BP7。

dbscSNV score

dbscSNV score:dbscSNV含有两个不同的算法,用来预测变异是否影响截切,一个是基于adaptive boostin,一个是基于Random Forest。当两种算法得分均小于0.6时,则认为不影响剪切。

参考

ClinGen及Zhang, J., Yao, Y., He, H., & Shen, J. (2020). Clinical Interpretation of Sequence Variants. Current Protocols in Human Genetics, 106(1), e98.

其他参数如pLI及 Haploinsufficiency score我们之前已经介绍过,请点击下文链接

http://www.zxzyl.com/archives/940

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

datapasta

https://github.com/MilesMcBain/datapasta/

还在手工的把excel的数据写成导到R里吗。不管横着还是竖着复制数据,datapasta可以自动、快速的把复制数据转成tibbles, data.frames, 或者 vectors格式。

更详细的参考https://cran.r-project.org/web/packages/datapasta/vignettes/how-to-datapasta.html

https://raw.githubusercontent.com/lorenzwalthert/some_raw_data/master/styler_0.1.gif

styler

https://github.com/r-lib/styler

还嫌自己写出的代码不够美吗,styler 可以把代码格式化成tidyverse规则的风格。

The goal of styler is to provide non-invasive pretty-printing of R source code while adhering to the tidyverse formatting rules. styler can be customized to format code according to other style guides too.

https://raw.githubusercontent.com/lorenzwalthert/some_raw_data/master/styler_0.1.gif

真香

links

https://github.com/MilesMcBain/datapasta/

https://github.com/r-lib/styler

https://github.com/daattali/addinslist

https://www.zhihu.com/question/398418315

新年快乐,21年的第一篇文章。

以前写过映射ENSEMBL ID 和 NCBI ID, http://www.zxzyl.com/archives/736

日常分析中,我们也会经常遇到其他的ID mapping的工作,这种工作不是基因ID转基因ID,而是转录本的ID转基因ID。

如果用的是refGene的注释,最简单了,直接用下面的命令即可

1
mysql --user=genome -N --host=genome-mysql.cse.ucsc.edu -A -D hg38 -e "select name,name2 from refGene"

不过我也经常通过解析gtf文件获得,因为gtf有转本的ID,也有基因的symbol或者ID,只要有gtf文件就可以提取。本着不造轮子的精神,我利用的是现成的R包

1
2
3
library(plyranges)
gr <- read_gff("/path/to/gtf/or/gff") %>% select(transcript_id, gene_id, gene_name)
gr <- unique(data.frame(gr))

我也在自己的包里面写了一个函数得到ensembl,refseq,hgnc,gene symbol的对应关系,biomaRt比较慢,可以把结果保存成文件

1
2
3
4
5
devtools::install_github("ProfessionalFarmer/loonR")
# 需要安装biomaRt和dplyr
mapping.table <- loonR::get_full_mapping_table()
# 保存mapping.table

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