作者归档:zzx

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

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

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
#####################################################################

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 500*2*0.5*0.5 = 250 (=   observed) 100*2*0.65*0.35 = 45.5 (observed has deficit of 15.5) 1,000*2*0.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 2*0.5*0.5=0.5 (Hexp 1) 2*0.65*3.5=0.455 (Hexp 2) 2*0.35*0.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.5*1000+0.65*200+0.35*2000)/3200=0.4156
q-bar (0.5*1000+0.35*200+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.5*500+0.3*100+0.5*1000)/1600=0.4875
HS (expected) (0.5*500+0.455*100+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 2)
     2) 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).
     3) 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),
     4) No excess or deficiency of   heterozygotes over the total population (FIT    is nearly zero).

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.

继续阅读