VCF檔如何在R中提取資料

[peggy]

雖然試過了vcfR這個package的extract.gt函數可以把INFO中的LP取出來,可是不知道如何將chrom跟id提取出來合併再一起,希望版上大大可以幫幫忙

程式碼:

#讀取VCF
library(vcfR)

vcf_path <- "path"

vcf <-read.vcfR(vcf_path)

#提取欄位
LP <- extract.gt(vcf, element ='LP’)

#生物資訊#.vcf

[WEPA ^_^]

Hi Peggy,

如果有幫助或是已經解決問題, 請再回覆 [已解決]

參考以下之內容: 除了 vcfR 套件, 另有 VariantAnnotation 套件可以參考, 以下說明以 vcfR 為主.

https://bioconductor.org/packages/release/bioc/html/VariantAnnotation.html

# 問題1 取出 vcf 物件的 genotype data(GT), 使用 strsplit 函數分割資料. 分割後使用 as.numeric 函數將 character 資料轉換為數值(numeric).

# 問題2 合併 CHROM 與 ID, 使用 cbind 函數.

# 使用 vcfR 套件
library(vcfR)

# 使用 pinfsc50 套件內建vcf, 並設定檔案路徑
vcf_file <- system.file("extdata", "pinf_sc50.vcf.gz", package = "pinfsc50")

# 匯入 VCF 檔案
vcf <- read.vcfR( vcf_file)

# VCF 摘要
head(vcf)

########################################
# 問題1 取出值, 分割, 轉換
########################################

########################################
# 方法1: 取出GT, 使用S4的@符號
########################################

bio_gt <- vcf@gt

# 類別: 矩陣, 陣列
class(bio_gt) # "matrix" "array"

# 維度
dim(bio_gt) # 22031*19

# 列名稱
rownames(bio_gt) # NULL

# 行名稱
colnames(bio_gt) # "FORMAT" "BL2009P4_us23" "DDR7602" "IN2009T1_us22" ....

# FORMAT: "GT:AD:DP:GQ:PL" 5種格式.
# head(vcf) 結果亦有顯示 "Unique GT formats:"
levels(factor(bio_gt[, "FORMAT"]))

# 取出前3列,前4行
bio_gt[1:3, 1:4]

# 本例 bio_gt[, 2] 其中2表示第2行: BL2009P4_us23, 其值是"1|1:0,7:7:21:283,21,0"
# 使用 strsplit 函數進行字串分割, 區隔符號是:
tmp_BL2009P4_us23 <- strsplit(bio_gt[, 2], split = ":")

# 分割結果 list
class(tmp_BL2009P4_us23)

# 前6筆資料
head(tmp_BL2009P4_us23)

# 使用 rbind 將list 合併與轉換為 matrix
matrix_gt <- do.call(rbind, tmp_BL2009P4_us23)

class(matrix_gt) # "matrix" "array"

# 取出第3行DP
matrix_gt[, 3]

# 使用 as.numeric 轉換為數值
DP_value1 <- as.numeric(matrix_gt[, 3])
DP_value1
head(DP_value1)

mean(DP_value1, na.rm = TRUE) # 29.37271

########################################
# 方法2: 使用 extract.gt 函數, 此方法保留欄位名稱較佳
########################################

# 預設取出GT值
gt_value2 <- extract.gt(vcf)
gt_value2[1:3, 1:10]

# 設定取出DP值, 此方法會保留行名稱
DP_value2 <- extract.gt(vcf, element = 'DP')
DP_value2[1:3, 1:4]

# 設定轉換資料型態為 numeric
DP_value2<- extract.gt(vcf, element = 'DP', as.numeric = TRUE)
DP_value2

mean(DP_value2[, "BL2009P4_us23"], na.rm = TRUE) # 29.37271, same as DP_value1

########################################
# 問題2 合併 CHROM 與 ID
########################################

########################################
# 方法1: 使用 S4 的@符號取出結果
########################################

bio_fix <- vcf@fix

colnames(bio_fix)
bio_fix[, "CHROM"]

# 建立合併的變數
tmp <- matrix(paste0(bio_fix[, "CHROM"], "-", bio_fix[, "ID"]),  dimnames = list(NULL, "CHROMID"))

# 使用 cbind 合併二個 matrix 為一個矩陣
bio_fix_new <- cbind(bio_fix, tmp)

dim(bio_fix_new)

########################################
# 方法2: 使用 getFIX函數 取出結果
########################################

# 取出前7欄變數, 預設值不包括INFO
myfix_without_INFO <- getFIX(vcf)
dim(myfix_without_INFO) # 22031*7

# 取出所有8欄變數, 包括INFO
myfix <- getFIX(vcf, getINFO = TRUE)
dim(myfix)

# 取出所有8欄變數名稱
# "CHROM" "POS" "ID" "REF" "ALT" "QUAL" "FILTER" "INFO"
colnames(myfix)

# 取出CHROM, 參考線上說明 ?getFIX
getCHROM(vcf)

# 取出ID
getID(vcf)

# 建立合併的變數
tmp <- matrix(paste0(getCHROM(vcf), "-", getID(vcf)),  dimnames = list(NULL, "CHROMID"))

# 使用 cbind 合併二個 matrix 為一個矩陣
myfix_new <- cbind(myfix, tmp)

dim(myfix_new )
# end

pegg...@gmail.com 在 2023年3月26日 星期日晚上9:05:41 [UTC+8] 的信中寫道：