好久不見，我們的直播又開始啦！今天，我們主要講的是人群分布，先用簡單的PCA來分析一下千人基因組的人群分布吧！
PCA分析，就是主成分分析，我博客有講過(點擊最底部的閱讀原文或複製連結http://www.bio-info-trainee.com/1232.html進行查看)。 PCA的原本目的是因為變量太多，想把它們合併成兩三個變量，從而簡化分析步驟。變量的多少代表維度的多少，一千維的數據已經無法想像了，但是二維和三維還是比較符合認知的。假設用PCA給千人基因組所有個體一個二維坐標，畫在圖上，就可以清清楚楚看到他們的分布了，是不是有規律的聚集成一個個種群。
主成分分析可以得到p個主成分，但是，由於各個主成分的方差是遞減的，包含的信息量也是遞減的，所以實際分析時，一般不是選取p個主成分，而是根據各個主成分累計貢獻率的大小選取前k個主成分。這裡貢獻率就是指某個主成分的方差佔全部方差的比重，實際也就是某個特徵值佔全部特徵值總和的比重。貢獻率越大，說明該主成分所包含的原始變量的信息越強。主成分個數k的選取，主要根據主成分的累積貢獻率來決定，即一般要求累計貢獻率達到85%以上，這樣才能保證綜合變量能包括原始變量的絕大多數信息。
先下載千人基因組計劃的數據吧(http://www.1000genomes.org/category/population/)， 涉及到2504個人，5大人種，26個亞人種！信息都存儲在integrated_call_samples_v3.20130502.ALL.panel文件裡面。
下載千人基因組計劃所有人的基因型的方法我以前講過
mkdir -p ~/annotation/variation/human/1000genomes
cd ~/annotation/variation/human/1000genomes
##ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/
nohup wget -c -r -nd -np -k -L -p ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502 &
我們隨便打開一個染色體看看裡面的數據是什麼~
zcat ALL.chr1.phase3_shapeit2_mvncall_integrated_v3plus_nounphased.rsID.genotypes.GRCh38_dbSNP_no_SVs.vcf.gz |cut -f 1-3
一般全基因組數據都是成千上萬的位點，我沒有看到教程告訴我如何挑選位點，比如http://online.cambridgecoding.com/notebooks/cca_admin/genetic-ancestry-analysis-python 我只能憑感覺挑選 allele frequency 在0.5附近的~
zcat ALL.chr1.phase3_shapeit2_mvncall_integrated_v3plus_nounphased.rsID.genotypes.GRCh38_dbSNP_no_SVs.vcf.gz |perl -alne '{/AF=(.*?);/;print join("\t",$F[2],@F[9..$#F]) if $1<0.51 && $1>0.49}' >chr1.alf.test
即使這樣，還是有9520個位點，這樣太多了，畫圖會比較卡，我就挑選了前面的一千個位點來做的。
dat[1:4,1:4]## 做好這個數據即可
# apply PCA - scale. = TRUE is highly
# advisable, but default is FALSE.
dat.pca <- prcomp(dat,
center = TRUE,
scale. = TRUE)
# print method
print(dat.pca)
# plot method
plot(dat.pca, type = "l")
summary(dat.pca)
biplot (dat.pca , scale =0,var.axes =F)
group_info <-read.table('integrated_call_samples_v3.20130502.ALL.panel',header =T,stringsAsFactors = F)
head(group_info)
pop = group_info[match(colname,group_info$sample ),'pop']
super_pop = group_info[match(colname,group_info$sample),'super_pop']
#library(devtools)
#install_github("ggbiplot", "vqv")
dat_tmp= as.data.frame(dat)
dat_tmp$sample = rownames(dat_tmp)
dat_df <- merge(dat_tmp,group_info,by='sample')
library(ggbiplot)
g <- ggbiplot(dat.pca, obs.scale = 1 ,var.scale = 1, var.axes=F,
groups = super_pop, ellipse = TRUE,
circle = TRUE)
g <- g + scale_color_discrete(name = '')
g <- g + theme(legend.direction = 'horizontal',
legend.position = 'top')
print(g)
library(ggfortify)
autoplot(prcomp(dat), data = dat_df, colour = 'super_pop')
還是老規矩，通過Google我找到了可視化方法！
用谷歌搜索來使用ggplot2做可視化（下）
就是上面代碼中的ggbiplot和ggfortify包，很容易就把千人基因組按照5個種群給分開了，當然，如果按照26個亞種會很難看，我就不秀圖片了！而且其實前兩個主成分的貢獻度都很低，說明它們兩個的把人群分開的作用力不夠。
首先是ggbiplot的圖片！
然後是ggfortify 圖片：
我只是隨機挑選的千人基因組計劃的1號染色體的1000個位點！我們看到，我們的數據區分的不是很明顯，我挑選的1000個位點沒辦法把人群清晰分開(前兩個主成分作用力太小了)，剛開始我選擇的是26個人種，更加麻煩，現在就標記5個超級人種，勉強還能看到規律。非洲人跟其他人區分挺好的。
AFR, African
AMR, Ad Mixed American
EAS, East Asian
EUR, European
SAS, South Asian
對於這個特定的分析，PCA一個最強大的部分是，一開始，我們無需指定尋找的集群數。
你們覺得哪一個好看呢？（投票ing）
參考文獻：
https://www.r-bloggers.com/computing-and-visualizing-pca-in-r/
https://cran.r-project.org/web/packages/ggfortify/vignettes/plot_pca.html
https://www.analyticsvidhya.com/blog/2016/03/practical-guide-principal-component-analysis-python/
http://stats.stackexchange.com/questions/72839/how-to-use-r-prcomp-results-for-prediction
文：Jimmy
圖文編輯：吃瓜群眾