最終拿到的數(shù)據(jù)格式如下
- 第一列是探針id
- 第二列是染色體編號(hào)
- 第三類是染色體位置
- 第四列是第一個(gè)樣本的基因型
- 接下來(lái)依次是每個(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)
第一步是讀入數(shù)據(jù)
因?yàn)槭莄sv格式存儲(chǔ)怀喉，直接使用read.csv()函數(shù)讀入
df1<-read.csv("chip_snp_example.csv",header=T)
接下來(lái)是使用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)來(lái)的數(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ù)就讀入了聊闯，但是接下來(lái)我想構(gòu)建一個(gè)分類樹(shù)，這個(gè)教程里是沒(méi)有的米诉，想起來(lái)之前重復(fù)過(guò)的一個(gè)教程里有這個(gè)內(nèi)容
教程的鏈接是
https://grunwaldlab.github.io/Population_Genetics_in_R/gbs_analysis.html
這個(gè)里的示例數(shù)據(jù)用到的是vcf格式文件菱蔬，讀入R語(yǔ)言后的數(shù)據(jù)對(duì)象是genlight，我們當(dāng)前讀入的數(shù)據(jù)是genind
那么這兩個(gè)數(shù)據(jù)格式能否相互轉(zhuǎn)化呢？
經(jīng)過(guò)搜素找到了一個(gè)R語(yǔ)言包dartR,對(duì)應(yīng)的函數(shù)是gi2gl() Converts a genind object to genlight object
第一次使用進(jìn)行安裝
install.packages("dartR")
加載的時(shí)候遇到報(bào)錯(cuò),提示沒(méi)有SNPRelate這個(gè)包拴泌，再單獨(dú)安裝一下就好了
BiocManager::install("SNPRelate")
library(dartR)
mydata1<-gi2gl(mydata)
接下來(lái)就是基于距離的upgma樹(shù)
library(poppr)
tree<-aboot(mydata1,tree = "upgma",
distance = bitwise.dist,
sample = 1000,
showtree = F)
接下來(lái)是使用ggtree進(jìn)行展示結(jié)果
library(ggtree)
ggtree(tree,layout = "circular")+
geom_tiplab()+
xlim(NA,0.12)
接下來(lái)是主成分分析PCA
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)造的沒(méi)有分組信息魏身，就畫一個(gè)簡(jiǎn)單的散點(diǎn)圖就好了
如果需要本文的示例數(shù)據(jù)的話點(diǎn)贊，然后點(diǎn)擊在看蚪腐，最后留言就好了
歡迎大家關(guān)注我的公眾號(hào)
小明的數(shù)據(jù)分析筆記本
