構建親緣樹的R語言套件-SNPRelate 1.34.1

簡介

SNPRelate 1.34.1是一個用於從VCF文件中直接計算樣本之間的距離矩陣,構建親緣樹或執行層次聚類的R語言套件。其主要目地是提供方便快捷的函數,幫助使用者從序列數據中快速生成親緣樹或聚類分析結果。

安裝SNPRelate

需要預先安裝R軟體以及RStudio接著再安裝SNPRelate 1.34.1套件,方法如下:

# 安裝套件
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("SNPRelate")

進行分析前的前置作業

使用SNPRelate 1.34.1分析不同樣本的資料時,須將所有VCF結果合併成一個VCF檔案才能進行分析,我們可以用vcftools將所有定序結果進行合併,詳細合併步驟請參考構建親緣樹的R語言套件-fastreeR 1.4.0中的第四小節。

進行VCF合併樣本分析

要使用任何 R 套件之前,需要先使用library()載入需要的套件。這裡會以作者提供的VCF文件為例子,文件可以從這裡下載

1. 載入所有要使用的R套件以及數據集

# 加載套件
library(gdsfmt)
library(SNPRelate)
library(ggplot2)

# 檢查版本
packageVersion("SNPRelate")

# 載入數據集
vcf.fn <- "C:/Users/USER/Desktop/VCFExampleData.vcf"

2. 解析數據集

將VCF文件轉換成一種數據密度較低的形式(GDS),以便進行更快速的計算。如果是加載了整個基因組的結果,這個指令會運行得越慢。

# 解析數據集
snpgdsVCF2GDS(vcf.fn,"data.gds",method ="biallelic.only", verbose=TRUE)

3. 創建IBS矩陣

創建一個Identity by State(IBS)矩陣。IBS是指狀態同源,如果一段DNA在兩個或多個個體中均有一致的核甘酸序列,則可以定義該DNA片段屬性為狀態同源。

# 創建IBS矩陣
genofile <- snpgdsOpen("data.gds")
set.seed(100)
ibs.hc<-snpgdsHCluster(snpgdsIBS(genofile,num.thread=2, autosome.only=FALSE))

4. 繪製樹狀圖

# 繪製樹狀圖
rv <- snpgdsCutTree(ibs.hc)
plot(rv$dendrogram,main="Dendrogram based on IBS")

5. 差異性矩陣分析(Dissimilarity matrices)

這個指令將創建一個包含所有樣本間差異性的矩陣。如果是做X染色體相關分析,請將autosome.only= TRUE代碼更改為autosome.only= FALSE

# 差異性矩陣分析

dissMatrix  =  snpgdsDiss(genofile , sample.id=NULL, autosome.only=TRUE, remove.monosnp=TRUE,  maf=NaN, missing.rate=NaN, num.thread=2, verbose=TRUE)

6. 聚類分析

# 聚類分析
snpHCluster =  snpgdsHCluster(dissMatrix, sample.id=NULL, need.mat=TRUE, hang=0.01)
cutTree = snpgdsCutTree(snpHCluster, z.threshold=15, outlier.n=5, n.perm = 5000, samp.group=NULL, 
                        col.outlier="red", col.list=NULL, pch.outlier=4, pch.list=NULL,label.H=FALSE, label.Z=TRUE, 
                        verbose=TRUE)
cutTree
snpgdsClose(genofile)

7. 繪製差異性樹狀圖

# 繪製差異性樹狀圖
snpgdsDrawTree(cutTree, main = "Dendrogram based on dissimilarity",edgePar=list(col=rgb(0.5,0.5,0.5,0.75),t.col="black"),
               y.label.kinship=T,leaflab="perpendicular",yaxis.kinship=F)

參考文獻

Benbow Lab – Phylogeny form VCF file (https://benbowlab.github.io/Phylogeny.html

您的分享,我非常感激!!!
MillionQuesn
MillionQuesn

一個旅居台灣的異鄉人,分享突然靈光一現的亮點。

文章: 46

發佈留言

發佈留言必須填寫的電子郵件地址不會公開。 必填欄位標示為 *

error: Alert: Content selection is disabled!!