親緣關(guān)系分析實(shí)操
前期準(zhǔn)備
給標(biāo)記加上ID
SNP data通常都是以VCF格式文件呈現(xiàn)冤荆，拿到VCF文件的第一件事情就是添加各個(gè)SNP位點(diǎn)的ID术瓮。
先看一下最開始生成的VCF文件：
可以看到溶弟，ID列都是"."沙廉，需要我們自己加上去扶叉。我用的是某不知名大神寫好的perl腳本嗓节，可以去我的github上下載荧缘，用法：
perl path2file/VCF_add_id.pl YourDataName.vcf YourDataName-id.vcf`
當(dāng)然也可以用excel手工添加。添加后的文件如下圖所示（格式：CHROMID__POS）：
SNP位點(diǎn)過濾（Missing rate and maf filtering）
SNP位點(diǎn)過濾前需要問自己一個(gè)問題拦宣，我的數(shù)據(jù)需要過濾嗎截粗？
一般要看后期是否做關(guān)聯(lián)分析（GWAS）；如果只是單純研究群體結(jié)構(gòu)建議不過濾鸵隧，因?yàn)檫^濾掉低頻位點(diǎn)可能會(huì)改變某些樣本之間的關(guān)系绸罗；如果需要和表型聯(lián)系其來做關(guān)聯(lián)分析，那么建議過濾豆瘫，因?yàn)樵诤笃诜治鲋械皖l位點(diǎn)是不在考慮范圍內(nèi)的珊蟀，需要保持前后一致。
如果過濾外驱，此處用到強(qiáng)大的plink軟件育灸，用法：
plink --vcf YourDataName-id.vcf --maf 0.05 --geno 0.2 --recode vcf-iid -out YourDataName-id-maf0.05 --allow-extra-chr
參數(shù)解釋：--maf 0.05：過濾掉次等位基因頻率低于0.05的位點(diǎn)；--geno 0.2：過濾掉有20%的樣品缺失的SNP位點(diǎn)昵宇；--allow-extra-chr：我的參考數(shù)據(jù)是Contig級(jí)別的磅崭，個(gè)數(shù)比常見分析所用的染色體多太多，所以需要加上此參數(shù)瓦哎。
格式轉(zhuǎn)換
將vcf文件轉(zhuǎn)換為bed格式文件砸喻。
這里注意一點(diǎn)！＝８畹骸！：應(yīng)該是軟件的問題羡铲，需要把染色體/contig名稱變成連續(xù)的數(shù)字（1 to n）蜂桶，不然會(huì)報(bào)錯(cuò)無法算出結(jié)果！（坑）
plink --vcf YourDataName-id-maf0.05.vcf --make-bed --out snp --chr-set 29 no-xy
參數(shù)解釋：--chr-set 給出染色體/contig的數(shù)目也切；no-xy 沒有xy染色體扑媚。
用gcta做親緣關(guān)系分析
gcta輸出grm陣列（genetic relationship matrix）
gcta64 --make-grm-gz --out snp.gcta --bfile snp --autosome-num 29
參數(shù)解釋：--autosome-num常染色體數(shù)目。
snp.gcta結(jié)果文件：
解讀：第一雷恃，第二列為樣品編號(hào)疆股；第三列為兩樣品間有多少個(gè)有效位點(diǎn)；第四列為兩樣品間的親緣關(guān)系的值倒槐。
將上述陣列轉(zhuǎn)化為矩陣形式
snp.gcta結(jié)果文件列舉了兩個(gè)樣品間的關(guān)系旬痹，我們需要把它變成常見的矩陣形式，這里用R可以輕松完成讨越，寫好的R包我放在了Github中：GRM2normal_format.R两残，大家自行下載使用。
如果不想下載把跨，可以復(fù)制如下代碼：
library(reshape2)
tmp <- read.table(gzfile("snp.gcta.grm.gz"), header = F, stringsAsFactors = F)
ids <- read.table("snp.gcta.grm.id", header = F, stringsAsFactors = F)
tmp <- tmp[,c(1,2,4)]
result_matrix <- acast(tmp, V1~V2, value.var="V4", drop = F)
makeSymm <- function(m) {
m[upper.tri(m)] <- t(m)[upper.tri(m)]
return(m)
}
result_full <- makeSymm(result_matrix)
result_df <- as.data.frame(result_full)
row.names(result_df) <- ids$V2
colnames(result_df) <- ids$V2
write.table(result_df, file = "ldak.weight.kinship.txt", row.names = T, col.names = NA, sep = "\t", quote = F)
需要用到上訴步驟生成的2個(gè)結(jié)果文件：snp.gcta.grm.gz和snp.gcta.grm.id人弓。
轉(zhuǎn)換后的結(jié)果文件：
解讀：第一列和第一行都是對(duì)應(yīng)的樣品名稱。
用LDAK做親緣關(guān)系分析
相比gcta着逐，能用LD對(duì)結(jié)果進(jìn)行校正崔赌，具體來說，就是先用LD計(jì)算每個(gè)SNP位點(diǎn)的權(quán)重耸别，根據(jù)權(quán)重再計(jì)算Kinship健芭，這樣的結(jié)果更接近真實(shí)情況。
LDAK輸出grm陣列（genetic relationship matrix）
- 在不考慮權(quán)重的情況下秀姐，方法如下：
ldak5.linux --calc-kins-direct snp.ldak --bfile snp --ignore-weights YES --kinship-gz YES --power -0.25
- 用LD計(jì)算每個(gè)SNP位點(diǎn)的權(quán)重慈迈，根據(jù)權(quán)重再計(jì)算Kinship
#切割
ldak5.linux --cut-weights snp.sections --bfile snp
#查看有多少個(gè)section
cat snp.sections/section.number
#根據(jù)自己的section個(gè)數(shù)分別計(jì)算權(quán)重（我這里是31個(gè)）
for section in {1..31}; do ldak5.linux --calc-weights snp.sections --bfile snp --section $section; done
#weight文件整合，給SNP賦權(quán)重值
ldak5.linux --join-weights snp.sections --bfile snp
#輸出grm陣列
ldak5.linux --calc-kins-direct snp.ldak.weight --bfile snp --weights snp.sections/weights.all --kinship-gz YES --power -0.25
結(jié)果文件和gcta類似省有，包含兩個(gè)文件：snp.ldak.weight.grm.gz和snp.ldak.weight.grm.id痒留，用上述同一個(gè)R包，同樣的方法轉(zhuǎn)化為矩陣形式锥咸。
同樣的狭瞎，如果不想下載，直接復(fù)制如下代碼：
library(reshape2)
tmp <- read.table(gzfile("snp.ldak.weight.grm.gz"), header = F, stringsAsFactors = F)
ids <- read.table("snp.ldak.weight.grm.id", header = F, stringsAsFactors = F)
tmp <- tmp[,c(1,2,4)]
result_matrix <- acast(tmp, V1~V2, value.var="V4", drop = F)
makeSymm <- function(m) {
m[upper.tri(m)] <- t(m)[upper.tri(m)]
return(m)
}
result_full <- makeSymm(result_matrix)
result_df <- as.data.frame(result_full)
row.names(result_df) <- ids$V2
colnames(result_df) <- ids$V2
write.table(result_df, file = "ldak.weight.kinship.txt", row.names = T, col.names = NA, sep = "\t", quote = F)
數(shù)據(jù)可視化
用R畫熱圖即可搏予，各種熱圖的畫法熊锭，我另外找個(gè)時(shí)間再詳細(xì)說明，先直接分享一下我作圖的代碼：
library(pheatmap)
kinship <- read.table("ldak.weight.kinship.txt", header = F, row.names = 1, skip=1)
colnames(kinship) <- row.names(kinship)
diag(kinship) <- NA
hist_data <- hist(as.matrix(kinship), xlab = "Kinship", col = "red", main = "Histogram of Kinship")
pheatmap(kinship, fontsize_row = 0.3, fontsize_col = 0.3)
color = colorRampPalette(c("white", "red","red4"),bias=0.5)(500)
#調(diào)整cell大小
par(mar=c(10,10,10,10))
pheatmap(kinship,color=color,border_color = F,fontsize_row = 0.3, fontsize_col = 0.3,cellwidth = 2,cellheight = 2)
#調(diào)整聚類樹高度
pheatmap(kinship,color=color,border_color = F,fontsize_row = 0.3, fontsize_col = 0.3,cellwidth = 2,cellheight = 2,treeheight_col= 40,treeheight_row = 40)
