qgg包：一款大规模数量遗传和基因组分析的R包
一、概述
该包基于：假设基因组特征可能会富集影响性状的因果变体。根据以往的研究和不同的信息来源，可以分成几种基因组特征,如基因、染色体或生物途径。
原文：http://psoerensen.github.io/qgg/
1、核心功能
拟合线性混合模型
构建基因组关系矩阵
估计遗传参数(遗传性和相关性)
基因预测
单标记关联分析
基因集合富集分析
2、qgg利用以下处理大规模数据
使用openMP的多核处理
在BLAS库(如OpenBLAS、ATLAS或MKL)中实现的多线程矩阵操作
对存储在二进制文件(如PLINK bedfiles)中的基因型数据进行快速、内存高效的批处理
qgg包提供了一系列基因组特征建模方法，包括基因组特征最佳线性无偏预测(GFBLUP)模型，使用似然或贝叶斯方法。这些模型可以包含多种特征和多种性状，并可以使用不同的遗传模型(如加性、显性、基因间的遗传和基因间的环境相互作用)。GFBLUP模型使用个体遗传标记关系的差异加权。可以执行计算非常快的标记集测试。这些标记集测试允许快速分析不同层次的基因组特征类，以发现基因组特征潜在丰富的因果变异。因此，标记集测试可以促进更准确的预测模型。
二、安装和使用
1、安装
install.packages("devtools")
library(devtools)
options(devtools.install.args=" --no-multiargs")
devtools::install_github("psoerensen/qgg")
出错时
library(remotes)
install_url(url="https://github.com/psoerensen/qgg/archive/master.zip", INSTALL_opt= "--no-multiarch")
2、使用示范
小数据集教程
包含获取和清理数据到基因组参数的估计、基因集富集分析和使用基因组最佳线性无偏预测(GBLUP)和基因组特征BLUP (GFBLUP)模型进行预测。
以表型和遗传数据来自公开的果蝇遗传参考图(DGRP)——一个由205个完全自交系的黑anogaster系组成的系统为例。
①数据准备
②数据分析
大数据集教程
这些关于两个大型数据集的简明教程举例说明了qgg的一些核心特性。从获取和清理数据到基因组参数的估计、基因集富集分析和使用基因组最佳线性无偏预测(GBLUP)和基因组特征BLUP (GFBLUP)模型进行预测。
3、所有函数
| function | explain |
|---|---|
| bgfm() | 使用贝叶斯方法(小数据)实现基因组特征模型分析 |
| gprep() | 为所有统计分析准备基因型数据(第一步) |
| greml() | 基因组REML分析 |
| grm() | 计算基因组关系矩阵(GRM) |
| gscore() | 基于单标记汇总统计的基因组预测 |
| gsea() | 基因集合富集分析 |
| gsolve() | 基于线性混合模型的基因组预测 |
| lma() | 单标记关联分析使用线性模型或线性混合模型 |
| qgg | 使用似然或贝叶斯方法实现基因组特征线性混合模型 |
续：一、小数据集教程
1、准备表型和协变量数据
在这个脚本中，我们准备了表型和协变量数据框架，以用于后续基因组分析。作为表型数据的一个例子，我们使用从果蝇遗传参考图(DGRP)的表型“耐饥饿”。数据可以在http://dgrp2.gnets.ncsu.edu/data.html的“表型文件”，“Mackay, et al.， Nature, 2012”下获得。使用了雄性和雌性的数据。倒置状态(染色体倒置)和沃尔巴克氏感染状态(沃尔巴克氏感染状态)可以在页面底部的“其他有用文件”标题下找到。
①使用plyr、dplyr和tidyr包编辑数据。readxl包用于读取.xlsx文件
安装这些包:
#install.packages("plyr")
#install.packages("dplyr")
#install.packages("tidyr")
#install.packages("readxl")
library(plyr)
library(dplyr)
library(tidyr)
library(readxl)
②下载表型和协变量数据
# Female
download.file("http://dgrp2.gnets.ncsu.edu/data/website/starvation.female.csv",
destfile = "./data/starvation.female.csv")
# Male
download.file("http://dgrp2.gnets.ncsu.edu/data/website/starvation.male.csv",
destfile = "./data/starvation.male.csv")
# Inverstion status
download.file("http://dgrp2.gnets.ncsu.edu/data/website/inversion.xlsx", mode = "wb",
destfile = "./data/inversion.xlsx")
# Wolbachia
download.file("http://dgrp2.gnets.ncsu.edu/data/website/wolbachia.xlsx", mode = "wb",
destfile = "./data/wolbachia.xlsx")
③读取和编辑表型数据
读取雌性数据
starF <- read.csv(file="./data/starvation.female.csv", header = FALSE)
head(starF)
## V1 V2
## 1 line_100 77.92000
## 2 line_101 57.76000
## 3 line_105 73.12000
## 4 line_109 53.44000
## 5 line_129 42.77551
## 6 line_136 104.32000
dim(starF) # 获取矩阵维度
## [1] 203 2
读取雄性数据
starM <- read.csv(file="./data/starvation.male.csv", header = FALSE)
head(starM)
## V1 V2
## 1 line_100 49.28000
## 2 line_101 47.20000
## 3 line_105 51.04000
## 4 line_109 44.96000
## 5 line_129 33.08475
## 6 line_136 63.04000
dim(starM)
## [1] 203 2
给列名称。“L”=lines，“F”=雌性，“M”=雄性
colnames(starF) <- c("L", "F")
colnames(starM) <- c("L", "M")
在dplyr中，必须将数据帧转换为tibble (tbl)。把starF和starM转换成tibbles。
starF <- tbl_df(starF)
starM <- tbl_df(starM)
查看tibble
starF
## # A tibble: 203 x 2
## L F
## <fct> <dbl>
## 1 line_100 77.9
## 2 line_101 57.8
## 3 line_105 73.1
## 4 line_109 53.4
## 5 line_129 42.8
## 6 line_136 104.
## 7 line_138 59.5
## 8 line_142 59.3
## 9 line_149 47
## 10 line_153 59.0
## # ... with 193 more rows
starM
## # A tibble: 203 x 2
## L M
## <fct> <dbl>
## 1 line_100 49.3
## 2 line_101 47.2
## 3 line_105 51.0
## 4 line_109 45.0
## 5 line_129 33.1
## 6 line_136 63.0
## 7 line_138 47.8
## 8 line_142 38.4
## 9 line_149 35.8
## 10 line_153 40.3
## # ... with 193 more rows
Join the tibbles for males and females.
starMF <- left_join(starM, starF, by= "L")
starMF
## # A tibble: 203 x 3
## L M F
## <fct> <dbl> <dbl>
## 1 line_100 49.3 77.9
## 2 line_101 47.2 57.8
## 3 line_105 51.0 73.1
## 4 line_109 45.0 53.4
## 5 line_129 33.1 42.8
## 6 line_136 63.0 104.
## 7 line_138 47.8 59.5
## 8 line_142 38.4 59.3
## 9 line_149 35.8 47
## 10 line_153 40.3 59.0
## # ... with 193 more rows
创建一列性别信息和一列表型(y)，在这里是耐饥饿。
starv <- gather(starMF, sex, y, -L)
head(starv)
## # A tibble: 6 x 3
## L sex y
## <fct> <chr> <dbl>
## 1 line_100 M 49.3
## 2 line_101 M 47.2
## 3 line_105 M 51.0
## 4 line_109 M 45.0
## 5 line_129 M 33.1
## 6 line_136 M 63.0
head(starv, 3)
## # A tibble: 3 x 3
## L sex y
## <fct> <chr> <dbl>
## 1 line_100 M 49.3
## 2 line_101 M 47.2
## 3 line_105 M 51.0
从“L”列的内容中删除前缀“line_”。
starv$L <- gsub("line_", "", starv$L, fixed = TRUE)
head(starv$L)
## [1] "100" "101" "105" "109" "129" "136"
读取和编辑倒置状态
所使用的缩写是:INV =倒置核型(纯合)，INV / ST =倒置的杂合子，ST =纯合形式的标准配置。
inv <- read_excel("./data/inversion.xlsx",
sheet = 1, col_names = TRUE)
head(inv)
## # A tibble: 6 x 17
## `DGRP Line` `In(2L)t` `In(2R)NS` `In(2R)Y1` `In(2R)Y2` `In(2R)Y3`
## <chr> <chr> <chr> <chr> <chr> <chr>
## 1 DGRP_21 ST ST ST ST ST
## 2 DGRP_26 INV ST ST ST ST
## 3 DGRP_28 ST INV ST ST ST
## 4 DGRP_31 ST ST ST ST ST
## 5 DGRP_32 INV ST ST ST ST
## 6 DGRP_38 ST ST ST ST ST
## # ... with 11 more variables: `In(2R)Y4` <chr>, `In(2R)Y5` <chr>,
## # `In(2R)Y6` <chr>, `In(2R)Y7` <chr>, `In(3L)P` <chr>, `In(3L)M` <chr>,
## # `In(3L)Y` <chr>, `In(3R)P` <chr>, `In(3R)K` <chr>, `In(3R)Mo` <chr>,
## # `In(3R)C` <chr>
dim(inv)
## [1] 205 17
从“DGRP行”列的内容中删除“DGRP_”前缀。将Inv的列名保存为vector Invcols。然后编辑Invcols向量中的列名:“DGRP行”变为“L”，倒置名称中的括号“()”变为下划线。
inv$`DGRP Line` <- gsub("DGRP_", "", inv$`DGRP Line`, fixed=TRUE)
head(inv, 3)
## # A tibble: 3 x 17
## `DGRP Line` `In(2L)t` `In(2R)NS` `In(2R)Y1` `In(2R)Y2` `In(2R)Y3`
## <chr> <chr> <chr> <chr> <chr> <chr>
## 1 21 ST ST ST ST ST
## 2 26 INV ST ST ST ST
## 3 28 ST INV ST ST ST
## # ... with 11 more variables: `In(2R)Y4` <chr>, `In(2R)Y5` <chr>,
## # `In(2R)Y6` <chr>, `In(2R)Y7` <chr>, `In(3L)P` <chr>, `In(3L)M` <chr>,
## # `In(3L)Y` <chr>, `In(3R)P` <chr>, `In(3R)K` <chr>, `In(3R)Mo` <chr>,
## # `In(3R)C` <chr>
invcols <- colnames(inv)
invcols[1] <- "L"
invcols[2:17] <- gsub("(", "", invcols[2:17], fixed = TRUE)
invcols[2:17] <- gsub(")", "", invcols[2:17], fixed = TRUE)
colnames(inv) <- invcols
inv
## # A tibble: 205 x 17
## L In2Lt In2RNS In2RY1 In2RY2 In2RY3 In2RY4 In2RY5 In2RY6 In2RY7
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 21 ST ST ST ST ST ST ST ST ST
## 2 26 INV ST ST ST ST ST ST ST ST
## 3 28 ST INV ST ST ST ST ST ST ST
## 4 31 ST ST ST ST ST ST ST ST ST
## 5 32 INV ST ST ST ST ST ST ST ST
## 6 38 ST ST ST ST ST ST ST ST ST
## 7 40 ST ST ST ST ST ST ST ST ST
## 8 41 ST ST ST ST ST ST ST ST ST
## 9 42 ST ST ST ST ST ST ST ST ST
## 10 45 ST ST ST ST ST ST ST ST ST
## # ... with 195 more rows, and 7 more variables: In3LP <chr>, In3LM <chr>,
## # In3LY <chr>, In3RP <chr>, In3RK <chr>, In3RMo <chr>, In3RC <chr>
阅读和编辑沃尔巴克状态
wo <- read_excel("./data/wolbachia.xlsx",
sheet = 1, col_names = TRUE)
wo
## # A tibble: 205 x 2
## `DGRP Line` `Infection Status`
## <chr> <chr>
## 1 DGRP__21 y
## 2 DGRP__26 n
## 3 DGRP__28 n
## 4 DGRP__31 n
## 5 DGRP__32 n
## 6 DGRP__38 n
## 7 DGRP__40 y
## 8 DGRP__41 n
## 9 DGRP__42 n
## 10 DGRP__45 n
## # ... with 195 more rows
dim(wo)
## [1] 205 2
更改wo的列名
colnames(wo) <- c("L", "wo")
wo$L <- gsub("DGRP__", "", wo$L, fixed=TRUE)
wo
## # A tibble: 205 x 2
## L wo
## <chr> <chr>
## 1 21 y
## 2 26 n
## 3 28 n
## 4 31 n
## 5 32 n
## 6 38 n
## 7 40 y
## 8 41 n
## 9 42 n
## 10 45 n
## # ... with 195 more rows
创建最终数据框
将表型数据与反转状态和沃尔巴克氏体感染状态合并。
starvInv <- left_join(starv, inv, by= "L")
starvInv
## # A tibble: 406 x 19
## L sex y In2Lt In2RNS In2RY1 In2RY2 In2RY3 In2RY4 In2RY5
## <chr> <chr> <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 100 M 49.3 INV/ST ST ST ST ST ST ST
## 2 101 M 47.2 INV/ST ST ST ST ST ST ST
## 3 105 M 51.0 ST ST ST ST ST ST ST
## 4 109 M 45.0 INV/ST ST ST ST ST ST ST
## 5 129 M 33.1 ST ST ST ST ST ST ST
## 6 136 M 63.0 ST ST ST ST ST ST ST
## 7 138 M 47.8 ST ST ST ST ST ST ST
## 8 142 M 38.4 ST ST ST ST ST ST ST
## 9 149 M 35.8 ST ST ST ST ST ST ST
## 10 153 M 40.3 ST ST ST ST ST ST ST
## # ... with 396 more rows, and 9 more variables: In2RY6 <chr>,
## # In2RY7 <chr>, In3LP <chr>, In3LM <chr>, In3LY <chr>, In3RP <chr>,
## # In3RK <chr>, In3RMo <chr>, In3RC <chr>
starvIW <- left_join(starvInv, wo, by="L")
starvIW[1:5,15:20]
## # A tibble: 5 x 6
## In3LY In3RP In3RK In3RMo In3RC wo
## <chr> <chr> <chr> <chr> <chr> <chr>
## 1 ST ST INV ST ST y
## 2 ST ST ST ST ST n
## 3 ST ST INV ST ST n
## 4 ST ST ST ST ST n
## 5 ST ST ST ST ST n
这里一个如何查看数据的示例。我们展示了前三个逆序的摘要(以表格形式)。
apply(starvIW[,4:19], 2, table)[1:3]
## $In2Lt
##
## INV INV/ST ST
## 38 50 318
##
## $In2RNS
##
## INV INV/ST ST
## 14 20 372
##
## $In2RY1
##
## INV/ST ST
## 2 404
将tibble转换为数据框并保存已编辑的表现型数据。
starv <- as.data.frame(starvIW)
save(starv, file="./phenotypes/starv_inv_wo.Rdata")
2.准备基因型数据
主要包括两个部分:1)下载并编辑原始基因型数据;2.)准备一个集中和缩放的基因型矩阵($用于qgg的后续统计分析。
①下载和编辑基因型数据
在这个脚本中，我们准备了用于后续基因组分析的基因型数据的数据框。以果蝇遗传参考面图(DGRP)中的基因型为例。基因型数据可以在http://dgrp2.gnets.ncsu.edu/data.html页面顶部“基因型文件”的标题下获得。使用“表格格式基因型”数据。
download.file("http://dgrp2.gnets.ncsu.edu/data/website/dgrp2.tgeno",
destfile = "./data/dgrp2.tgeno")
genotypes <- read.table( "./data/dgrp2.tgeno",
header = TRUE, sep="")
dim(genotypes)
## [1] 4438427 214
genotypes[1:5,1:12]
## chr pos id ref alt refc altc qual cov line_21 line_26 line_28
## 1 2L 4998 2L_4998_SNP G A 117 5 999 12 0 - 0
## 2 2L 5002 2L_5002_SNP G T 127 1 999 13 0 - 0
## 3 2L 5039 2L_5039_SNP C T 1 118 999 21 2 - 2
## 4 2L 5040 2L_5040_SNP G A 1 118 999 21 2 - 2
## 5 2L 5092 2L_5092_SNP C T 6 119 999 22 2 - 2
数据以空格分隔，等位基因编码为‘0’作为参考等位基因，‘2’作为替代等位基因(不包括小等位基因)。杂合或不能可靠地称为基因型(低质量的碱基呼叫、低质量的定位或浅读深度)的基因型用“-”标注。
编辑数据框
编辑列名称:删除每一行的“line_”，使用id列(SNP id)作为行名，并将类更改为字符。
colnames(genotypes) <- gsub("line_", "",colnames(genotypes))
rownames(genotypes) <- as.character(genotypes$id)
②子集基因型数据-只包括SNP变异
创建一个载体包含变异类型的遗传标记，即。、单核苷酸多聚性(SNP)、插入(INS)、删除(DEL)和微卫星(MNP)。id列的第三个组件(由“_”分隔)表示变体类型。
vtype <- sapply(as.character(genotypes$id), function(x){
strsplit(x,split="_")[[1]][3]
} )
unique(vtype)
## [1] "SNP" "INS" "DEL" "MNP"
准备一个只包含SNP变异的基因型数据子集。插入、删除和微卫星变体的子集也可以用类似的方法制备。
gsnp <- genotypes[vtype=="SNP",]
dim(gsnp)
## [1] 3963420 214
③为SNP信息和基因型创建数据框
将基因型数据框拆分成2个数据框:
1). 一个数据框(snpI)，包含关于每个变体的信息(例如，变体ID、染色体、在染色体上的位置、参考(ref)/替代(alt)等位基因、ref/alt等位基因计数、phred质量分数等)。
2). 包含205个DGRP株系每个变种的基因型的数据框(snpG)。
snpI <- gsnp[,1:9]
snpG <- gsnp[,-c(1:9)]
从snpI和snpG数据帧的行名中删除“SNP”后缀。
rownames(snpI) <- gsub("_SNP","",rownames(snpI))
rownames(snpG) <- gsub("_SNP","",rownames(snpG))
head(rownames((snpG)))
## [1] "2L_4998" "2L_5002" "2L_5039" "2L_5040" "2L_5092" "2L_5095"
dim(snpI)
## [1] 3963420 9
head(snpI)
## chr pos id ref alt refc altc qual cov
## 2L_4998 2L 4998 2L_4998_SNP G A 117 5 999 12
## 2L_5002 2L 5002 2L_5002_SNP G T 127 1 999 13
## 2L_5039 2L 5039 2L_5039_SNP C T 1 118 999 21
## 2L_5040 2L 5040 2L_5040_SNP G A 1 118 999 21
## 2L_5092 2L 5092 2L_5092_SNP C T 6 119 999 22
## 2L_5095 2L 5095 2L_5095_SNP T A 4 115 999 22
dim(snpG)
## [1] 3963420 205
str(snpG[1:5])
## 'data.frame': 3963420 obs. of 5 variables:
## $ 21: Factor w/ 3 levels "-","0","2": 2 2 3 3 3 3 2 2 2 2 ...
## $ 26: Factor w/ 3 levels "-","0","2": 1 1 1 1 1 1 2 2 2 2 ...
## $ 28: Factor w/ 3 levels "-","0","2": 2 2 3 3 3 3 2 2 2 2 ...
## $ 31: Factor w/ 3 levels "-","0","2": 2 2 3 3 3 3 2 2 2 3 ...
## $ 32: Factor w/ 3 levels "-","0","2": 1 1 1 1 1 1 2 2 2 2 ...
snpG[1:5,1:15]
## 21 26 28 31 32 38 40 41 42 45 48 49 57 59 69
## 2L_4998 0 - 0 0 - 0 0 0 - 0 0 0 - - 0
## 2L_5002 0 - 0 0 - 0 0 0 - 0 0 0 - - 0
## 2L_5039 2 - 2 2 - 2 2 2 - 2 - 2 - - 2
## 2L_5040 2 - 2 2 - 2 2 2 - 2 - 2 - - 2
## 2L_5092 2 - 2 2 - 2 2 2 - 2 2 2 - - 2
到目前为止，运行脚本需要一些时间，而且它还填充了R环境。在这个阶段，保存数据框snpI和snpG可能比较方便。这一步是可选的。
save(snpG,file="./genotypes/snpG.Rdata")
save(snpI,file="./genotypes/snpI.Rdata")
这里我们清理R环境，因为R在执行一个已填充的环境后可能会遇到问题。
rm(list = ls(all=TRUE))
gc()
## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 527583 28.2 11196372 598.0 13995465 747.5
## Vcells 9353552 71.4 1333548262 10174.2 2070267005 15794.9
④对SNP信息和基因型数据框进行过滤
编辑SNP基因型数据(snpG)和SNP信息(snpI):根据Phred缩放的变异质量、基因型召唤率和小等位基因频率标准对数据进行过滤。
加载snpG和snpI数据框。如果没有按照上述建议清除R环境，则忽略此步骤。
load(file = "./genotypes/snpG.Rdata")
load(file = "./genotypes/snpI.Rdata")
⑤phr缩放质量
列出Phred缩放质量(qual变量)大于或等于500的snp。
phred500 <- snpI$qual>=500
⑥基因型检出率
选择基因型检出率为80%的snp。得分为snpI$altc的备选等位基因数和得分为snpI$refc的参考等位基因数(保存为nC)之和，就得到了可以在所有系中为一个SNP得分的基因型数。在可能的205 = 80%的检出中，164个基因型得分。
nC80是一个逻辑向量，其中snpI$altc和snpI$refc之和大于或等于164的snp为真。
nC <- snpI$altc + snpI$refc
nC80 <- nC>=164
⑦次要等位基因频率
去除次要等位基因频率(MAF) < 0.05的snp。mafA 为替代等位基因的频率，mafR 为参考等位基因的频率。生成DGRP 20代全同胞交配产生的近似近亲系数0.986,导致大多数隔离网站是纯合子(编码为0或2)。杂合的基因型,或不能被称为可靠(低质量的基础,低质量的映射,或者浅阅读深度)由“-”注释。
mafR <- snpI$refc/nC
mafA <- snpI$altc/nC
head(mafR)
## [1] 0.959016393 0.992187500 0.008403361 0.008403361 0.048000000 0.033613445
head(mafA)
## [1] 0.04098361 0.00781250 0.99159664 0.99159664 0.95200000 0.96638655
创建逻辑，如果maf是>= 0.05为真，否则为假:替代等位基因为mafA05，参考等位基因为mafR05。
mafA05 <- mafA>=0.05
mafR05 <- mafR>=0.05
maf05 <- 0.05<=mafA和mafA<=0.95也可以选择频率大于0.05的等位基因。
⑧子集相关的单核苷酸多态性
创建逻辑向量，这些逻辑向量对于snp来说是真实的，满足所有关于phred缩放质量，呼叫率和MAF的标准。基于此向量对包含所选snp子集的snpI和snpG数据框进行子集划分。
keep <- phred500 & nC80 & mafA05 & mafR05
sum(keep)
## [1] 1725755
head(rownames(snpI)[keep])
## [1] "2L_5317" "2L_5372" "2L_5390" "2L_5403" "2L_5465" "2L_5598"
snpG <- snpG[keep,]
snpI <- snpI[keep,]
将数据框snpG转换为矩阵。矩阵的内容是字符。
snpG <- as.matrix(snpG)
head(snpG[,1:20])
## 21 26 28 31 32 38 40 41 42 45 48 49 57 59 69 73
## 2L_5317 "0" "0" "0" "2" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
## 2L_5372 "0" "0" "0" "0" "0" "2" "2" "0" "2" "-" "2" "0" "0" "-" "-" "2"
## 2L_5390 "2" "0" "2" "0" "0" "2" "2" "2" "2" "-" "2" "-" "-" "-" "0" "2"
## 2L_5403 "0" "0" "0" "0" "0" "0" "0" "0" "0" "2" "2" "0" "0" "0" "0" "2"
## 2L_5465 "0" "0" "0" "0" "0" "0" "0" "0" "0" "2" "2" "0" "0" "0" "0" "2"
## 2L_5598 "0" "0" "0" "0" "0" "0" "0" "0" "0" "2" "0" "0" "0" "0" "0" "2"
## 75 83 85 88
## 2L_5317 "0" "0" "-" "0"
## 2L_5372 "2" "2" "-" "2"
## 2L_5390 "2" "2" "-" "2"
## 2L_5403 "2" "2" "0" "0"
## 2L_5465 "2" "2" "0" "0"
## 2L_5598 "-" "2" "0" "0"
将snpG的矩阵内容转换为数值。apply中的第二个参数(数字2)指定按列而不是按行工作。
snpG <- apply(snpG,2,as.numeric)
rownames(snpG) <- rownames(snpI)
head(snpG[,1:20])
## 21 26 28 31 32 38 40 41 42 45 48 49 57 59 69 73 75 83 85 88
## 2L_5317 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 NA 0
## 2L_5372 0 0 0 0 0 2 2 0 2 NA 2 0 0 NA NA 2 2 2 NA 2
## 2L_5390 2 0 2 0 0 2 2 2 2 NA 2 NA NA NA 0 2 2 2 NA 2
## 2L_5403 0 0 0 0 0 0 0 0 0 2 2 0 0 0 0 2 2 2 0 0
## 2L_5465 0 0 0 0 0 0 0 0 0 2 2 0 0 0 0 2 2 2 0 0
## 2L_5598 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 2 NA 2 0 0
将编辑过的SNP数据保存为.rdata和.rds文件。在compute_W脚本中加载.rds文件。更多信息在?saveRDS.
save(snpG,file="./genotypes/snpGE.Rdata")
save(snpI,file="./genotypes/snpIE.Rdata")
saveRDS(snpG, file = "./genotypes/snpGE.rds")
saveRDS(snpI,file="./genotypes/snpIE.rds")
准备集中和缩放的基因型矩阵W
这里提供了两种准备W的方法。
方法1
加性基因组关系矩阵G 使用所有遗传标记构建如下:
G=WW′/m ,其中W为居中缩放的基因型矩阵,m是标记的总数,W的各列向量计算如下:
式中pi为第i个遗传标记的小等位基因频率，mi为等位基因计数矩阵的第i列向量，M，包含编码为0、1或2的基因型，计数小等位基因的数量。(对于几乎纯合的DGRP系，基因型编码为0或2。)
加载编辑好的SNP基因型数据文件snpGE。Rdata，之前在qgg用户指南中创建的。基因型数据帧snpG由readRDS()函数直接加载到对象W中。
W <- readRDS(file="./genotypes/snpGE.rds")
计算每个SNP的小等位基因数(nMinor)和总等位基因数(n等位基因)。实际上，对于n等位基因，我们计算每行/单核苷酸多态性测定的基因型的数量(因此不是“NA’s”)。这相当于每个基因型计数一个等位基因。因此等位基因的总数是n等位基因的2倍。
nMinor <- rowSums(W, na.rm=TRUE)
nAlleles <- rowSums(!is.na(W))
计算小等位基因频率:
p<-nMinor/(2*nAlleles)
min(p)
max(p)
使用观察到的等位基因频率得到W:
for ( i in 1:205) {
W[,i] <- (W[,i]-2*p)/sqrt(2*p*(1-p))
isNA <- is.na(W[,i])
W[isNA,i] <- 0
}
W <- t(W
保存居中和缩放W计算方法1为dgrp2_W1.Rdata。
save(W, file="./genotypes/dgrp2_W1.Rdata")
方法2
在这种方法中，使用scale()函数对W中的列进行缩放，该函数的默认方法为中心和/或对数字矩阵的列进行缩放。
rm(list=ls(all=TRUE))
加载编辑好的SNP基因型数据文件:
W <- readRDS(file="./genotypes/snpGE.rds")
W被转置，因为scale()函数默认情况下对数值矩阵的列进行缩放。
W <- t(W)
W <- scale(W)
dim(W)
设置缺失值为0。
for ( i in 1:205) {
isNA <- is.na(W[i,])
W[i, isNA] <- 0
}
保存居中和缩放W计算方法2为dgrp2_w . rdata。
save(W, file="./genotypes/dgrp2_W2.Rdata")
3.准备注释数据和标记集
以下说明如何根据基因组特征将遗传标记(例如snp)分组。遗传标记集可用于基因集富集分析或作为GFBLUP预测模型中的“特征”。这个例子展示了如何通过序列本体、flybase (fb)基因id、染色体和基因本体术语来定义snp集。
遗传标记注释数据(即遗传标记与基因、染色体和序列本体之间的映射)来自公开的果蝇遗传参考面板(DGRP)。基因本体论注释(基因与基因本体论术语之间的映射)来自生物导体包，org.Dm.eg.db，内含黑腹果蝇全基因组注释。全基因组注释主要基于使用Entrez基因标识符的绘图。
①下载并读取注释数据
变体注释(基于FB5.49)可以在http://dgrp2.gnets.ncsu.edu/data.html页面底部的“其他有用文件”下找到。
download.file("http://dgrp2.gnets.ncsu.edu/data/website/dgrp.fb549.annot.txt",
destfile = "./data/dgrp.fb549.annot.txt")
annotation <- read.table(file = "./data/dgrp.fb549.annot.txt", sep="\t",
colClasses="character", quote="")
该文件以制表符分隔。在site class栏(第3栏)中使用的分隔器在这里解释:http://dgrp2.gnets.ncsu.edu/faq.html 在项目 “3.输出文件是什么?GeneAnnotation”。
dim(annotation)
## [1] 4438427 4
str(annotation)
## 'data.frame': 4438427 obs. of 4 variables:
## $ V1: chr "2L_10000016_SNP" "2L_10000023_SNP" "2L_10000029_SNP" "2L_10000033_SNP" ...
## $ V2: chr "C" "C" "G" "G" ...
## $ V3: chr "SiteClass[FBgn0051875|CG31875|INTRON|0;FBgn0051755|SoYb|NON_SYNONYMOUS_CODING|0],TranscriptAnnot[INTRON(MODIFIE"| __truncated__ "SiteClass[FBgn0051875|CG31875|INTRON|0;FBgn0051755|SoYb|NON_SYNONYMOUS_CODING|0],TranscriptAnnot[INTRON(MODIFIE"| __truncated__ "SiteClass[FBgn0051875|CG31875|INTRON|0;FBgn0051755|SoYb|NON_SYNONYMOUS_CODING|0],TranscriptAnnot[INTRON(MODIFIE"| __truncated__ "SiteClass[FBgn0051875|CG31875|INTRON|0;FBgn0051755|SoYb|SYNONYMOUS_CODING|0],TranscriptAnnot[INTRON(MODIFIER|||"| __truncated__ ...
## $ V4: chr "-" "-" "-" "-" ...
②编辑注释数据框
注释数据框架的第一列用作行名称。as.character()用于避免水平因素。
rownames(annotation) <- as.character(annotation[,1])
使用sapply将每个行名(例如2L_10000016_SNP)拆分为下划线“_”。得到2L 10000016 SNP。行名的第三个元素([3]，这里是“SNP”)被保存在一个名为vtype(变体类型)的向量中。
vtype <- sapply(rownames(annotation), function(x){
strsplit(x,split="_")[[1]][3]
} )
head(vtype)
## 2L_10000016_SNP 2L_10000023_SNP 2L_10000029_SNP 2L_10000033_SNP
## "SNP" "SNP" "SNP" "SNP"
## 2L_10000089_SNP 2L_10000133_SNP
## "SNP" "SNP"
注释数据框架(SNP id)的Rownames作为vtype向量中每个元素的名称保存。
head(names(vtype))
## [1] "2L_10000016_SNP" "2L_10000023_SNP" "2L_10000029_SNP" "2L_10000033_SNP"
## [5] "2L_10000089_SNP" "2L_10000133_SNP"
Look at unique terms in vtype.
unique(vtype)
## [1] "SNP" "DEL" "INS" "MNP"
snpA_raw是一个包含数据框注释的第3列(site类)信息的向量，只包含snp。看一下snpA_raw向量的第一个元素。
snpA_raw <- as.character(annotation[vtype=="SNP",3])
length(snpA_raw)
## [1] 3963420
snpA_raw[1]
## [1] "SiteClass[FBgn0051875|CG31875|INTRON|0;FBgn0051755|SoYb|NON_SYNONYMOUS_CODING|0],TranscriptAnnot[INTRON(MODIFIER||||161|CG31875|protein_coding|CODING|FBtr0303883|3);INTRON(MODIFIER||||161|CG31875|protein_coding|CODING|FBtr0303884|3);INTRON(MODIFIER||||161|CG31875|protein_coding|CODING|FBtr0308121|3);NON_SYNONYMOUS_CODING(MODERATE|MISSENSE|gaG/gaT|E471D|1476|SoYb|protein_coding|CODING|FBtr0303882|);NON_SYNONYMOUS_CODING(MODERATE|MISSENSE|gaG/gaT|E471D|1476|SoYb|protein_coding|CODING|FBtr0308122|)]"
给snpA_raw矢量命名:注释数据框的行名，即SNP id。从矢量名称中删除冗余的后缀“_SNP”，这样它就只包含SNP id。
names(snpA_raw) <- gsub("_SNP","",rownames(annotation)[vtype=="SNP"])
只使用“SiteClass”信息(即直接在“SiteClass”后面的方括号中的信息)。这包括flybase基因id、基因符号、映射序列本体(位点类)和到基因的碱基对距离，与转录本注释信息以“，”分隔.
(see: http://dgrp2.gnets.ncsu.edu/faq.html).
通过将snpA_raw拆分为“，”来选择站点类信息，并保留第一个字符串。从这个字符串中删除“SiteClass”和方括号，并在分号处拆分字符串。
snpA_raw <- lapply(snpA_raw, function(x) {
x <- strsplit(x,",")[[1]][1]
x <- gsub("SiteClass","",x)
x <- gsub("[","",x, fixed=TRUE)
x <- gsub("]","",x, fixed=TRUE)
x <- strsplit(x,";")[[1]]
x
})
head(snpA_raw)
## $`2L_10000016`
## [1] "FBgn0051875|CG31875|INTRON|0"
## [2] "FBgn0051755|SoYb|NON_SYNONYMOUS_CODING|0"
##
## $`2L_10000023`
## [1] "FBgn0051875|CG31875|INTRON|0"
## [2] "FBgn0051755|SoYb|NON_SYNONYMOUS_CODING|0"
##
## $`2L_10000029`
## [1] "FBgn0051875|CG31875|INTRON|0"
## [2] "FBgn0051755|SoYb|NON_SYNONYMOUS_CODING|0"
##
## $`2L_10000033`
## [1] "FBgn0051875|CG31875|INTRON|0"
## [2] "FBgn0051755|SoYb|SYNONYMOUS_CODING|0"
##
## $`2L_10000089`
## [1] "FBgn0051875|CG31875|INTRON|0"
## [2] "FBgn0051755|SoYb|NON_SYNONYMOUS_CODING|0"
##
## $`2L_10000133`
## [1] "FBgn0051875|CG31875|INTRON|0"
## [2] "FBgn0051755|SoYb|SYNONYMOUS_CODING|0"
length(snpA_raw)
## [1] 3963420
创建一个向量(nA_raw)，其中包含列表snpA_raw中每一行的长度。
nA_raw <- sapply(snpA_raw,length)
table(nA_raw)显示在flybase中包含1、2、3或最多12个单独注释的SNP位点的数量。
table(nA_raw)
## nA_raw
## 1 2 3 4 5 6 7 8 9
## 3241772 602432 100328 13818 2828 963 494 440 170
## 10 11 12
## 130 27 18
看看nA_raw中长度为12的18个元素中的前3个。
nA_raw[nA_raw==12][1:3]
## 2L_20418974 2L_20418976 2L_20418995
## 12 12 12
矢量snpNames包含snpA_raw的名称，这些名称重复出现在nA_raw中多次，即每个flybase注释一次。
snpNames <- rep(names(snpA_raw), times=nA_raw)
length(snpNames)
## [1] 4833131
取消snpA_raw列表并删除向量的名称，使向量只包含分隔的段。取消snpA_raw列表将创建一个与nA_raw长度相同的向量。
snpA_raw <- unlist(snpA_raw, use.names=FALSE)
head(snpA_raw)
## [1] "FBgn0051875|CG31875|INTRON|0"
## [2] "FBgn0051755|SoYb|NON_SYNONYMOUS_CODING|0"
## [3] "FBgn0051875|CG31875|INTRON|0"
## [4] "FBgn0051755|SoYb|NON_SYNONYMOUS_CODING|0"
## [5] "FBgn0051875|CG31875|INTRON|0"
## [6] "FBgn0051755|SoYb|NON_SYNONYMOUS_CODING|0"
length(snpA_raw)
## [1] 4833131
在“|”符号处分割snpA_raw的段。结果是一个列表，其中元素名称对应于上面准备的分隔段。“GBgn” (flybase gene id), gene id, mapped sequence ontology terms (“SeqOnt”) and base pair distance to gene (“distance”) 是每个元素的内容。如前所述，向量nA_raw给出了snpA_raw中每个元素的长度(每个元素分隔的段的数量)。table(nA_raw)显示，snpA_raw中的919806个元素包含3个片段，3913325个元素包含4个片段。
snpA_raw <- sapply(snpA_raw,function(x) { unlist(strsplit(x,split="|",fixed=TRUE)) })
nA_raw <- sapply(snpA_raw,length)
table(nA_raw)
## nA_raw
## 3 4
## 919806 3913325
创建一个空矩阵snpA，行数等于snpNames的长度,为4列。snpNames向量用作行名称。这对应于前面准备的snpA_raw列表的元素名。
snpA <- matrix(NA,nrow=length(snpNames),ncol=4)
rownames(snpA) <- snpNames
head(snpA)
## [,1] [,2] [,3] [,4]
## 2L_10000016 NA NA NA NA
## 2L_10000016 NA NA NA NA
## 2L_10000023 NA NA NA NA
## 2L_10000023 NA NA NA NA
## 2L_10000029 NA NA NA NA
## 2L_10000029 NA NA NA NA
Unlist snpA_raw,这创建了向量“FBgn”, gene names, mapped sequence ontology terms (“SeqOnt”) and base pair distance to gene (“distance”).通过对向量的转置，四种不同类型的信息分别放入到snpA的一列中。只包含包含4段的snpA_raw元素(因此长度为4的nA_raw元素)。
snpA[nA_raw==4,] <- t(sapply(snpA_raw[nA_raw==4],function(x) { unlist(x) }))
head(snpA)
## [,1] [,2] [,3] [,4]
## 2L_10000016 "FBgn0051875" "CG31875" "INTRON" "0"
## 2L_10000016 "FBgn0051755" "SoYb" "NON_SYNONYMOUS_CODING" "0"
## 2L_10000023 "FBgn0051875" "CG31875" "INTRON" "0"
## 2L_10000023 "FBgn0051755" "SoYb" "NON_SYNONYMOUS_CODING" "0"
## 2L_10000029 "FBgn0051875" "CG31875" "INTRON" "0"
## 2L_10000029 "FBgn0051755" "SoYb" "NON_SYNONYMOUS_CODING" "0"
最终的注释数据框将只包括居中和缩放的基因型矩阵(W)中的snp。
load(file="./genotypes/dgrp2_W2.Rdata")
准备一个逻辑向量，inW.TRUE取决于这个数据集中的snp是否也在W中发现。
inW <- rownames(snpA)%in%colnames(W)
只保留逻辑向量inW中正确的snp。
snpA <- snpA[inW,]
给出数据框snpA相关列名
colnames(snpA) <- c("FBid","GeneName","SeqOnt","distance")
snpA[1:5,]
## FBid GeneName SeqOnt distance
## 2L_10000016 "FBgn0051875" "CG31875" "INTRON" "0"
## 2L_10000016 "FBgn0051755" "SoYb" "NON_SYNONYMOUS_CODING" "0"
## 2L_10000033 "FBgn0051875" "CG31875" "INTRON" "0"
## 2L_10000033 "FBgn0051755" "SoYb" "SYNONYMOUS_CODING" "0"
## 2L_10000089 "FBgn0051875" "CG31875" "INTRON" "0"
保存snpA矩阵
save(snpA, file="./annotation/snpA_W2.Rdata")
基于序列本体的SNP集
删除无注释的snp。summary(snpA)显示一些snp存在NA 's。
snpA <- na.omit(snpA)
通过拆分snpA数据框架的第3列(“SeqOnt”)，创建序列本体集(seqSets)列表，使列表中的每个元素(即序列本体术语)都包含相关的snp。
seqSets <- split( rownames(snpA),f=as.factor( snpA[,3]) )
str(seqSets, vec.len=3)
## List of 14
## $ DOWNSTREAM : chr [1:258237] "2L_10008898" "2L_10008930" "2L_10008944" ...
## $ EXON : chr [1:11020] "2L_10506345" "2L_10506362" "2L_10506377" ...
## $ INTRON : chr [1:884290] "2L_10000016" "2L_10000033" "2L_10000089" ...
## $ NON_SYNONYMOUS_CODING: chr [1:52142] "2L_10000016" "2L_10000089" "2L_10000135" ...
## $ NON_SYNONYMOUS_START : chr [1:14] "2L_11591491" "2L_7857121" "2R_12186085" ...
## $ START_GAINED : chr [1:4090] "2L_10060881" "2L_10060894" "2L_10060983" ...
## $ START_LOST : chr [1:35] "2L_12111149" "2L_15895897" "2L_1791822" ...
## $ STOP_GAINED : chr [1:196] "2L_10236254" "2L_10363747" "2L_10527805" ...
## $ STOP_LOST : chr [1:71] "2L_10013272" "2L_10279121" "2L_11519730" ...
## $ SYNONYMOUS_CODING : chr [1:178949] "2L_10000033" "2L_10000487" "2L_10000538" ...
## $ SYNONYMOUS_STOP : chr [1:203] "2L_10055872" "2L_1166360" "2L_1199578" ...
## $ UPSTREAM : chr [1:258001] "2L_10002466" "2L_10002474" "2L_10002491" ...
## $ UTR_3_PRIME : chr [1:49454] "2L_10008023" "2L_10008437" "2L_10008502" ...
## $ UTR_5_PRIME : chr [1:24115] "2L_10002418" "2L_10002429" "2L_10003939" ...
save(seqSets, file="./annotation/seqSets2.Rdata")
基于基因id的SNP集
创建逻辑值inGene, inGene500bp and inGene1kbp对于基因距离0,500和1000碱基对(bp)的snp。
inGene <- as.numeric(snpA[,4])==0
inGene500bp <- as.numeric(snpA[,4])<=500
inGene1kbp <- as.numeric(snpA[,4])<=1000
通过拆分snpA数据框的第1列(FBgn)创建fb基因集(fbset)，使列表中的每个元素(即fb基因)都由与flybase基因相关的snp组成，且距离转录起始位点0 bp。
fbSets <- split( rownames(snpA[inGene,]),f=as.factor(snpA[inGene,1]) )
str(fbSets[1:5], vec.len=3)
## List of 5
## $ FBgn0000008: chr [1:853] "2R_18024616" "2R_18024682" "2R_18024761" ...
## $ FBgn0000014: chr [1:297] "3R_12633508" "3R_12633542" "3R_12633571" ...
## $ FBgn0000015: chr [1:450] "3R_12753398" "3R_12753600" "3R_12753792" ...
## $ FBgn0000017: chr [1:286] "3L_16615697" "3L_16615763" "3L_16616747" ...
## $ FBgn0000018: chr [1:21] "2L_10973623" "2L_10973738" "2L_10973758" ...
save(fbSets, file="./annotation/fbSets2.Rdata")
像以前一样使用分裂功能根据fb基因id对snp进行分组。本例包括转录起始位点SNPs 500 (fbSets500bp)和1000 (fbSets1kbp) bp。
fbSets500bp <- split( rownames(snpA[inGene500bp,]),f=as.factor(snpA[inGene500bp,1]) )
fbSets1kbp <- split( rownames(snpA[inGene1kbp,]),f=as.factor(snpA[inGene1kbp,1]) )
save(fbSets500bp, file="./annotation/fbSets2_500bp.Rdata")
save(fbSets1kbp, file="./annotation/fbSets2_1kbp.Rdata")
基于染色体的SNP集
创建染色体组。使用sapply将W的每个列名(例如2L_10000016)拆分下划线“_”为2L 10000016。结果行名的第一个元素(这里是[1]，例如“2L”)被保存在一个向量chr中。unique(chr)表示不同的染色体。
load(file="./genotypes/dgrp2_W2.Rdata")
chr <- sapply(colnames(W), function(x){ strsplit(x,split="_")[[1]][1] } )
unique(chr)
## [1] "2L" "2R" "3L" "3R" "4" "X"
chrSets <- unstack( colnames(W), colnames(W) ~ as.factor( chr ))
str(chrSets)
## List of 6
## $ 2L: chr [1:406577] "2L_5317" "2L_5372" "2L_5390" "2L_5403" ...
## $ 2R: chr [1:327967] "2R_10037" "2R_10468" "2R_10959" "2R_12079" ...
## $ 3L: chr [1:390711] "3L_39998" "3L_40145" "3L_40202" "3L_40635" ...
## $ 3R: chr [1:368096] "3R_1339" "3R_1651" "3R_2158" "3R_2318" ...
## $ 4 : chr [1:2686] "4_61790" "4_62622" "4_62905" "4_62908" ...
## $ X : chr [1:229718] "X_19380" "X_19797" "X_20390" "X_20491" ...
save(chrSets, file="./annotation/chrSets2.Rdata")
基于基因本体论(GO)的SNP集
生物导入包包含黑腹果蝇全基因组注释。注释主要基于使用Entrez基因标识符的映射。安装生物注释包，有关生物导子包装的更多信息，请访问:https://bioconductor.org/packages/release/data/annotation/html/org.Dm.eg.db.html.
source("https://bioconductor.org/biocLite.R")
biocLite("org.Dm.eg.db")
library(org.Dm.eg.db, verbose = FALSE, quietly = TRUE)
从转录起始位点将基因id载于0 bp处。
load(file="./annotation/fbSets2.Rdata")
有关基因本体论的信息可以在这里找到:http://www.geneontology.org/.
因为我们的数据集由fb基因id组成，所以fb基因id必须链接到entrez基因id，才能链接到GO项。
链接fb基因id到Entrez基因id。
fb2eg <- org.Dm.egFLYBASE2EG
mapped_genes <- mappedkeys(fb2eg)
fb2eg <- as.list(fb2eg[mapped_genes])
str(fb2eg[1:5])
## List of 5
## $ FBgn0040373: chr "30970"
## $ FBgn0040372: chr "30971"
## $ FBgn0261446: chr "30972"
## $ FBgn0000316: chr "30973"
## $ FBgn0005427: chr "30975"
链接Entrez基因id到fb基因id。
eg2fb <- org.Dm.egFLYBASE
mapped_genes <- mappedkeys(eg2fb)
eg2fb <- as.list(eg2fb[mapped_genes])
str(eg2fb[1:5])
## List of 5
## $ 10178776: chr "FBgn0261702"
## $ 10178777: chr "FBgn0262024"
## $ 10178779: chr "FBgn0058354"
## $ 10178780: chr "FBgn0053929"
## $ 10178781: chr "FBgn0262141"
Link GO ids to Entrez Gene ids
go2eg <- as.list(org.Dm.egGO2EG)
str(go2eg[1:5])
## List of 5
## $ GO:0000001: Named chr [1:2] "36309" "53567"
## ..- attr(*, "names")= chr [1:2] "IMP" "IMP"
## $ GO:0000002: Named chr [1:8] "33567" "34307" "35236" "38046" ...
## ..- attr(*, "names")= chr [1:8] "IBA" "IMP" "IMP" "IBA" ...
## $ GO:0000003: Named chr [1:5] "31153" "31309" "36619" "36619" ...
## ..- attr(*, "names")= chr [1:5] "IBA" "NAS" "IMP" "NAS" ...
## $ GO:0000011: Named chr "34110"
## ..- attr(*, "names")= chr "IBA"
## $ GO:0000012: Named chr [1:5] "31451" "33530" "41872" "42322" ...
## ..- attr(*, "names")= chr [1:5] "IEA" "IBA" "ISS" "IBA" ...
链接fb基因id到GO id。
go2fb <- lapply(go2eg,function(x){
fb <- na.omit(unlist(eg2fb[x]))
m <- match(fb,names(fbSets))
fb <- fb[!is.na(m)]
fb <- fb[!duplicated(fb)]
return(fb)
})
str(go2fb[1:5], vec.len = 3)
## List of 5
## $ GO:0000001: Named chr "FBgn0261618"
## ..- attr(*, "names")= chr "53567"
## $ GO:0000002: Named chr [1:7] "FBgn0031540" "FBgn0032154" "FBgn0040268" ...
## ..- attr(*, "names")= chr [1:7] "33567" "34307" "35236" ...
## $ GO:0000003: Named chr [1:4] "FBgn0023509" "FBgn0000479" "FBgn0003742" ...
## ..- attr(*, "names")= chr [1:4] "31153" "31309" "36619" ...
## $ GO:0000011: Named chr "FBgn0261064"
## ..- attr(*, "names")= chr "34110"
## $ GO:0000012: Named chr [1:4] "FBgn0026751" "FBgn0260817" "FBgn0026737" ...
## ..- attr(*, "names")= chr [1:4] "31451" "33530" "41872" ...
只包括包含5个或更多fb基因的GO。保存结果GO SNP集。
go2fb <- go2fb[sapply(go2fb,length)>=5]
save(go2fb,file="./annotation/go2fb.Rdata")
创建GO集:链接到GO项目的SNPs。
goSets <- lapply(go2fb,function(x){ unique(unlist(fbSets[x])) })
goSets <- goSets[sapply(goSets,length)>199]
str(goSets[1:5])
## List of 5
## $ GO:0000022: chr [1:315] "2R_2636861" "2R_2637023" "2R_2637293" "2R_2637711" ...
## $ GO:0000027: chr [1:644] "X_1771396" "X_1772294" "X_1772337" "X_1772796" ...
## $ GO:0000028: chr [1:283] "X_1376457" "X_1376829" "X_1376835" "X_9448856" ...
## $ GO:0000045: chr [1:1602] "X_7208847" "X_7209453" "X_7209470" "X_7209504" ...
## $ GO:0000060: chr [1:1816] "X_2056263" "X_2056399" "X_2056417" "X_2056810" ...
save(goSets,file="./annotation/goSets2.Rdata")
1764
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
