R語言處理基因芯片測(cè)序得到的SNP數(shù)據(jù)的分析是怎樣的

 

• 第一列是探針id
• 第二列是染色體編號(hào)
• 第三類是染色體位置
• 第四列是第一個(gè)樣本的基因型
• 接下來依次是每個(gè)樣本的基因型

相當(dāng)于是每行是一個(gè)位點(diǎn)，每列是一個(gè)樣本

這個(gè)數(shù)據(jù)集下載自 鏈接 https://popgen.nescent.org/StartSNP.html

 

今天推文的主要內(nèi)筒也是參考自這個(gè)鏈接，這個(gè)鏈接的數(shù)據(jù)集是500多個(gè)樣本，3000多個(gè)位點(diǎn)，我只選取了其中30個(gè)樣本，996個(gè)位點(diǎn)

因?yàn)槭莄sv格式存儲(chǔ)，直接使用read.csv()函數(shù)讀入

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

接下來是使用adegenet包中的df2genind()函數(shù)對(duì)數(shù)據(jù)進(jìn)行整合

要求的是

• 位點(diǎn)
• 倍性
• 樣本名稱
• 種群
• 后面這個(gè)分隔符起到什么作用我暫時(shí)還不知道(好像突然知道了，sep應(yīng)該指的是AT之間的分隔符把，有的數(shù)據(jù)可能是A/T這種，那么sep就需要指定斜線了)

這個(gè)要求樣本是行，位點(diǎn)是列，所以要對(duì)讀進(jìn)來的數(shù)據(jù)進(jìn)行轉(zhuǎn)置

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

到這里數(shù)據(jù)就讀入了，但是接下來我想構(gòu)建一個(gè)分類樹，這個(gè)教程里是沒有的，想起來之前重復(fù)過的一個(gè)教程里有這個(gè)內(nèi)容

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

這個(gè)里的示例數(shù)據(jù)用到的是vcf格式文件，讀入R語言后的數(shù)據(jù)對(duì)象是genlight，我們當(dāng)前讀入的數(shù)據(jù)是genind那么這兩個(gè)數(shù)據(jù)格式能否相互轉(zhuǎn)化呢？經(jīng)過搜素找到了一個(gè)R語言包dartR,對(duì)應(yīng)的函數(shù)是gi2gl() Converts a genind object to genlight object

第一次使用進(jìn)行安裝

install.packages("dartR")

加載的時(shí)候遇到報(bào)錯(cuò),提示沒有SNPRelate這個(gè)包，再單獨(dú)安裝一下就好了

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

 

因?yàn)閿?shù)據(jù)是隨便構(gòu)造的沒有分組信息，就畫一個(gè)簡單的散點(diǎn)圖就好了