芬兰网站后缀,开源社区源码,企业关键词排名优化网址,学生个人网页制作图片全基因组关联#xff08;GWA#xff09;研究扫描整个物种基因组#xff0c;寻找多达数百万个SNPs与特定感兴趣特征之间的关联。值得注意的是#xff0c;感兴趣的性状实际上可以是归因于群体的任何类型的表型#xff0c;无论是定性的#xff08;例如疾病状态#xff09;还…全基因组关联GWA研究扫描整个物种基因组寻找多达数百万个SNPs与特定感兴趣特征之间的关联。值得注意的是感兴趣的性状实际上可以是归因于群体的任何类型的表型无论是定性的例如疾病状态还是定量的例如身高。本质上给定p个SNP和n个样本或个体GWA分析将拟合p个独立的单变量线性模型每个模型基于n个样本使用每个SNP的基因型作为感兴趣特征的预测因子。每个P检验中的关联显著性P值由相应SNP的系数估计确定从技术上讲关联的显著性为。请注意由于这些测试是独立的且数量众多因此在设置并行GWA分析时有很大的计算优势我们将很快进行。非常合理的是有必要使用多种假设检验方法如Bonferroni、Benjamini Hochberg或错误发现率FDR来调整所得的P值。GWA研究现在在许多不同物种的遗传学中很常见。 读取数据
让我们首先从三个民族中的每一个导入PLINK转换的.ded、.fam和.bim Illumina文件。我们将使用snpStats包中的函数read.plink并在本教程的其余部分中处理生成的对象。此函数读取.bed、.fam和.bim并创建一个包含三个元素的列表——$genesis、$fam和$map。第一个包含从所有样本中确定的所有SNP第二个包含关于谱系和性别的信息第三个包含SNP的基因组坐标。目前我们共有323人110名中国人、105名印度人和108名马来人和2527458个SNPs。接下来我们将把$genetype中存储的Illumina SNP ID更改为更传统的rs ID这将使我们能够在USCS基因组浏览器中放大候选SNP周围的基因组区域。我准备了一个表来建立两者之间的对应关系这样我们就可以轻松地切换ID。
library(snpStats)
load(conversionTable.RData)pathM - paste(Genomics/108Malay_2527458snps, c(.bed, .bim, .fam), sep )
SNP_M - read.plink(pathM[1], pathM[2], pathM[3])pathI - paste(Genomics/105Indian_2527458snps, c(.bed, .bim, .fam), sep )
SNP_I - read.plink(pathI[1], pathI[2], pathI[3])pathC - paste(Genomics/110Chinese_2527458snps, c(.bed, .bim, .fam), sep )
SNP_C - read.plink(pathC[1], pathC[2], pathC[3])# Merge the three SNP datasets
SNP - rbind(SNP_M$genotypes, SNP_I$genotypes, SNP_C$genotypes)# Take one bim map (all 3 maps are based on the same ordered set of SNPs)
map - SNP_M$map
colnames(map) - c(chr, SNP, gen.dist, position, A1, A2)# Rename SNPs present in the conversion table into rs IDs
mappedSNPs - intersect(map$SNP, names(conversionTable))
newIDs - conversionTable[match(map$SNP[map$SNP %in% mappedSNPs], names(conversionTable))]
map$SNP[rownames(map) %in% mappedSNPs] - newIDs
接下来我们导入并合并三个脂质数据集存储为.txt并确定SNP和脂质数据集中都存在哪些样本。在以下分析中我们将使用两个平台中介绍的样本子集共计319个。最后我们创建一个列表genData存储合并的SNP数据$SNP这是三个$map中的一个因为它们都是相同的$map和合并的脂质数据$LIP。最后让我们为后续步骤保存RAM并在将genData保存到之后删除所有文件。RData文件。
# Load lipid datasets match SNP-Lipidomics samples
lipidsMalay - read.delim(Lipidomic/117Malay_282lipids.txt, row.names 1)
lipidsIndian - read.delim(Lipidomic/120Indian_282lipids.txt, row.names 1)
lipidsChinese - read.delim(Lipidomic/122Chinese_282lipids.txt, row.names 1)
all(Reduce(intersect, list(colnames(lipidsMalay),
colnames(lipidsIndian),
colnames(lipidsChinese))) colnames(lipidsMalay)) # TRUE
lip - rbind(lipidsMalay, lipidsIndian, lipidsChinese)matchingSamples - intersect(rownames(lip), rownames(SNP))
SNP - SNP[matchingSamples,]
lip - lip[matchingSamples,]genData - list(SNP SNP, MAP map, LIP lip)
save(genData, file PhenoGenoMap.RData)# Clear memory
rm(list ls()) 预处理
让我们重新加载genData并对其进行清理。简而言之数据的预处理包括 丢弃呼叫率1或MAF0.1的SNPs丢弃呼叫率100%、IBD亲缘系数0.1或近交系数|F|0.1的样本呼叫率是进行基因分型的SNPs或样本的比例。例如特定SNP样本的呼叫率为0.95意味着丢失了5%的值。因为我们有这么多的SNP我们可以通过对SNP和样本施加1的呼叫率阈值在$SNP矩阵中绝对没有遗漏值。如果你想放宽门槛并容忍缺失的值请记住你需要运行一个完整的程序来对这些值进行指责。Reed等人2015描述了一种基于PCA的插补方法该方法利用1000基因组计划作为代理以防您感兴趣。 次要等位基因频率MAF表示每个SNP的最不常见等位基因的比例。当然很难检测出与罕见变异的关联这就是为什么我们选择低MAF值。我读过的大多数GWA研究通常报告MAF阈值为0.05。在这里我选择了更严格的0.1因为我们有大量的数据而且由于所有这些都是为了说明目的我们希望进行快速的GWA分析。
library(snpStats)
library(doParallel)
library(SNPRelate)
library(GenABEL)
library(dplyr)
source(GWASfunction.R)
load(PhenoGenoMap.RData)# Use SNP call rate of 100%, MAF of 0.1 (very stringent)
maf - 0.1
callRate - 1
SNPstats - col.summary(genData$SNP)maf_call - with(SNPstats, MAF maf Call.rate callRate)
genData$SNP - genData$SNP[,maf_call]
genData$MAP - genData$MAP[maf_call,]
SNPstats - SNPstats[maf_call,]
接下来我们需要考虑表现出过度杂合性的样本——从技术上讲偏离Hardy-Weinberg平衡HWE # Sample call rate heterozygosity
callMat - !is.na(genData$SNP)
Sampstats - row.summary(genData$SNP)
hetExp - callMat %*% (2 * SNPstats$MAF * (1 - SNPstats$MAF)) # Hardy-Weinberg heterozygosity (expected)
hetObs - with(Sampstats, Heterozygosity * (ncol(genData$SNP)) * Call.rate)
Sampstats$hetF - 1-(hetObs/hetExp)
# Use sample call rate of 100%, het threshold of 0.1 (very stringent)
het - 0.1 # Set cutoff for inbreeding coefficient;
het_call - with(Sampstats, abs(hetF) het Call.rate 1)
genData$SNP - genData$SNP[het_call,]
genData$LIP - genData$LIP[het_call,]
最后我们将使用基于血统认同的亲属系数IBD来研究样本之间的相关性。请注意SNPRelate包中的这些函数需要GDS文件。出于这个原因我们首先需要将来自三个群体的.bed、.fam和.bim文件聚合到convertGDS中。函数snpgdsBED2GDS2创建这部分分析所需的GDS。为了确定样本对之间的亲缘关系系数我们将使用不相关SNPs的子集以便进行无偏估计。为此我们将使用连锁不平衡LD作为标记之间相关性的衡量标准。LD的范围从0到1其值越高两个SNP就越有可能共同分离并因此相关。在这里我们将利用LD0.2p~12000的SNPs子集来确定IBD亲缘系数。我用函数snpgdsLDpruning计算LD大约花了2个小时所以要耐心。接下来根据Reed等人2015我们将排除所有亲缘系数0.1的样本。 # LD and kinship coeff
ld - .2
kin - .1
snpgdsBED2GDS(bed.fn convertGDS.bed, bim.fn convertGDS.bim,
fam.fn convertGDS.fam, out.gdsfn myGDS, cvt.chr char)
genofile - snpgdsOpen(myGDS, readonly F)
gds.ids - read.gdsn(index.gdsn(genofile, sample.id))
gds.ids - sub(-1, , gds.ids)
add.gdsn(genofile, sample.id, gds.ids, replace T)
geno.sample.ids - rownames(genData$SNP)
# First filter for LD
snpSUB - snpgdsLDpruning(genofile, ld.threshold ld,
sample.id geno.sample.ids,
snp.id colnames(genData$SNP))
snpset.ibd - unlist(snpSUB, use.names F)
# And now filter for MoM
ibd - snpgdsIBDMoM(genofile, kinship T,
sample.id geno.sample.ids,
snp.id snpset.ibd,
num.thread 1)
ibdcoef - snpgdsIBDSelection(ibd)
ibdcoef - ibdcoef[ibdcoef$kinship kin,]# Filter samples out
related.samples - NULL
while (nrow(ibdcoef) 0) {
# count the number of occurrences of each and take the top one
sample.counts - arrange(count(c(ibdcoef$ID1, ibdcoef$ID2)), -freq)
rm.sample - sample.counts[1, x]
cat(Removing sample, as.character(rm.sample), too closely related to, sample.counts[1, freq],other samples.\n)# remove from ibdcoef and add to list
ibdcoef - ibdcoef[ibdcoef$ID1 ! rm.sample ibdcoef$ID2 ! rm.sample,]
related.samples - c(as.character(rm.sample), related.samples)
}
genData$SNP - genData$SNP[!(rownames(genData$SNP) %in% related.samples),]
genData$LIP - genData$LIP[!(rownames(genData$LIP) %in% related.samples),]
预处理后我们留下316个样本110个中国人、100个印度人和106个马来人其特征在于795668个SNP标记和282个脂质物种。请注意由于LD修剪过程是随机的因此您的样本量可能略有不同。 分析
主成分分析
既然我们已经完成了预处理那么使用主成分分析PCA检查基因型数据中最大的变异源并找出异常值或聚类模式可能是个好主意。因为我们使用的是S4对象所以我们将使用SNPRelate中的PCA函数snpgdsPCA。让我们绘制前两个主要组件PC。
# PCA
pca - snpgdsPCA(genofile, sample.id geno.sample.ids, snp.id snpset.ibd, num.thread 1)
pctab - data.frame(sample.id pca$sample.id,
PC1 pca$eigenvect[,1],
PC2 pca$eigenvect[,2],
stringsAsFactors F)origin - read.delim(countryOrigin.txt, sep \t)
origin - origin[match(pca$sample.id, origin$sample.id),]pcaCol - rep(rgb(0,0,0,.3), length(pca$sample.id)) # Set black for chinese
pcaCol[origin$Country I] - rgb(1,0,0,.3) # red for indian
pcaCol[origin$Country M] - rgb(0,.7,0,.3) # green for malaypng(PCApopulation.png, width 500, height 500)
plot(pctab$PC1, pctab$PC2, xlab PC1, ylab PC2, col pcaCol, pch 16)
abline(h 0, v 0, lty 2, col grey)
legend(top, legend c(Chinese, Indian, Malay), col 1:3, pch 16, bty n)
dev.off()
不出所料795668个SNP标记清楚地描绘了这三个群体。研究结果还表明与印度人相比中国人和马来人更接近彼此这一观察结果可以通过分层聚类等方法更好地解决。此外前两个PC没有发现明显的异常值。一切都很好可以继续进行GWA。
全基因组联合
最后抵抗力。从282种脂质中我选择了最熟悉的胆固醇作为感兴趣的特征之一。您可以完全自由地选择不同的脂质种类并继续操作。我们将使用Reed等人2015年提供的GWA函数并进行一些小的修改。我建议您打开脚本GWAS函数。R然后浏览一遍。这是一个优秀的、文档丰富的并行化实现。注意glm函数用于确定每个SNP与感兴趣的性状之间的关联的显著性。如果你想知道glm比lm更通用因为它可以进行高斯、泊松、二项和多项回归/分类这取决于你感兴趣的特征是如何分布的表型文件中的所有脂质都是高斯的。该GWA函数不会创建变量而是编写一个.txt汇总表列出每个SNP的系数估计\βt和相应的P值以及相应的基因组坐标。2012年年中我使用MacBook Pro8Gb2.9 GHz英特尔酷睿i7运行GWA功能花了大约1.5个小时。
# Choose trait for association analysis, use colnames(genData$LIP) for listing
# NOTE: Ignore the first column of genData$LIP (gender)
target - Cholesterolphenodata - data.frame(id rownames(genData$LIP),
phenotype scale(genData$LIP[,target]), stringsAsFactors F)# Conduct GWAS (will take a while)
start - Sys.time()
GWAA(genodata genData$SNP, phenodata phenodata, filename paste(target, .txt, sep ))
Sys.time() - start # benchmark
一旦完成我们可以使用所谓的曼哈顿图来可视化结果。我们所需要的只是加载在前一步中编写的.txt汇总表添加一个带有-\log_{10}P的列并根据所有SNP的基因组坐标绘制这些显著性估计。
# Manhattan plot
GWASout - read.table(paste(target, .txt, sep ), header T, colClasses c(character, rep(numeric,4)))
GWASout$type - rep(typed, nrow(GWASout))
GWASout$Neg_logP - -log10(GWASout$p.value)
GWASout - merge(GWASout, genData$MAP[,c(SNP, chr, position)])
GWASout - GWASout[order(GWASout$Neg_logP, decreasing T),]png(paste(target, .png, sep ), height 500,width 1000)
GWAS_Manhattan(GWASout)
dev.off()
此外一条完整的分别为虚线线表示1000000次分别为10000次测试中Bonferroni调整的\alpha0.05的水平。 我们发现共有四个SNP通过了“宽松”的Bonferroni阈值没有一个通过“硬”阈值。这些是SNPs rs7527051Chr.1、rs12140539与第一个Chr.1共定位、rs9509213Chr.13和rs2250402Chr.15。 在进行这四次点击之前将所得P值的分布与偶然预期的分布进行对比是有帮助的以确保没有混淆的系统性偏差。这很容易用分位数-分位数Q-Q图来解决。顾名思义它根据分位数比较两种分布。在目前的情况下我们想将我们的t统计数据的分布与偶然获得的分布进行对比。如果两个分布的形状相同Q-Q图将显示一条xy线。因此来自可靠GWA分析的Q-Q图将显示一条xy线其中只有很少的偏差值表明存在关联。否则如果线向上或向下移动GWA分析可能会受到混杂因素的影响。我们将使用GenABEL包中的函数estlambda绘制Q-Q图。 # QQ plot using GenABEL estlambda function
png(paste(target, _QQplot.png, sep ), width 500, height 500)
lambda - estlambda(GWASout$t.value**2, plot T, method median)
dev.off() 由此得到的Q-Q图清楚地描绘了一条与xy黑色重叠的趋势线λ0.99红色和右尾部的轻微偏差因此我们可以对我们的结果充满信心。 候选标记的功能
见解最后我们将尝试使用USCS基因组浏览器输入基因组浏览器在文本框中插入SNP ID输入并缩小通过搜索其附近的基因来找到这四个候选SNP的功能相关性。我发现rs9509213Chr.13正好落在CRYL1晶体λ1内含子序列上使其成为后续研究的有趣候选者。 SNP rs9509213在图的底部显示为一个黑色文本框其坐标由一条垂直的黄线突出显示。CRYL1基因模型显示在染色体模型最上面下方的顶部。 有趣的是最近有一份出版物发现“POE、HP和CRYL1都与阿尔茨海默病有关阿尔茨海默病的病理学涉及脂质和胆固醇途径。”。 最后我想谈谈对候选基因进行实验验证的必要性以验证由此产生的SNP特征关联。在这种特殊的情况下可以测试在小鼠的胆固醇水平中过表达或敲除CRYL1的直系同源物的效果。
