使用R语言用DNA序列做主成分分析（PCA）的简单小例子

之前也有人在我的公众号 小明的数据分析笔记本 留言问过如何用DNA序列做主成分分析，当时我也不知道，但是大体有一个思路 就是先比对，然后把比对的数据转换成通常用的snp数据应该就可以了，但是也仅限于思路，完全不知道如何操作，今天坐车回家，路上无聊，翻了一下电脑上保存的一些资料，发现了一个办法：可以借助R语言的adegenet包，用到的函数是fasta2genlight()

fasta2genlight()函数的只要作用

The function fasta2genlight extracts SNPs from alignments with fasta format. 从比对好的fasta文件中提取snp数据

下面开始实际操作

adegenet这个包第一使用需要先安装，直接运行如下命令

install.packages("adegenet")

今天的推文使用的数据集是这个包的内置数据集，首先是获取这个数据集的存储路径

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

加载包读入数据

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

数据读入以后做一些分析就比较容易了

首先是看一下snp位点在染色体上的分布密度

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

image.png

还可以划分不同的密码子位置

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

image.png

这个图如果分面画成山脊图的形式可能会更好看，但是自己目前还不知道如何实现

还能够检测snp在染色体上是否分布均匀

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

这一步的时间可能会比较长

image.png

接下来是做主成分分析了

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

自己随便构造一个分组信息，然后用散点图加置信椭圆的方式展示结果

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")

image.png

欢迎大家关注我的公众号
小明的数据分析笔记本

今天参考的资料是一份pdf文档，里面还有好多其他的内容，大家如果需要的话可以直接留言呀！