搜索
查看: 3024|回复: 0

[software] SNPrelate包

[复制链接]

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-7-31 15:11:15 | 显示全部楼层 |阅读模式
跟这运行的时候注意里面有些东西要换
[AppleScript] 纯文本查看 复制代码
---

GDS也被R / Bioconductor包GWASTools用作其数据存储格式之一2,3。 GWASTools提供许多功能,用于GWAS的质量控制和分析,包括SNP或扫描统计,批次质量,染色体异常,关联检测等。
```{r setup, include=FALSE}
library(gdsfmt)
library(SNPRelate)
```

 
典型的GDS文件
```{r GDS}
snpgdsSummary(snpgdsExampleFileName())
```

snpgdsExampleFileName()返回在SNPRelate中作为示例使用的GDS文件的文件名,它是HapMap项目的数据的一个子集,样本由约翰霍普金斯大学的继承性疾病研究中心(CIDR)进行基因分型, 麻省理工学院和哈佛大学广泛研究所(Broad)。 snpgdsSummary()总结存储在GDS文件中的基因型。 “个人主要模式”表示在列出下一个人的SNP之前列出个人的所有SNP等。相反,“SNP主要模式”表示列出第一个SNP的所有个体之前列出第一个SNP的所有个体等 有时候“SNP主要模式”比“个人主要模式”更具有计算效率。 例如,遗传协方差矩阵的计算由SNP处理基因型数据SNP,然后“SNP主模式”应该更有效。
```{r  , echo=FALSE}
# Open a GDS file
(genofile <- snpgdsOpen(snpgdsExampleFileName()))
```

Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot.
```{r  , echo=FALSE}
# Get the attributes of chromosome coding
get.attr.gdsn(index.gdsn(genofile, "snp.chromosome"))
```

autosome.start是常数形式的起始数字代码,autosome.end是常数的最后一个数字代码。 put.attr.gdsn()可用于添加新属性或修改现有属性。

可变基因型中存在四个可能的值:0,1,2和3.对于双等位基因SNP位点,“0”表示两个B等位基因,“1”表示一个A等位基因和一个B等位基因,“2”表示 两个A等位基因,“3”是缺失的基因型。 对于多等位基因,它是参考等位基因的计数(3表示无呼叫)。 “Bit2”表示每个字节最多编码四个SNP基因型,因为一个字节由8位组成。

```{r  , echo=FALSE}
# Take out genotype data for the first 3 samples and the first 5 SNPs
(g <- read.gdsn(index.gdsn(genofile, "genotype"), start=c(1,1), count=c(5,3)))
```
或者取出具有样本和SNP ID的基因型数据,并返回0,1,2和NA(4个可替换为NA)的四个可能值:
```{r  , echo=FALSE}
# Get the attribute of genotype
get.attr.gdsn(index.gdsn(genofile, "genotype"))
```
返回的值可以是“snp.order”或“sample.order”,分别表示个别主模式(snp是第一维)和SNP主模式(样本是第一维)
```{r  , echo=FALSE}
# Take out snp.id
head(read.gdsn(index.gdsn(genofile, "snp.id")))
```

```{r  , echo=FALSE}
#  snp.rs.id
head(read.gdsn(index.gdsn(genofile, "snp.rs.id")))
```
有两个额外的和可选的变量:

snp.rs.id,一个可能不是唯一的参考SNP ID的字符串。
snappallele,没有必要进行分析,但是合并来自不同平台的基因型是必要的。  snp.allele的形式是“A等位基因/ B等位基因”,如“T / G”,T是A等位基因,G是B等位基因。
样本注释的信息可以通过相同的函数read.gdsn()获得。 例如,opulation information。 “VStr8”表示字符型变量。
```{r  , echo=FALSE}
# Read population information
pop <- read.gdsn(index.gdsn(genofile, path="sample.annot/pop.group"))
table(pop)
```
```{r  , echo=FALSE}
# Close the GDS file
snpgdsClose(genofile)
```

###建立一个snp的GDS文件
snpgdsCreateGeno
函数snpgdsCreateGeno()可用于创建GDS文件。 第一个参数应该是SNP基因型的数字矩阵。 输入基因型矩阵中存在可能的值:0,1,2和其他值。 “0”表示两个B等位基因,“1”表示一个A等位基因和一个B等位基因,“2”表示两个A等位基因,其他值表示缺失的基因型。 SNP矩阵可以是nsample×nsnpnsample×nsnp(snpfirstdim = FALSE,snpgdsCreateGeno()中的参数)或nsnp×nsamplensnp×nsample(snpfirstdim = TRUE)。
```{r  , echo=FALSE}
# Load data
data(hapmap_geno)

# Create a gds file
snpgdsCreateGeno("test.gds", genmat = hapmap_geno$genotype,
    sample.id = hapmap_geno$sample.id, snp.id = hapmap_geno$snp.id,
    snp.chromosome = hapmap_geno$snp.chromosome,
    snp.position = hapmap_geno$snp.position,
    snp.allele = hapmap_geno$snp.allele, snpfirstdim=TRUE)

# Open the GDS file
(genofile <- snpgdsOpen("test.gds"))
```

