如何利用R將多個定序檔案合併起來

前言

由於定序技術(Sequencing)已經廣泛用於個人化醫療、遺傳疾病診斷、癌症基因變異的鑒定以及病原體的追蹤和研究等,再加上近年來的定序價格也越來越低,所以很多實驗室都會以次世代定序(Next Generation Sequencing)、RNA-定序等方式來探討基因組結構、變異、表達和調控機制。但是當樣本多達數十上百個時,在前端的資料整理就足以讓人頭昏眼睛裂,再加上實驗室錢錢不足無法購買分析軟體的情況下,老闆大多都只能犧牲分析者的眼睛去手動完成這些事情。既然實驗室沒錢,又不想要眼睛裂掉,R語言就是我們的最佳朋友,以下是將所有定序檔案合併起來的方法,提供大家參考。

前置作業

本篇文章會使用tidyversejanitor這兩個套件來進行資料處理和清理,教學檔案可以在這裡下載

# 下載檔案
install.packages("tidyverse")
install.packages("janitor")
# 檢查版本
packageVersion("tidyverse")
packageVersion("janitor")

下載完tidyversejanitor套件後可以使用packageVersion()函數檢查版本,本篇所使用的tidyverse套件版本為v2.0.0,janitor套件版本為v2.2.0,接著我們便直接進入操作流程囉!

多個定序檔案合併範例

以下是檔案合併步驟:

  1. 載入所需的套件:首先載入tidyversejanitor這兩個套件。
  2. 列出CSV檔案:建立了一個檔案路徑列表,包含指定目錄中所有檔名結尾為.csv的CSV檔案。
  3. 資料處理:使用purrr套件中的map_df()函數來處理列表中的每個CSV檔案,對於每個檔案,執行以下操作:
    • a. 使用readr套件中的read_csv()函數讀取CSV檔案。
    • b. 使用janitor套件中的clean_names()函數來清理欄位名稱。
    • c. 添加一個名為symbol的新欄位,包含來自gene_homo_sapiens_refseq_gr_ch37_p13_genes欄位的值。
    • d. 篩選reference_allele欄位等於No的資料列。
    • e. 使用dplyr套件中的across()函數將所有欄位轉換為字元類型。
  4. 修改file_name:將file_name欄位修改為只保留檔案的基本名稱(不包含目錄路徑和.csv副檔名)。
  5. 合併結果:使用map_df()函數將處理過的每個CSV檔案的資料框結合成一個大的資料框。file_name欄位用於追蹤每一行資料來自哪個原始檔案。
  6. 檢視結果:使用View()函數檢視合併且處理過的資料框(Merge)。
  7. 匯出資料:將處理過的Merge資料框分別以CSV檔案和RDS檔案(R Data Storage)的形式保存在指定的路徑上,注意需更換成自己電腦的路徑
    • CSV檔案保存為”C:/Users/USER/Desktop/merge.csv”。
    • RDS檔案保存為”C:/Users/USER/Desktop/merge.rds”。
# 輸入套件
library(tidyverse)
library(janitor)

# 輸入資料夾
output_list <- list.files(path="C:/Users/USER/Desktop/WES_sample", pattern="*.csv$",
                          full.names = T)
# 整合所有.CSV檔案
Merge <- output_list %>%
  setNames(nm = .) %>% 
  map_df(~read_csv(.x)%>% 
           clean_names() %>%
           mutate(symbol=gene_homo_sapiens_refseq_gr_ch37_p13_genes)%>%
           subset(reference_allele=="No") %>%
           mutate(across(.fns = as.character)), .id = "file_name") %>%
  mutate(file_name=gsub(".*?/", "",file_name)) %>%
  mutate(file_name=gsub(".csv", "",file_name))

# 檢視整合結果
View(Merge)

# 儲存整合結果
write.csv(file = "C:/Users/USER/Desktop/merge.csv", x = Merge, row.names = FALSE)
saveRDS(Merge, file = "C:/Users/USER/Desktop/merge.rds")

將上述指令執行完畢後,我們就能得到一個合併好65筆資料的結果,生成的CSV檔案包含指定目錄中所有CSV檔案的合併和處理過的CNV資料,而RDS檔案將以一種可以在後續時間重新讀取到R中的格式保存相同的資料。

要如何從檔案中抽取所需的數據

當我們將所有資料整合到一個檔案後,接下來我們就可以在Merge檔案中進行資料的提取,以下是其中一種提取方法:

  1. 定義需要提取的基因列表:將要提取的基因名稱以管道符號(|)分隔,並賦值給變數gene.list。例如,gene.list="NOC2L|PERM1|OR4F5"表示要提取NOC2L、PERM1和OR4F5這三個基因的結果。
  2. 處理結果:使用dplyr套件中的mutate()函數對Merge資料框進行資料處理。
  3. 篩選結果:使用subset()函數,基於symbol欄位是否包含在gene.list中的條件並進行篩選,只保留符合條件的資料。
  4. 檢視結果:使用View()函數以表格形式檢視處理後的結果資料框Results
# 選擇需要提取的基因
gene.list="NOC2L|PERM1|OR4F5"

# 結果
Results <- Merge %>% mutate(amino_acid_change_in_longest_transcript=gsub(".*?:", "",amino_acid_change_in_longest_transcript)) %>%
  mutate(coding_region_change_in_longest_transcript=gsub(".*?:", "",coding_region_change_in_longest_transcript)) %>%
  subset(str_detect(symbol,gene.list))

# 檢視結果
View(Results)

結語

結果為從Merge數據集中的65筆資料提取了其中的35筆出來。實話實說,使用這種方法提取出來的資料,有時候會參雜其他基因名稱相近的資料,需要自行篩選檢查。雖然有點小缺點,但這也讓我們省下了非常多的時間,如果各位先進有更好提取資料的方法,歡迎在下方留言處分享,我將不勝感激!!!

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

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

文章: 46

發佈留言

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

error: Alert: Content selection is disabled!!