0%

我们常用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
# rms包可以画calibration plot
library(rms)
val.prob(predicted.probability, TRUE.Label)

#我自己也参考别人写了一个函数,结果如下图,欢迎使用
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
#####################################################################

Fixation index (FST)

整理来源 http://www.uwyo.edu/dbmcd/popecol/maylects/fst.html,把这个课程的步骤用表格来表示

Subpopulation 1 Subpopulation 2 Subpopulation 3 Total
Genotype AA 125 50 100
Aa 250 30 500
aa 125 20 400
Number of individual 500 100 1000 1600
Number of alleles 1000 200 2000 3200
Step 1. Calculate the gene (allele) frequencies
Observed allele frequency A (p) (125*2+250)/1000=0.5 (50*2+30)/200=0.65 (2*100+500)/2000=0.35
a (q) 0.5 0.35 0.65
Step 2. Calculate the expected genotypic counts under Hardy-Weinberg Equilibrium, and then calculate the excess or deficiency of homozygotes in each subpopulation. Summary of homozygote deficiency or excess relative to HWE: Pop. 1. Observed = Expected: perfect fit Pop. 2. Excess of 15.5 homozygotes: some inbreeding Pop. 3. Deficiency of 45 homozygotes: outbred or experiencing a Wahlund effect (isolate breaking).
Expected allele frequency AA 500*0.5^2 = 125 (= observed) 100*0.65^2 = 42.25 (observed has excess of 7.75) 1,000*0.35^2 = 122.5 (observed has deficiency of 22.5)
Aa 50020.5*0.5 = 250 (= observed) 10020.65*0.35 = 45.5 (observed has deficit of 15.5) 1,00020.65*0.35 = 455 (observed has excess of 45)
aa 500*0.5^2 = 125 (= observed) 100*0.35^2 = 12.25 (observed has excess of 7.75) 1,000*0.35^2 = 422.5 (observed has deficiency of 22.5)
Step 3. Calculate the local observed heterozygosities of each subpopulation (we will call them Hobs s, where the s subscript refers to the sth of n populations – 3 in this example).
Local observed heterozygosities 250/500 = 0.5 (Hobs 1) 30/100 = 0.3 (Hobs 2) 500/1000 = 0.5(Hobs 3)
Step 4. Calculate the local expected heterozygosity, or gene diversity, of each subpopulation Hexp = 2pq
Local expected heterozygosity 20.50.5=0.5 (Hexp 1) 20.653.5=0.455 (Hexp 2) 20.350.65=0.455 (Hexp 2)
Step 5. Calculate the local inbreeding coefficient of each subpopulation F = (Hexps -Hobs)/Hexp [positive F means fewer heterozygotes than expected indicates inbreeding] [negative F means more heterozygotes than expected means excess outbreeding]
F1=(0.5—0.5)/0.5=0 F2=(0.455—0.3)/0.455=0.341 F3=(0.455—0.5)/0.455=-0.099
Step 6. and 7. Calculate p-bar (p-bar, the frequency of allele A) over the total population. Calculate q-bar (q-bar, the frequency of allele a) over the total population. Check: p-bar + q-bar = 1.0
the frequency of allele over the total population p-bar (0.51000+0.65200+0.35*2000)/3200=0.4156
q-bar (0.51000+0.35200+0.65*2000)/3200=0.5844
Step 8. Calculate the global heterozygosity indices (over Individuals, Subpopulations and Total population) HI based on observed heterozygosities in individuals in subpopulations HS based on expected heterozygosities in subpopulations HT based on expected heterozygosities for overall total population
HI (observed) (0.5500+0.3100+0.5*1000)/1600=0.4875
HS (expected) (0.5500+0.455100+0.455*1000)/1600=0.4691
HT (in overall total population) 2*p-bar *q-bar = 2 * 0.4156 * 0.5844 = 0.4858
Step 9. Calculate the global F-statistics Compare and contrast the global FISbelow with the “local inbreeding coefficient” Fs of Step 5. Here we are using a weighted average of the individual heterozygosities over all the subpopulations. Both FIS and Fs are, however, based on the observed heterozygosities, whereas FST and FIT are based on expected heterozygosities.
FIS (Hs-Hi)/Hs=(0.4691-0.4875)/0.4691=-0.0393
FST (Ht-Hs)/Ht=(0.4858-0.4691)/0.4858=-0.0344
FIT (Ht-Hi)/Ht=(0.4858-0.4875)/0.4858=-0.0036
Step 10 conclusions

