分类目录归档:Copyrighted

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

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

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: 帅旸
#####################################################################

推荐两个Rstudio Addins

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

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.

真香

links

https://github.com/MilesMcBain/datapasta/
https://github.com/r-lib/styler
https://github.com/daattali/addinslist
https://www.zhihu.com/question/398418315

最后推荐两个可以将ggplot导出成ppt的R包

1,export的graph2ppt函数

https://github.com/tomwenseleers/export

export虽然从CRAN下架了,但依然可以通过github的库来安装,devtools::install_github(“tomwenseleers/export”)

2,eoffice的topptx函数

https://github.com/guokai8/eoffice

3, officer

https://github.com/davidgohel/officer

用法可以参考:https://www.brodrigues.co/blog/2018-10-05-ggplot2_purrr_officer/

Map NM ID to Gene Symbol

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

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

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

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

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包

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

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

Parse gtf

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.

/wp/f4w/2020/2020-10-02-gtf-parser.gif

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.

>>> 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")

# Example gtf
Here is an simple example of gtf file. You can use to test. A subset from refSeq.hg38.gtf.

继续阅读

整合多组学数据的通路富集分析-ActivePathways

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

/wp/f4w/2020/2020-09-12-ActivePathways-1.png

方法

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的图,来分析更广泛的生物学主题。点为通路,边表示有共有基因。

/wp/f4w/2020/2020-09-12-ActivePathways-2.png

例子

自然发到了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的强大。

/wp/f4w/2020/2020-09-12-ActivePathways-3.png

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

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

参考

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

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