R语言处理基因芯片测序得到的SNP数据的分析是怎样的

 

• 第一列是探针id
• 第二列是染色体编号
• 第三类是染色体位置
• 第四列是第一个样本的基因型
• 接下来依次是每个样本的基因型

相当于是每行是一个位点，每列是一个样本

这个数据集下载自 链接 https://popgen.nescent.org/StartSNP.html

 

今天推文的主要内筒也是参考自这个链接，这个链接的数据集是500多个样本，3000多个位点，我只选取了其中30个样本，996个位点

因为是csv格式存储，直接使用read.csv()函数读入

df1<-read.csv("chip_snp_example.csv",header=T)

接下来是使用adegenet包中的df2genind()函数对数据进行整合

要求的是

• 位点
• 倍性
• 样本名称
• 种群
• 后面这个分隔符起到什么作用我暂时还不知道(好像突然知道了，sep应该指的是AT之间的分隔符把，有的数据可能是A/T这种，那么sep就需要指定斜线了)

这个要求样本是行，位点是列，所以要对读进来的数据进行转置

df2<-t(df1)
dim(df2)
df3<-df2[4:33,1:996]
df3[1:6,1:6]

 

mydata<-df2genind(df3,ploidy = 2,ind.names = rownames(df3),
                  sep="")

到这里数据就读入了，但是接下来我想构建一个分类树，这个教程里是没有的，想起来之前重复过的一个教程里有这个内容

教程的链接是 https://grunwaldlab.github.io/Population_Genetics_in_R/gbs_analysis.html

这个里的示例数据用到的是vcf格式文件，读入R语言后的数据对象是genlight，我们当前读入的数据是genind那么这两个数据格式能否相互转化呢？经过搜素找到了一个R语言包dartR,对应的函数是gi2gl() Converts a genind object to genlight object

第一次使用进行安装

install.packages("dartR")

加载的时候遇到报错,提示没有SNPRelate这个包，再单独安装一下就好了

BiocManager::install("SNPRelate")

library(dartR)
mydata1<-gi2gl(mydata)

library(poppr)
tree<-aboot(mydata1,tree = "upgma", 
      distance = bitwise.dist, 
      sample = 1000, 
      showtree = F)

library(ggtree)
ggtree(tree,layout = "circular")+
  geom_tiplab()+
  xlim(NA,0.12)

 

df.pca<-glPca(mydata1,nf=3)  
df.pca.scores<-as.data.frame(df.pca$scores)  
df.pca.scores 
library(ggplot2)
ggplot(df.pca.scores,aes(x=PC1,y=PC2))+
  geom_point(size=2,color="blue")+
  theme_bw()

 

因为数据是随便构造的没有分组信息，就画一个简单的散点图就好了