```{r  , echo=FALSE}
  # Close the GDS file
snpgdsClose(genofile)
```
#软件包gdsfmt中函数的使用
在以下代码中,函数createfn.gds(),add.gdsn(),put.attr.gdsn(),write.gdsn()和index.gdsn()在包gdsfmt中定义:
```{r  , echo=FALSE}
# Create a new GDS file
newfile <- createfn.gds("your_gds_file.gds")

# add a flag
put.attr.gdsn(newfile$root, "FileFormat", "SNP_ARRAY")

# Add variables
add.gdsn(newfile, "sample.id", sample.id)
add.gdsn(newfile, "snp.id", snp.id)
add.gdsn(newfile, "snp.chromosome", snp.chromosome)
add.gdsn(newfile, "snp.position", snp.position)
add.gdsn(newfile, "snp.allele", c("A/G", "T/C", ...))

#####################################################################
# Create a snp-by-sample genotype matrix

# Add genotypes
var.geno <- add.gdsn(newfile, "genotype",
    valdim=c(length(snp.id), length(sample.id)), storage="bit2")

# Indicate the SNP matrix is snp-by-sample
put.attr.gdsn(var.geno, "snp.order")

# Write SNPs into the file sample by sample
for (i in 1:length(sample.id))
{
    g <- ...
    write.gdsn(var.geno, g, start=c(1,i), count=c(-1,1))
}

#####################################################################
# OR, create a sample-by-snp genotype matrix

# Add genotypes
var.geno <- add.gdsn(newfile, "genotype",
    valdim=c(length(sample.id), length(snp.id)), storage="bit2")

# Indicate the SNP matrix is sample-by-snp
put.attr.gdsn(var.geno, "sample.order")

# Write SNPs into the file sample by sample
for (i in 1:length(snp.id))
{
    g <- ...
    write.gdsn(var.geno, g, start=c(1,i), count=c(-1,1))
}

# Get a description of chromosome codes
#   allowing to define a new chromosome code, e.g., snpgdsOption(Z=27)
option <- snpgdsOption()
var.chr <- index.gdsn(newfile, "snp.chromosome")
put.attr.gdsn(var.chr, "autosome.start", option$autosome.start)
put.attr.gdsn(var.chr, "autosome.end", option$autosome.end)
for (i in 1:length(option$chromosome.code))
{
    put.attr.gdsn(var.chr, names(option$chromosome.code)[i],
        option$chromosome.code[[i]])
}

# Add your sample annotation
samp.annot <- data.frame(sex = c("male", "male", "female", ...),
    pop.group = c("CEU", "CEU", "JPT", ...), ...)
add.gdsn(newfile, "sample.annot", samp.annot)

# Add your SNP annotation
snp.annot <- data.frame(pass=c(TRUE, TRUE, FALSE, FALSE, TRUE, ...), ...)
add.gdsn(newfile, "snp.annot", snp.annot)

# Close the GDS file
closefn.gds(newfile)
```
 从PLINK文本/二进制文件格式转换
SNPRelate包提供了一个函数snpgdsPED2GDS()和snpgdsBED2GDS(),用于将PLINK文本/二进制文件转换为GDS文件:
```{r  , echo=FALSE}
# The PLINK BED file, using the example in the SNPRelate package
bed.fn <- system.file("extdata", "plinkhapmap.bed.gz", package="SNPRelate")
fam.fn <- system.file("extdata", "plinkhapmap.fam.gz", package="SNPRelate")
bim.fn <- system.file("extdata", "plinkhapmap.bim.gz", package="SNPRelate")
```

```{r  , echo=FALSE}
#bed.fn <- "C:/your_folder/your_plink_file.bed"
#fam.fn <- "C:/your_folder/your_plink_file.fam"
#bim.fn <- "C:/your_folder/your_plink_file.bim"
# Convert
snpgdsBED2GDS(bed.fn, fam.fn, bim.fn, "test.gds")
```