Finally, draw some conclusions about the genetic structure of the population and its subpopulations. 1) One of the possible HWE conclusions we could make: Pop. 1 is consistent with HWE (results of Step

  1. Two of the possible “local inbreeding” conclusions we could make from Step 5:

Pop. 2 is inbred (results of Step 5), and

Pop. 3 may have disassortative mating or be experiencing a Wahlund effect (more heterozygotes than expected).

  1. Conclusion concerning overall degree of genetic differentiation (FST)

Subdivision of populations, possibly due to genetic drift, accounts for approx. 3.4% of the total genetic variation (result of Eqn FST.8 FST calculation in Step 9),

  1. No excess or deficiency of heterozygotes over the total population (FIT is nearly zero).

I always use gtf file and retrieve gene information. There isn’t a highly flexible tool to solve my demand. I modified the code from “https://github.com/Jverma/GFF-Parser”, thanks Jverma. This tool will be easier to use.

Usage

Basically, there are three parameters.

id: either transcript id or gene id.

attType: attribute defined in gtf file. E.g. feature (column 3), gene_name, exon_number, transcript_id in column 9

attValue: the attribute value you want to search for.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
>>> import sys
>>> from gtfParser import gtfParser

>>> gtf = gtfParser("example.gtf")

>>> # Get all exons in CDK7
>>> gtf.getRecordsByID("CDK7", "feature", "exon")

>>> # Get all features of transcript_id defined as "NM_001324069" in gene "CDK7"
>>> gtf.getRecordsByID("CDK7", "transcript_id", "NM_001324069")

>>> # Get start codon where feature was defined as "start_codon" in transcript "NM_001324069"
>>> gtf.getRecordsByID("NM_001324069", "feature", "start_codon")

>>> # Get a exon where its id is "NM_001324078.1" in "NM_001324078" transcript
>>> gtf.getRecordsByID("NM_001324078", "exon_id", "NM_001324078.1")
Read more »

我是在这篇文章(Integrative pathway enrichment analysis of multivariate omics data
)中遇到的合并多个p-value的操作。这篇文章是今年发表在NC上。所有的组学或者大规模的数据分析,都需要探索数据背后相关的生物学功能,所以通路富集分析非常普遍。通常的做法是基于单一组学、单一数据集的数据进行分析,随着生物学数据的爆发,大规模多组学数据变得普遍,这篇文章介绍了基于整合的多组学或多数据集的数据进行通路分析的工具ActivePathways。

方法

ActivePathways的方法,如下图:

(a) 需要的输入文件

(1) 基于多组学数据集的基因P-value,传统的富集分析是单组学,只有一列,现在是多组学,对应多列P-value
(2) 基因集,这个和其他的通路富集分析一样,用来表示生物学过程和通路

(b)

(1) 用Brown method合并基因的P-values,并且排序,用一个宽松的阈值来过滤检阳性的基因。
(2) 对每个通路,用排序的基因(从第一个开始从少到多作为sub-list)进行超几何检验,并找到最优的sub-list长度。
(3) 基因单一组学的数据进行富集分析,找到支持每个通路的证据。

(c) ActivePathways 提供整合之后的富集分析结果,相关的Brown P-value,支持通路的证据。还可以在Cytoscape中画Enrichmentmap的图,来分析更广泛的生物学主题。点为通路,边表示有共有基因。

例子

自然发到了NC这种水平的期刊上,出了豪华的团队外,ActivePathways一定有很厉害的功能才行。这篇文章提了好几个case study,我挑了一个,稍微讲一下,如下图。

昨天分析了乳腺癌中与预后相关的通路和过程,其中整合了三种数据集,TC(Tumor cell mRNA),TAC(Tumor adjacent cell mRNA)和CNA(拷贝数变异),这里也挺有意思,把TAC也纳入到P-values的矩阵中,并不是三个组学。

(a) Enrichment map,其中蓝色的表示仅由整合数据发现的通路,比如凋亡等,显示出其强大的功能。
(b) 乳腺癌Basel和HER2亚型中,与预后相关的免疫基因的Hazard Ratio(HR),显示出两种亚型间不同的免疫pattern。
(c) HER2亚型中,与预后相关的基因(凋亡过程负调控)在每个数据集上的P-value和整合后的Brown P-value。其中DUSP1在三个数据集中都非常显著,如果作者聚焦到了这个基因。
(d) 这个基因在肿瘤细胞mRNA,癌旁细胞mRNA和拷贝数变异三种数据集中,分成high和low两组,做生存分析,可以发现DUSP1低表达,预后显著好于高表达,以此证明ActivePathway的强大。

ActivePaths的下载地址:https://github.com/reimandlab/ActivePathways。为一个R包,整理好P-values的数据框之后,一步命令即可分析,此外结果还可以在Cytoscape上用Enrichment map展示。

1
2
# scores是P-values的数据框,GMT是基因集
ActivePathways(scores, fname_GMT)

参考

https://www.nature.com/articles/s41467-019-13983-9

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