如何進(jìn)行R語(yǔ)言用DNA序列做主成分的分析

The function fasta2genlight extracts SNPs from alignments with fasta format. 從比對(duì)好的fasta文件中提取snp數(shù)據(jù)

adegenet這個(gè)包第一使用需要先安裝，直接運(yùn)行如下命令

install.packages("adegenet")

今天的推文使用的數(shù)據(jù)集是這個(gè)包的內(nèi)置數(shù)據(jù)集，首先是獲取這個(gè)數(shù)據(jù)集的存儲(chǔ)路徑

dfpath<-system.file("files/usflu.fasta",package="adegenet")
dfpath

library(adegenet)
flu<-fasta2genlight(dfpath,chunkSize = 10,parallel = F)
flu

數(shù)據(jù)讀入以后做一些分析就比較容易了

library(ggplot2)
snpposi.plot(position(flu),genome.size = 1700,codon = F)+
  theme_bw()

 

snpposi.plot(position(flu),genome.size = 1700,codon = T)+
  theme_bw()

 

這個(gè)圖如果分面畫(huà)成山脊圖的形式可能會(huì)更好看，但是自己目前還不知道如何實(shí)現(xiàn)

snpposi.test(position(flu),genome.size = 1700)

這一步的時(shí)間可能會(huì)比較長(zhǎng)

df.pca<-glPca(flu,nf=3)  
df.pca.scores<-as.data.frame(df.pca$scores)  
df.pca.scores

自己隨便構(gòu)造一個(gè)分組信息，然后用散點(diǎn)圖加置信橢圓的方式展示結(jié)果

df.pca.scores$population<-ifelse(df.pca.scores$PC1>0,"pop1",
                                 ifelse(df.pca.scores$PC2>1,"pop2","pop3"))
library(ggplot2)
ggplot()+
  geom_point(data=df.pca.scores,
             size=2,
             aes(x=PC1,y=PC2,
                 color=population))+
  theme_bw()+
  stat_ellipse(data=df.pca.scores,
               aes(x=PC1,y=PC2,fill=population),
               geom = "polygon",alpha=0.2,lty="dashed",color="black")