```{r  , echo=FALSE}
# Summary
snpgdsSummary("test.gds")
```
从VCF文件格式转换
SNPRelate包提供了一个函数snpgdsVCF2GDS()来重新格式化一个VCF文件。 从VCF文件中提取标记物用于下游分析有两种选择:1.仅提取和存储参考等位基因的双重SNP的剂量2.提取和存储所有变异位点的参考等位基因的剂量,包括双等位基因 SNP,多等位基因SNP,indel和结构变体。
```{r  , echo=FALSE}
# The VCF file, using the example in the SNPRelate package
vcf.fn <- system.file("extdata", "sequence.vcf", package="SNPRelate")
# Reformat
snpgdsVCF2GDS(vcf.fn, "test.gds", method="biallelic.only")
# Summary
snpgdsSummary("test.gds")
```
##数据分析
我们开发了gdsfmt和SNPRelate(用于多核对称多处理计算机架构的高性能计算R包),以加速GWAS中的两个关键计算:主成分分析(PCA)和使用身份下降(IBD)度量的相关性分析。
```{r  , echo=FALSE}
# Open the GDS file
genofile <- snpgdsOpen(snpgdsExampleFileName())
# Get population information
#   or pop_code <- scan("pop.txt", what=character())
#   if it is stored in a text file "pop.txt"
pop_code <- read.gdsn(index.gdsn(genofile, path="sample.annot/pop.group"))

table(pop_code)
# Display the first six values
head(pop_code)
```
基于LD的SNP修剪
建议使用剪切的一组SNP,它们彼此大致连接平衡,以避免SNP簇在主成分分析和相关性分析中的强烈影响。
```{r  , echo=FALSE}
set.seed(1000)

# Try different LD thresholds for sensitivity analysis
snpset <- snpgdsLDpruning(genofile, ld.threshold=0.2)
names(snpset)
head(snpset$chr1)  # snp.id
# Get all selected snp id
snpset.id <- unlist(snpset)
```
主成分分析(PCA)
用于PCA的SNPRelate中的功能包括从基因型计算遗传协方差矩阵,计算每个SNP的样本载荷和基因型之间的相关系数,计算SNP特征向量(载荷),以及从指定的SNP特征向量估计新数据集的样本加载。
```{r  , echo=FALSE}
# Run PCA
pca <- snpgdsPCA(genofile, snp.id=snpset.id, num.thread=2)
```
下面的代码显示了如何计算最主要组成部分所占的变化百分比。 很明显,前两个特征向量在人口中占有最大的差异百分比,尽管总方差仍然少于总数的四分之一。
```{r  , echo=FALSE}
# variance proportion (%)
pc.percent <- pca$varprop*100
head(round(pc.percent, 2))
#In the case of no prior population information,

# make a data.frame
tab <- data.frame(sample.id = pca$sample.id,
    EV1 = pca$eigenvect[,1],    # the first eigenvector
    EV2 = pca$eigenvect[,2],    # the second eigenvector
    stringsAsFactors = FALSE)
head(tab)
# Draw
plot(tab$EV2, tab$EV1, xlab="eigenvector 2", ylab="eigenvector 1")

```
If there are population information,
```{r  , echo=FALSE}
# Get sample id
sample.id <- read.gdsn(index.gdsn(genofile, "sample.id"))

# Get population information
#   or pop_code <- scan("pop.txt", what=character())
#   if it is stored in a text file "pop.txt"
pop_code <- read.gdsn(index.gdsn(genofile, "sample.annot/pop.group"))

# assume the order of sample IDs is as the same as population codes
head(cbind(sample.id, pop_code))

# Make a data.frame
tab <- data.frame(sample.id = pca$sample.id,
    pop = factor(pop_code)[match(pca$sample.id, sample.id)],
    EV1 = pca$eigenvect[,1],    # the first eigenvector
    EV2 = pca$eigenvect[,2],    # the second eigenvector
    stringsAsFactors = FALSE)
head(tab)

# Draw
plot(tab$EV2, tab$EV1, col=as.integer(tab$pop), xlab="eigenvector 2", ylab="eigenvector 1")
legend("bottomright", legend=levels(tab$pop), pch="o", col=1:nlevels(tab$pop))
```
绘制前四个PC的主要组件对:
```{r  , echo=FALSE}
lbls <- paste("PC", 1:4, "\n", format(pc.percent[1:4], digits=2), "%", sep="")
```
主要组成部分的平行坐标图:
```{r  , echo=FALSE}
library(MASS)

datpop <- factor(pop_code)[match(pca$sample.id, sample.id)]
```
计算本征实体和SNP基因型之间的SNP相关性:
```{r  , echo=FALSE}
# Get chromosome index
chr <- read.gdsn(index.gdsn(genofile, "snp.chromosome"))
CORR <- snpgdsPCACorr(pca, genofile, eig.which=1:4)
savepar <- par(mfrow=c(2,1), mai=c(0.45, 0.55, 0.1, 0.25))
for (i in 1:2)
{
    plot(abs(CORR$snpcorr[i,]), ylim=c(0,1), xlab="", ylab=paste("PC", i),
        col=chr, pch="+")
}
par(savepar)
parcoord(pca$eigenvect[,1:16], col=datpop)
pairs(pca$eigenvect[,1:4], col=tab$pop, labels=lbls)

```
#Fst估计

