構建親緣樹的R語言套件-fastreeR 1.4.0

簡介

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

1. 安裝fastreeR 1.4.0

需要預先安裝R軟體以及RStudio,接著再安裝fastreeR (version 1.4.0)方法如下:

if (!requireNamespace("BiocManager", quietly=TRUE))
install.packages("BiocManager")
BiocManager::install("fastreeR")

2. 額外需安裝的套件

由於fastreeR (version 1.4.0)套件是需要在JAVA的環境下執行,所以還需另外安裝搭配JAVA套件。需求套件如下:

install.packages("rJava")
install.packages("xlsxjars")
install.packages("xlsx")

此外,還需下載JAVA,並將JAVA設定到R的環境中,JAVA_HOME路徑請查詢各自電腦相應JAVA儲存位置

#for 64-bit version #rJAVA must run in JAVA environment
Sys.setenv(JAVA_HOME='C:/Program Files/Java/jre-1.8')

3. 載入所有要使用的R套件

要使用任何的 R 套件之前,需要先使用 library ()將套件載入

library(fastreeR)
library(rJava)

4. 進行分析前的前置作業

使用fastreeR分析不同樣本的資料時,須將所有檔案合併成一個檔案才能進行分析,這裡以.vcf檔案當例子。假設今天有三套全外顯子定序 (whole exome sequencing, WES)的結果,我們必須要用vcftools將它們進行合併,程式碼如下:

vcftools --gzvcf sample1.vcf.gz --chr chr16 --recode --out sample1.vcf
vcftools --gzvcf sample2.vcf.gz --chr chr16 --recode --out sample2.vcf
vcftools --gzvcf sample3.vcf.gz --chr chr16 --recode --out sample2.vcf

--recode 表示輸出篩選的內容 –chr表示需要抽取的染色體名字 –out表示產出結果

接著進一步將文件壓縮:

bgzip sample1.vcf.recode.vcf
bgzip sample2.vcf.recode.vcf
bgzip sample3.vcf.recode.vcf

再用tabix指令對文件建立索引:

tabix -p vcf sample1.vcf.recode.vcf.gz
tabix -p vcf sample2.vcf.recode.vcf.gz
tabix -p vcf sample3.vcf.recode.vcf.gz

最後將所有檔案合併成一個文件:

vcf-merge sample1.vcf.recode.vcf.gz sample2.vcf.recode.vcf.gz sample3.vcf.recode.vcf.gz  > out.vcf

到這裡樣本的前置作業就完成囉!!!

5. 用FastreeR進行樣本的分析

fastreeR的主要功能有四

  1. Sample Statistics: 用於計算VCF文件中的樣本數據統計信息,如樣本數、SNP數、缺失值率等。
  2. Calculate distances from vcf: 根據VCF文件中的基因型數據,計算樣本之間的距離矩陣,以便後續構建親緣樹或進行層次聚類。
  3. Histogram of distances: 將距離矩陣轉化為直方圖,用於直觀展示樣本間的距離分布情況。
  4. Plot tree from fastreeR::dist2tree / Plot tree from stats::hclust: 使用fastreeR套件自帶的dist2tree函數或R基礎包stats中的hclust函數,將距離矩陣轉化為親緣樹並進行可視化。
# 載入套件
library(fastreeR)
library(rJava)
# 載入數據
tempVcf <- (inputFile = "C:/Users/USER/Desktop/out.vcf")
myVcfIstats <- fastreeR::vcf2istats(inputFile = tempVcf)
head(myVcfIstats)
plot(myVcfIstats[,7:9])
# Histogram of distances
myVcfDist <- fastreeR::vcf2dist(inputFile = tempVcf, threads = 2)
graphics::hist(myVcfDist, breaks = 100, main=NULL,
xlab = "Distance", xlim = c(0,max(myVcfDist)))
# Plot tree
myVcfTree <- fastreeR::dist2tree(inputDist = myVcfDist)
plot(ape::read.tree(text = myVcfTree), direction = "down", cex = 0.3)
ape::add.scale.bar()
ape::axisPhylo(side = 2)
# stats::hclust
myVcfTreeStats <- stats::hclust(myVcfDist)
plot(myVcfTreeStats, ann = FALSE, cex = 0.3)

運行完上述程式碼後,就可以得到類似於下方的樹狀圖囉:

圖示為使用stats::hclust指令畫出來的樹狀圖,擷取至fastreeR教學流程內文。

結語

本篇所編寫的R程式語言於fastreeR (version 1.4.0)套件中已經確認完全可行,如果大家覺得寫得還可以,就請轉發給同是生資苦命人的小夥伴們吧!

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

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

文章: 46

發佈留言

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

error: Alert: Content selection is disabled!!