文章來(lái)源：Local PCA Shows How the Effect of Population Structure Differs Along the Genome
群體結(jié)構(gòu)是衡量大型基因組數(shù)據(jù)集中個(gè)體之間相關(guān)性度量的系統(tǒng)化模式姥敛，通常使用降維技術(shù)，如主成分分析（PCA）來(lái)發(fā)現(xiàn)和可視化這些模式雏逾。平均關(guān)聯(lián)性是特定于位點(diǎn)的系譜樹(shù)之間關(guān)系的平均值座每，這種關(guān)聯(lián)性可能會(huì)受到相關(guān)選擇和其他因素的強(qiáng)烈影響慧库，如結(jié)構(gòu)變異等。通過(guò)使用區(qū)域PCA(localPCA)在關(guān)聯(lián)模式中描述這種小規(guī)模的異質(zhì)性，用于發(fā)現(xiàn)影響僅在幾個(gè)兆堿基的種群結(jié)構(gòu)變化劲适。
1. 安裝
install.packages("devtools")
install.packages("data.table")
devtools::install_github("petrelharp/local_pca/lostruct")
library(lostruct)
2. 所需軟件
（1）讀取VCF或BCF文件：bcftools
（2）實(shí)現(xiàn)測(cè)試分析：下載模版
3. 運(yùn)行
（1）使用 eigen_windows() 計(jì)算局部PCA坐標(biāo)
函數(shù) eigen_windows() 將數(shù)據(jù)存儲(chǔ)在矩陣中，每個(gè)SNP一行厢蒜，每個(gè)樣本一列（因此x [ i霞势，j ]是樣本 j 在位置 i 處的等位基因數(shù)）。
#兩種方式從標(biāo)準(zhǔn)格式(tped和vcf)中讀取數(shù)據(jù)
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)
#計(jì)算localPCA坐標(biāo)
pcs <- eigen_windows(snps,k=2)
#計(jì)算距離矩陣
pcdist <- pc_dist(pcs,npc=2)
這將提供一個(gè)矩陣pcs斑鸦，行是每個(gè)窗口的前兩個(gè)特征值和特征向量愕贡，以及pcdist，即這些窗口之間的成對(duì)距離矩陣巷屿。
（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()
