1. 
dNdScv R包是一组最大似然 dN/dS 方法，旨在量化癌症和体细胞进化中的选择。该包包含量化错义、无义和必需剪接突变的 dN/dS 比率的函数，可在单个基因、基因组或全外显子组水平上进行。包中的dndscv函数旨在检测全外显子组/基因组或靶向测序研究中从几个样本到数千个样本的数据集上的癌症驱动基因（即癌症中受到正向选择的基因）。
虽然该软件包最初是为癌症基因组研究而设计的，但它也可用于量化其他重测序研究中的选择，例如 SNP 分析、细菌中的突变积累研究或使用来自人类三重奏的数据发现导致发育障碍的突变。
通过结合局部信息（基因中的同义突变）和全局信息（基因间突变率的变化，利用表观基因组协变量），并控制基因的序列组成和突变特征，可以估计每个基因的背景突变率。与传统的 dN/dS 实现不同，dNdScv使用三核苷酸上下文相关替换矩阵来避免影响 dN/dS 的常见突变偏差.
1.1 安装并加载所需的R包
# 安装并加载所需的R包
# install.packages("devtools")
# library(devtools)
# install_github("im3sanger/dndscv")
library(vcfR)
library(dndscv)
library(IRanges)
library(dplyr)
# 下载所需的参考文件(dndscv 默认hg19)
$ wget https://github.com/im3sanger/dndscv_data/blob/master/data/RefCDS_human_GRCh38.p12_dNdScv.0.1.0.rda
$ wget https://github.com/im3sanger/dndscv_data/blob/master/data/RefCDS_human_GRCh38_GencodeV18_recommended.rda
1.2 VCF文件处理
vcf <- read.vcfR("/data/shumin/software/snpEff/Somatic/CC56tissueA_snpEff_annotated.vcf")
# 提取位置信息
pos_info <- as.data.frame(getFIX(vcf))
# 提取INFO字段中的ANN信息
ann_info <- extract.info(vcf, element = "ANN")
# 提取必要的列
variants_1 <- data.frame(
chr = pos_info$CHROM,
pos = pos_info$POS,
ref = pos_info$REF,
mut = pos_info$ALT
)
# 解析注释信息条目
parsed_annotations <- lapply(ann_info, function(ann_info) {
# 根据 | 分割注释信息
fields <- unlist(strsplit(ann_info, "|", fixed = TRUE))
# 返回需要的字段，例如基因名和变异类型
gene_name <- fields[4] # 基因名在第四个位置
variant_type <- fields[2] # 变异类型在第二个位置
impact <- fields[3] # 影响程度通常在第三个位置
return(list(Gene = gene_name, Variant_Type = variant_type, Impact = impact))
})
# 将列表转换为数据框
variants_2 <- do.call(rbind, parsed_annotations)
# 合并数据
variants <- cbind(variants_1, variants_2)
# 将属性为list的列转换为字符串
variants$Gene <- sapply(variants$Gene, function(x) paste(x, collapse = ""))
variants$Variant_Type <- sapply(variants$Variant_Type, function(x) paste(x, collapse = ", "))
variants$Impact <- sapply(variants$Impact, function(x) paste(x, collapse = ""))
# 处理缺失值：这里假设直接删除含有NA值的行
variants <- na.omit(variants)
# 确保每列的数据类型正确
variants$chr <- as.character(variants$chr) # 染色体编号
variants$pos <- as.integer(variants$pos) # 位置，转换为整数类型
variants$ref <- as.character(variants$ref) # 参考碱基
variants$mut <- as.character(variants$mut) # 突变碱基
variants$SampleID <- "CC56tissue" # 添加样本名一列
# 移除 "chr" 前缀
variants$chr <- gsub("^chr", "", variants$chr)
# 过滤出单碱基突变（SNPs）
variants <- variants %>%
filter(nchar(ref) == 1 & nchar(mut) == 1)
str(variants)
## 'data.frame': 1624817 obs. of 8 variables:
## $ chr : chr "1" "1" "1" "1" ...
## $ pos : int 13273 14653 14677 14717 16495 16875 17365 17538 17588 17614 ...
## $ ref : chr "G" "C" "G" "G" ...
## $ mut : chr "C" "T" "A" "A" ...
## $ Gene : chr "WASH7P" "DDX11L1" "DDX11L1" "DDX11L1" ...
## $ Variant_Type: chr "downstream_gene_variant" "downstream_gene_variant" "downstream_gene_variant" "downstream_gene_variant" ...
## $ Impact : chr "MODIFIER" "MODIFIER" "MODIFIER" "MODIFIER" ...
## $ SampleID : chr "CC56tissue" "CC56tissue" "CC56tissue" "CC56tissue" ...
# 筛选影响程度“HIGH"的数据进行后续分析
variants_sel <- subset(x=variants, variants$Impact==c('HIGH')) %>% .[,c(8,1:4)]
1.3 运行dndscv
dndsout <- dndscv(variants_sel, refdb="/usr/local/lib/R/library/dndscv/data/RefCDS_human_GRCh38_GencodeV18_recommended.rda", cv=NULL)
1）全局结果（$globaldnds）：这个表包含了所有基因全局的dN/dS值，用于评估正选择、负选择或中性进化
dndsout[["globaldnds"]]
## name mle cilow cihigh
## wmis wmis 0.6447645 0.3704856 1.122098
## wnon wnon 32.2041125 18.5571570 55.887055
## wspl wspl 17.5535865 9.1980014 33.499495
## wtru wtru 26.4194841 15.9302538 43.815318
## wall wall 3.4408418 2.1853716 5.417565
? name: 选择模式的名称（如wmis、wnon、wspl、wtru、wall）。
? mle: 最大似然估计（maximum likelihood estimate）的dN/dS值。
? cilow: dN/dS值的下限95%置信区间。
? cihigh: dN/dS值的上限95%置信区间
2）基因结果（$sel_cv）：这个表显示了每个基因的选择信号，用于确定哪些基因可能在癌症中起作用
dndsout[["sel_cv"]]
## gene_name n_syn n_mis n_non n_spl wmis_cv wnon_cv wspl_cv pmis_cv ptrunc_cv pallsubs_cv qmis_cv qtrunc_cv qallsubs_cv
## 15213 SLC9B1 0 0 3 0 0 10531.6440 10531.6440 0.9443649 2.092424e-10 1.611936e-09 0.9855767 4.020802e-06 3.097496e-05
## 17932 USP17L22 0 0 2 0 0 9729.2562 9729.2562 0.9446037 6.287852e-08 4.273784e-07 0.9855767 5.973446e-04 3.980201e-03
## 11565 PABPC3 0 0 2 0 0 8019.2737 8019.2737 0.9377189 9.325738e-08 6.213885e-07 0.9855767 5.973446e-04 3.980201e-03
## 18196 WDPCP 0 0 2 0 0 3805.0805 3805.0805 0.9349721 4.351865e-07 2.744960e-06 0.9855767 2.090636e-03 1.318679e-02
## 9776 MROH2A 0 0 1 1 0 1761.9124 1761.9124 0.8956018 2.096382e-06 1.176643e-05 0.9855767 8.056814e-03 4.522073e-02
## 8301 KMT2C 0 0 1 1 0 716.5142 716.5142 0.8386098 1.257247e-05 5.867674e-05 0.9855767 2.906275e-02 1.557393e-01
## ......
? gene_name: 基因的名称。
? n_syn: 同义突变的数量。
? n_mis: 错义突变的数量。
? n_non: 无义突变的数量。
? n_spl: 剪接位点突变的数量。
? n_ind: 插入/缺失突变的数量。
? wmis_cv: 错义突变的dN/dS比值（最大似然估计）。
? wnon_cv: 无义突变的dN/dS比值（最大似然估计）。
? wspl_cv: 剪接位点突变的dN/dS比值（最大似然估计）。
? pmis_cv: 错义突变的P值，表示该基因中错义突变的选择压力的显著性。
? ptrunc_cv: 截断突变（无义和剪接位点突变）的P值，表示该基因中截断突变的选择压力的显著性。
? pallsubs_cv: 所有非同义突变（错义、无义和剪接位点突变）的P值，表示该基因中所有非同义突变的选择压力的显著性。
? qmis_cv: 错义突变的q值（FDR校正后的P值）。
? qtrunc_cv: 截断突变（无义和剪接位点突变）的q值（FDR校正后的P值）。
? qallsubs_cv: 所有非同义突变（错义、无义和剪接位点突变）的q值（FDR校正后的P值）。
3）突变注释（$annotmuts）：这个表包含了所有输入的突变及其注释信息
dndsout[["annotmuts"]]
## sampleID chr pos ref mut geneind gene ref_cod mut_cod strand ref3_cod mut3_cod aachange ntchange impact pid
## 204 CC56tissue 1 924949 G A 14271 SAMD11 G A 1 GGT GAT - - Essential_Splice ENSP00000480678
## 10550 CC56tissue 1 12716215 C A 12 AADACL3 C A 1 GCG GAG C13* C39A Nonsense ENSP00000352268
## 12672 CC56tissue 1 15626753 C T 4093 DDI2 C T 1 ACG ATG R75* C223T Nonsense ENSP00000417748
## 22914 CC56tissue 1 26553429 T C 14086 RPS6KA1 T C 1 CTG CCG A169A T507C Synonymous ENSP00000363283
## 34888 CC56tissue 1 44941527 G A 4826 EIF2B3 C T -1 TCA TTA Q145* C433T Nonsense ENSP00000353575
## 35452 CC56tissue 1 46032703 A G 9269 MAST2 A G 1 TAG TGG V1031V A3093G Synonymous ENSP00000501318
## ......
? sampleID: 样本的ID。
? chr: 突变所在的染色体。
? pos: 突变的位置。
? ref: 参考碱基。
? mut: 突变碱基。
? gene: 突变所在的基因。
? ref_cod: 参考序列的三联密码子（codon）。
? mut_cod: 突变后的三联密码子（codon）。
? strand: 基因所在的链（正链或负链）。
? ref3_cod: 参考序列的三联密码子（codon）及其上下游核苷酸序列。
? mut3_cod: 突变后的三联密码子（codon）及其上下游核苷酸序列。
? aachange: 氨基酸变化（从参考氨基酸到突变后的氨基酸）。
? ntchange: 核苷酸变化（从参考核苷酸到突变后的核苷酸）。
? impact: 突变的影响（如synonymous、missense、nonsense等）。
4）局部中性测试（$sel_loc）：包含了局部选择分析的结果，用于识别基因中特定位置的突变是否显示出显著的选择信号。这在识别可能具有功能性影响的热点区域时非常有用
dndsout[["sel_loc"]]
## gene_name n_syn n_mis n_non n_spl wmis_loc wnon_loc wspl_loc pmis_loc pall_loc qmis_loc qall_loc
## 15213 SLC9B1 0 0 8 0 0 10000 0 0.9999991 8.466339e-12 0.9999998 1.626892e-07
## 11565 PABPC3 0 0 6 0 0 10000 0 0.9999989 5.799150e-08 0.9999998 5.571824e-04
## 2129 CALHM4 0 0 0 3 0 0 10000 0.9999993 1.035322e-06 0.9999998 6.631585e-03
## 6152 GCNT7 0 0 0 3 0 0 10000 0.9999991 7.769189e-06 0.9999998 3.732318e-02
## 18819 ZNF419 0 0 0 3 0 0 10000 0.9999991 3.321557e-05 0.9999998 1.276541e-01
## 17350 TRIM51GP 0 0 0 3 0 0 10000 0.9999991 7.640284e-05 0.9999998 1.503908e-01
## ......
signif_genes_localmodel = as.vector(dndsout$sel_loc$gene_name[dndsout$sel_loc$qall_loc<0.1])
print(signif_genes_localmodel)
## [1] "SLC9B1" "PABPC3" "CALHM4" "GCNT7"
? 该模型检测出SLC9B1, PABPC3, CALHM4, GCNT7有显著突变。