给定两个或更多的群体 ,Fst可以通过Weir&Cockerham(1984)的方法来估计。
```{r  , echo=FALSE}
# Get sample id
sample.id <- read.gdsn(index.gdsn(genofile, "sample.id"))

# Get population information
#   or pop_code <- scan("pop.txt", what=character())
#   if it is stored in a text file "pop.txt"
pop_code <- read.gdsn(index.gdsn(genofile, "sample.annot/pop.group"))

# Two populations: HCB and JPT
flag <- pop_code %in% c("HCB", "JPT")
samp.sel <- sample.id[flag]
pop.sel <- pop_code[flag]
```
相关性分析

对于相关性分析,SNPRelate中的身份下降(IBD)估计可以通过时间方法(MoM)(Purcell等,2007)或最大似然估计(MLE)(Milligan,2003; Choi et al ,,2009)。 对于这两种方法,优先使用LD修剪的SNP集合。
```{r  , echo=FALSE}
v <- snpgdsFst(genofile, sample.id=samp.sel, population=as.factor(pop.sel),
    method="W&C84")
#估计IBD使用PLINK方法(MoM)
# Estimate IBD coefficients
ibd <- snpgdsIBDMoM(genofile, sample.id=YRI.id, snp.id=snpset.id,
    maf=0.05, missing.rate=0.05, num.thread=2)
# Make a data.frame
ibd.coeff <- snpgdsIBDSelection(ibd)
head(ibd.coeff)
plot(ibd.coeff$k0, ibd.coeff$k1, xlim=c(0,1), ylim=c(0,1),
    xlab="k0", ylab="k1", main="YRI samples (MoM)")
lines(c(0,1), c(1,0), col="red", lty=2)
#使用最大似然估计(MLE)估计IBD
# Estimate IBD coefficients
set.seed(100)
snp.id <- sample(snpset.id, 1500)  # random 1500 SNPs
ibd <- snpgdsIBDMLE(genofile, sample.id=YRI.id, snp.id=snp.id,
    maf=0.05, missing.rate=0.05, num.thread=2)
# Make a data.frame
ibd.coeff <- snpgdsIBDSelection(ibd)

plot(ibd.coeff$k0, ibd.coeff$k1, xlim=c(0,1), ylim=c(0,1),
    xlab="k0", ylab="k1", main="YRI samples (MLE)")
lines(c(0,1), c(1,0), col="red", lty=2)
#Relationship inference Using KING method of moments
#家庭间和家庭间的关系可以通过在人口分层的情况下采用KING-robust方法来推断。
# Incorporate with pedigree information
family.id <- read.gdsn(index.gdsn(genofile, "sample.annot/family.id"))
family.id <- family.id[match(YRI.id, sample.id)]
table(family.id)
ibd.robust <- snpgdsIBDKING(genofile, sample.id=YRI.id,
    family.id=family.id, num.thread=2)
names(ibd.robust)
# Pairs of individuals
dat <- snpgdsIBDSelection(ibd.robust)
head(dat)
plot(dat$IBS0, dat$kinship, xlab="Proportion of Zero IBS",
    ylab="Estimated Kinship Coefficient (KING-robust)")
#Identity-By-State Analysis
#对于nn个个体,snpgdsIBS()可用于创建全基因组平均IBS成对特征的n×nn×n矩阵:
ibs <- snpgdsIBS(genofile, num.thread=2)
# individulas in the same population are clustered together
pop.idx <- order(pop_code)

image(ibs$ibs[pop.idx, pop.idx], col=terrain.colors(16))
#对全基因组IBS成对距离的n×nn×n矩阵执行多维缩放分析:
loc <- cmdscale(1 - ibs$ibs, k = 2)
x <- loc[, 1]; y <- loc[, 2]
race <- as.factor(pop_code)

plot(x, y, col=race, xlab = "", ylab = "",
    main = "Multidimensional Scaling Analysis (IBS)")
legend("topleft", legend=levels(race), pch="o", text.col=1:nlevels(race))
#对全基因组IBS成对距离的n×nn×n矩阵进行聚类分析,并通过置换评分确定组:
set.seed(100)
ibs.hc <- snpgdsHCluster(snpgdsIBS(genofile, num.thread=2))
# Determine groups of individuals automatically
rv <- snpgdsCutTree(ibs.hc)
plot(rv$dendrogram, leaflab="none", main="HapMap Phase II")
table(rv$samp.group)
# Determine groups of individuals by population information
rv2 <- snpgdsCutTree(ibs.hc, samp.group=as.factor(pop_code))
plot(rv2$dendrogram, leaflab="none", main="HapMap Phase II")
legend("topright", legend=levels(race), col=1:nlevels(race), pch=19, ncol=4)
# Close the GDS file
snpgdsClose(genofile)
```

回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-10-14 22:22 , Processed in 0.039932 second(s), 24 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.