文章来源：Local PCA Shows How the Effect of Population Structure Differs Along the Genome
群体结构是衡量大型基因组数据集中个体之间相关性度量的系统化模式，通常使用降维技术，如主成分分析（PCA）来发现和可视化这些模式。平均关联性是特定于位点的系谱树之间关系的平均值，这种关联性可能会受到相关选择和其他因素的强烈影响，如结构变异等。通过使用区域PCA(localPCA)在关联模式中描述这种小规模的异质性，用于发现影响仅在几个兆碱基的种群结构变化。
1. 安装
install.packages("devtools")
install.packages("data.table")
devtools::install_github("petrelharp/local_pca/lostruct")
library(lostruct)
2. 所需软件
（1）读取VCF或BCF文件：bcftools
（2）实现测试分析：下载模版
3. 运行
（1）使用 eigen_windows() 计算局部PCA坐标
函数 eigen_windows() 将数据存储在矩阵中，每个SNP一行，每个样本一列（因此x [ i，j ]是样本 j 在位置 i 处的等位基因数）。
#两种方式从标准格式(tped和vcf)中读取数据
snps <- read_tped("mydata.tped")
snps <- read_vcf("mydata.vcf")
#或者
bcftools convert -O b mydata.vcf > mydata.bcf
bcftools index mydata.bcf
snps <- vcf_windower("my_data.bcf",size=1e3,type='snp')
snps(10)
#计算localPCA坐标
pcs <- eigen_windows(snps,k=2)
#计算距离矩阵
pcdist <- pc_dist(pcs,npc=2)
这将提供一个矩阵pcs，行是每个窗口的前两个特征值和特征向量，以及pcdist，即这些窗口之间的成对距离矩阵。
（2）使用 cmdscale() 可视化MDS
#可视化代码示例
a<-read.table("quick_method_pairwise_distance_between_win_100_chr1")
k=which(is.na(a),arr.ind=TRUE)
a[k]=0
a=a+t(a)
a <- as.matrix(a)
a=sqrt(a)
fit2d<-cmdscale(a,eig=TRUE, k=2) # MDS_2D#
x <- fit2d$points[,1]
y <- fit2d$points[,2]
pdf(file="MDS_2D_win100_chr1.pdf")
layout(1:2)
par(mar=c(4,3,1,1)+0.1)
plot(x, y, xlab="Coordinate 1", ylab="Coordinate 2", main="MDS_2D_chr1", col=rainbow(2*nrow(a)))
plot(1:nrow(a),col=rainbow(2*nrow(a)))
dev.off()
fit1d<-cmdscale(a,eig=TRUE, k=1) # MDS_1D#
x <- fit1d$points[,1]
pdf(file="MDS_1D_win100_chr1.pdf")
plot(x,ylab="Coordinate 1", main="MDS_1D_chr2")
dev.off()
网友评论