論文

The population structure and recent colonization history of Oregon threespine stickleback determined using restriction-site associated DNA-sequencing

本地pdf檔案

nihms465650。pdf

這個論文對應的資料是可以公開下載的

找到了一本電子書

https：//

bookdown。org/hhwagner1/

LandGenCourse_book/

裡面用到這篇文章的資料做了群體PCA，今天的推文我們試著重複一下這本電子書中的程式碼

如果要用這個資料的話首先得安裝R包

devtools：：install_github（“hhwagner1/LandGenCourse”）

devtools：：install_github（“hhwagner1/LandGenCourseData”）

讀取資料集

require（“LandGenCourse”）

example_df<-system。file（“extdata”， “

stickleback_data。txt

”，

package = “LandGenCourseData”）

data <- read。table（example_df，

sep=“\t”，

as。is=T，

check。names=F）

dim（data）

data［1：5，1：10］

看到了書裡介紹

library（）

和

require（）

函式的區別：

library（）

載入的包即使之前已經載入過了還是會載入一遍

require（）

如果之前載入過就不會再載入了

資料集應該是行是樣本，列是位點，總共571個樣本，10000個位點

生成每個樣本屬於哪個群體

pops <- unique（unlist（lapply（rownames（data），

function（x）{y<-c（）；y<-c（y，unlist（strsplit（x，“_”）［［1］］［1］））}）））

pops

sample_sites <- rep（NA，nrow（data））

for （i in 1：nrow（data））{

sample_sites［i］ <- strsplit（rownames（data），“_”）［［i］］［1］}

N <- unlist（lapply（pops，function（x）{length（which（sample_sites==x））}））

names（N） <- pops

N

主成分分析

pcaS<-prcomp（data，center = T）

獲取每個主成分所解釋的變異佔比

perc <- round（100*（pcaS$sdev^2 / sum（pcaS$sdev^2））［1：10］，2）

names（perc） <- apply（array（seq（1，10，1））， 1， function（x）{paste0（“PC”， x）}）

perc

用ggplot2來做個圖

首先是生成9個顏色

colors <- RColorBrewer：：brewer。pal（9， “Paired”）

作圖

library（ggplot2）

ggplot（data=pca。df，aes（x=PC1，y=PC2））+

geom_point（aes（color=pop））+

scale_color_manual（values = colors）+

theme_bw（）+

labs（x=“PC1 （37。47%）”，

y=“PC2 （3。78%）”）

用主成分3,4再做一個圖

pca。df34

<-data。frame（pcaS$x［，3：4］，

pop=sample_sites）

head（pca。df34）

ggplot（data=pca。df34，aes（x=PC3，y=PC4））+

geom_point（aes（color=pop），size=5）+

scale_color_manual（values = colors）+

theme_bw（）+

labs（x=“PC1 （2。55%）”，

y=“PC2 （0。88%）”）+

theme（legend。position = “top”，

legend。title =

element_text

（hjust=0。5））+

guides（color=guide_legend（nrow=1，

title。position = “top”））

最後是拼圖

library（patchwork）

p1+p2+plot_layout（guides = “collect”）&theme（legend。position = “top”）

推文的示例資料和程式碼大家可以直接找到開頭提到的線上電子書，或者關注公眾號直接給推文打賞

1元

，如果打賞了沒有收到我的回覆，可以加我的微信

mingyan24

催我

歡迎大家關注我的公眾號

小明的資料分析筆記本

小明的資料分析筆記本 公眾號 主要分享：1、R語言和python做資料分析和資料視覺化的簡單小例子；2、園藝植物相關轉錄組學、

基因組學

、群體遺傳學文獻閱讀筆記；3、

生物資訊學

入門學習資料及自己的學習筆記！

後續還要仔細看看這篇論文，看看其中對結果的解釋