用於自動處理高通量測序（RNA-seq）資料的R指令碼

標籤：

反饋方式：

1. 本文的任何錯誤，請在留言中指正；也可發郵件至[email protected]，歡迎交流；
2. 對於任何關於新功能的建議，也可按上一步交流；

本程式待改進地方：

1. 想著，在運行程式的同時，程式會將自身複製一份到輸出檔案夾用於備份（current_file_path_getter）；但是該函數的可移植性很差，暫時無法識別以R CMD方式運行該指令碼，但是通過source("")和“R --file=”方式運行時沒有問題的；也請有更好方法的牛人不吝賜教，謝謝先；
2. 期望在每一步完成後，指令碼自動發郵件到指定郵箱提醒或告知使用者；但是，“sendmailR”這個程式包，不太會用，研究中；
3. 該指令碼，寫得臭長臭長的，很多地方可以更簡練；之所以這麼長，主要是個人看著清晰；

腳步說明：

1. 本指令碼可以免費、自由使用；
2. 本指令碼目前只使用於單端測序；

使用方法：

1. 安裝必要的軟體，這個應該不是問題的；主要有Tophat，Cufflinks，FastQC等；
2. 建立一個工作目錄，並將腳步的“working_directory”設定為該建立的工作目錄；
3. 在工作目錄下，建立“input_files”夾；
4. 在“input_files”檔案夾下，建立“sequenced_reads”檔案夾。並將且僅將所有的測序資料檔案放入該檔案夾下；
5. 根據樣本的物種，在“input_files”檔案夾，建立以物種命名的檔案夾（例“Homo_sapiens”），同時將指令碼中的“specie_name”設定為物種名。
6. 在以物種命名的檔案夾（Homo_sapiens）下，建立用於存放Bowtie2 index檔案的檔案夾“Bowtie2Index”和用於存放基因注釋的“Genes”檔案夾；
7. 到IGenomes首頁下載感興趣的對應物種的參考檔案和注釋（Reference Sequences and Annotations）。並將iGenomes壓縮包中對應的檔案放入上一步建立的檔案中。
8. 根據所啟動並執行系統選擇相應的方法，運行該指令碼；
9. 運行約一兩分鐘後，本指令碼會在工作目錄下產生“output_files_××××××××_××××××”檔案夾，進入該檔案夾找到“accepted_hits_bam.txt”檔案，按以下步驟編輯（考慮到不同機器的運行速度，更適合編輯該檔案的時機是，“output_files_××××××××_××××××”檔案夾下的“tophat_output”檔案夾不為空白）：
   1. 在檔案開頭插入一新行，輸入：“sample_id　　group_label”（注，不包括引號，中間已定位字元tab分割；
   2. 在其他行後，添加一個定位字元tab，然後輸入改行所對應樣本的標籤（label，用於分組），組內不同重複請使用相同的標籤；
   3. 儲存該檔案。
10. 等待機器運行完成，就好了；
11. 使用cummeRbund進行下遊分析；

 

代碼如下：

#! /data02/software/R/bin/R## Created by Roger Young, 2015-02-05## nohup /data02/software/R/bin/R --file=~/mRNA_workflow.R &> ./nohup.log## nohup /data02/software/R/bin/R CMD BATCH ~/mRNA_workflow.R &> ./nohup.log##-------Those lines should be changed accordingly-----------------------------## 正常情況下，需要根據需要更改以下幾行：output_folder和output_folder；## directory: 代表一個目錄的絕對路徑, 以"/"結尾；另有說明的除外！## foler: 代表一個檔案夾的檔案名稱，以"/"結尾；另有說明的除外！## path: 代表一個檔案的絕對路徑；另有說明的除外！DEBUG_FLAG <- FALSE;## R路徑: /data02/software/R/bin/R## source("~/Documents/mRNA_profiles/mRNA_workflow.R")working_directory <- "~/Documents/mRNA_profiles/";## 如果在Linux系統下測試，請酌情修改。## C:\Program Files\R\R-3.1.2\bin\i386\R.exe --file=~/mRNA_workflow.Rif(DEBUG_FLAG){    working_directory <- "D:/roger/data/mRNA/";}output_folder <- paste("output_files_",                       format(Sys.time(), "%Y%m%d_%I%M%S"),"/", sep="");output_directory <- paste(working_directory, output_folder, sep ="");dir.create(path = output_directory, recursive = TRUE, showWarnings = FALSE);##-------所用到的軟體目錄，需根據具體裝置更改------------------------------------------------fastqc_path <- "~/bin/FastQC/fastqc";tophat_app_path <- "tophat";cufflinks_app_path <- "cufflinks";cuffmerge_app_path <- "cuffmerge";cuffdiff_app_path <- "cuffdiff";## Parameters used by all those softwares!threads_params <- "-p 10";debug_flag <- "    DEBUG FLAG: ";##-------輸出檔案夾，一般無需修改----------------------------------------------------log_folder <- "log/";log_directory <- paste(output_directory, log_folder,sep = "");dir.create(path = log_directory, recursive = TRUE, showWarnings = FALSE);sink(file = paste(log_directory, "log_",                  format(Sys.time(), "%Y%m%d_%I%M%S"), ".log", sep=""),      split = TRUE);##-----擷取本指令碼的路徑，並儲存在輸出檔案夾中（注，尚需修改以保證可移植性）--------------current_file_path_getter <- function() {    cmdArgs <- commandArgs(trailingOnly = FALSE);        print(paste("cmdArgs: ",paste(commandArgs(), collapse = " ")));        R_CMD_indicator <- "[[:blank:]]CMD[[:blank:]]";        CMD_MODE_matched <- grep(pattern = R_CMD_indicator,                              x = cmdArgs);        if(length(CMD_MODE_matched) > 0){            }else{                file_arg_indicator <- "--file=";                matched <- grep(pattern = file_arg_indicator, x = cmdArgs);                if (length(matched) > 0) {            # Rscript            return(normalizePath(sub(pattern = file_arg_indicator,                                      replacement = "",                                     x = cmdArgs[matched])));        } else {            # ‘source‘d via R console            return(normalizePath(sys.frames()[[1]]$ofile));        }    }}current_script_path <- current_file_path_getter();print(paste("current_script_path: ", current_script_path));file.copy(from = current_file_path_getter(),          to = paste(output_directory,"/RScript.R",sep = ""),          overwrite = TRUE, copy.mode = TRUE, copy.date = TRUE);##-------定義輸入檔案和檔案夾-----------------------------------------------input_folder <- "input_files/";input_directory <- paste(working_directory, input_folder, sep = "");sequenced_reads_folder <- "sequenced_reads"; ## 沒有"/"sequenced_reads_directory <- paste(input_directory,                                    sequenced_reads_folder, sep = "");##-------Genome References-----------------------------------------------------specie_name <- "Homo_sapiens/";annotation_file_name <- "genes.gtf";annotation_file_folder <- "Genes/";annotation_index_directory <- paste(input_directory, specie_name,                                     annotation_file_folder , sep = "");annotation_index_path <-paste(annotation_index_directory,                               annotation_file_name, sep = "");bowtie_index_name <- "genome";bowtie_index_folder <- "Bowtie2Index/";bowtie_index_directory <- paste(input_directory, specie_name,                                 bowtie_index_folder , sep = "");bowtie_index_path <-paste(bowtie_index_directory,                           bowtie_index_name, sep = "");##-------用於確認所依賴的包已正確安裝------------------------------------------check_bioconductor_pkg <- function(pkgs) {        if(require(pkgs , character.only = TRUE)){        print(paste("******", pkgs,"is loaded correctly", "******"));    }     else {        print(paste("******Trying to install ", pkgs, "******"));        ## Update R packages        print("******Update R packages******");        update.packages(ask = FALSE);                source("http://bioconductor.org/biocLite.R");        biocLite();        ## Install desired Packages        biocLite(pkgs);         ## Re-check        if(require(pkgs, character.only = TRUE)){            print(paste("******",pkgs, "installed and loaded!", "******"));        }         else {            stop(paste("??????","Error: could not install",pkgs, "??????"));        }    }} check_pkg <- function(pkgs) {        if(require(pkgs , character.only = TRUE)){        print(paste("******", pkgs,"is loaded correctly", "******"));    }     else {        print(paste("******Trying to install ", pkgs, "******"));        #         ## Update R packages        print("******Update R packages******");        update.packages(ask = FALSE);        install.packages(pkgs);        ## Re-check        if(require(pkgs, character.only = TRUE)){            print(paste("******",pkgs, "installed and loaded!", "******"));        }         else {            stop(paste("??????????","Error: could not install",                       pkgs, "??????????"));        }    }} ## (check_pkg("sendmailR"));## (check_bioconductor_pkg("cummeRbund"));##-------擷取所需處理的測序檔案清單------------------------------------------(sequenced_files <- list.files(path = sequenced_reads_directory,                               pattern = "*[^(.md5)|(.txt)]$",                               all.files = TRUE,                               full.names = TRUE,                               recursive = TRUE));##-------測序品質評估--------------------------------(fastqc_output_folder <- "fastqc_output/");(fastqc_output_directory <- paste(output_directory,                                  fastqc_output_folder,sep = ""));(dir.create(path = fastqc_output_directory, recursive = TRUE, showWarnings = FALSE));(sequenced_file_lists <- sequenced_files);print(paste("sequenced_file_lists: ", sequenced_file_lists));for (files_to_test in sequenced_file_lists){    print(paste("files_to_test: ", files_to_test));    fastqc_params <- "--outdir=";    fastqc_output_directory;    ## Tedious lines for output folder construction    temp_output_directory_dot_repaced <- gsub(pattern = "\\.",                                              replacement = "_",                                              x = files_to_test);        fastqc_output_directory <- gsub(pattern = paste(input_folder,                                                    sequenced_reads_folder,                                                     "/", sep=""), ## specical                                    replacement = paste(output_folder,                                                        fastqc_output_folder,                                                        sep = ""),,                                    x = temp_output_directory_dot_repaced);        print(paste(debug_flag, "fastqc_output_directory: ", fastqc_output_directory));        dir.create(fastqc_output_directory, recursive = TRUE);        ## 構建用於運行fastqc的命令    fastqc_cmd <- paste(fastqc_path, files_to_test,                        paste( fastqc_params, fastqc_output_directory,                               sep = ""),                        sep = " ");    print(paste(debug_flag,"fasqc_cmd: ", fastqc_cmd));    ##品質評估與後續部分無聯絡，所以在後台運行，以提高速度！    if(!DEBUG_FLAG){        print("Running fastqc!")        system(command = fastqc_cmd, wait = FALSE);     } else {        print("In Debug mode!")    }}##-------Align the RNA-seq reads to the genomes--------------------------------## Map the reads for each sample to the reference genome(tophat_output_folder <- "tophat_output/");(tophat_output_directory <- paste(output_directory,                                   tophat_output_folder, sep = ""));(dir.create(path = tophat_output_directory, recursive = TRUE, showWarnings = FALSE));(cufflinks_output_folder <- "cufflinks_output/");(cufflinks_output_directory <- paste(output_directory,                                      cufflinks_output_folder, sep = ""));(dir.create(path = cufflinks_output_directory, recursive = TRUE, showWarnings = FALSE));(sequenced_file_lists <- sequenced_files);(cuffmerge_assemblies_file <- "assemblies.txt")(cuffmerge_assemblies_file_text <- "");(cuffmerge_assemblies_file_path <- paste(output_directory,                                          cuffmerge_assemblies_file,                                         sep = ""));(accepted_hits_bam_file <- "accepted_hits_bam.txt")(accepted_hits_bam_file_text <- "");(accepted_hits_bam_file_path <- paste(output_directory,                                       accepted_hits_bam_file,                                      sep = ""));## 該部分用於構建cuffdiff所需要用的accepted_hits_bam_file_path檔案；## 該部分運行完後需要按要求修改accepted_hits_bam_file_path檔案。for(sequenced_read in sequenced_file_lists){        print(paste(debug_flag, "The file to process: ", sequenced_read));        temp_output_directory_dot_repaced <- gsub(pattern = "\\.",                                              replacement = "_",                                              x = sequenced_read);        tophat_output_directory <- gsub(pattern = paste(input_folder,                                                    sequenced_reads_folder,                                                     "/", sep=""), ## specical                                    replacement = paste(output_folder,                                                        tophat_output_folder,                                                        sep = ""),,                                    x = temp_output_directory_dot_repaced);        accepted_hits_bam_file_text <- paste(accepted_hits_bam_file_text,                                         tophat_output_directory, "/",                                         "accepted_hits.bam",                                         "\n",                                         sep = "");        cufflinks_output_directory <- gsub(pattern = paste(input_folder,                                                       sequenced_reads_folder,                                                        "/", sep=""), ## specical                                       replacement = paste(output_folder,                                                           cufflinks_output_folder,                                                           sep = ""),,                                       x = temp_output_directory_dot_repaced);        cuffmerge_assemblies_file_text <- paste(cuffmerge_assemblies_file_text,                                            cufflinks_output_directory, "/",                                            "transcripts.gtf",                                            "\n",                                            sep = "");    }write(cuffmerge_assemblies_file_text,       file = cuffmerge_assemblies_file_path);write(accepted_hits_bam_file_text,       file = accepted_hits_bam_file_path);accepted_hits_bam_file_modify_time <- file.info(accepted_hits_bam_file_path);## 給使用者發郵件，告知accepted_hits_bam_file_path檔案已建立，請修改；print(paste("Please edit file: ", accepted_hits_bam_file_path));print(paste("First, insert the following line (seperated by tab): \"sample_id group_label\";"));print("The, Append each other lines with a tab, followed by the sample‘s label (group)")print("Finally, Save the file!!")number_of_files_had_processed <- 0; ##length(sequenced_file_lists);print(paste(debug_flag, "number_of_files_had_processed: ", number_of_files_had_processed));for(sequenced_read in sequenced_file_lists){        print(paste(debug_flag, "The file to process: ", sequenced_read));        temp_output_directory_dot_repaced <- gsub(pattern = "\\.",                                              replacement = "_",                                              x = sequenced_read);        tophat_output_directory <- gsub(pattern = paste(input_folder,                                                    sequenced_reads_folder,                                                     "/", sep=""), ## specical                                    replacement = paste(output_folder,                                                        tophat_output_folder,                                                        sep = ""),,                                    x = temp_output_directory_dot_repaced);        print(paste(debug_flag, "Tophat output folder: ",tophat_output_directory));        (dir.create(path = tophat_output_directory, recursive = TRUE, showWarnings = FALSE));        ##構建運行tophat的命令    (tophat_cmd <- paste(tophat_app_path, threads_params, "-G",                         annotation_index_path,                          "-o",                         tophat_output_directory,                          bowtie_index_path,                          sequenced_read, sep = " "));    print(paste(debug_flag, "tophat_cmd: ",tophat_cmd));        if(!DEBUG_FLAG){        system(command = tophat_cmd, wait = TRUE);    }            ## 構建運行cufflinks的命令    cufflinks_output_directory <- gsub(pattern = paste(input_folder,                                                       sequenced_reads_folder,                                                        "/", sep=""), ## specical                                       replacement = paste(output_folder,                                                           cufflinks_output_folder,                                                           sep = ""),,                                       x = temp_output_directory_dot_repaced);        print(paste(debug_flag, "cufflinks output directory", cufflinks_output_directory))        (dir.create(path = cufflinks_output_directory, recursive = TRUE, showWarnings = FALSE));        (cufflinks_cmd <- paste(cufflinks_app_path, threads_params, "-o",                            cufflinks_output_directory,                             paste(tophat_output_directory,"/",                                  "accepted_hits.bam", sep = ""),                            sep = " "));    print(paste(debug_flag, "cufflinks_cmd: ",cufflinks_cmd));        ##同步運行cufflinks！    number_of_files_had_processed <- number_of_files_had_processed + 1;    print(paste(debug_flag, "Number of files had been handled: ", number_of_files_had_processed));        if(!DEBUG_FLAG){        system(command = cufflinks_cmd, wait = TRUE);    }}warnings()##-------Run Cuffmerge to create a single merged transcriptome annotation------(cuffmerge_output_folder <- "cuffmerge_output"); ##沒有“/”(cuffmerge_output_directory <- paste(output_directory,                                      cuffmerge_output_folder, sep = ""));dir.create(path = cuffmerge_output_directory, recursive = TRUE, showWarnings = FALSE);cuffmerge_cmd <- paste(cuffmerge_app_path, "--ref-gtf", annotation_index_path,                        "-o", cuffmerge_output_directory,                       "--ref-sequence", paste(bowtie_index_directory,"genome.fa",                                               sep = ""),                        threads_params,                        cuffmerge_assemblies_file_path,                       sep = " ");print(paste(debug_flag, "cuffmerge_cmd: ",cuffmerge_cmd));if(!DEBUG_FLAG){    system(command = cuffmerge_cmd, wait = TRUE);}##-------Identify differentially expressed genes and transcripts------(cuffdiff_output_folder <- "cuffdiff_output/");(cuffdiff_output_directory <- paste(output_directory,                                     cuffdiff_output_folder, sep = ""));(dir.create(path = cuffdiff_output_directory, recursive = TRUE, showWarnings = FALSE));cuffdiff_cmd <- paste(cuffdiff_app_path, "-use-sample-sheet",                      "-o", cuffdiff_output_directory,                      "-b", paste(bowtie_index_directory,"genome.fa",                                  sep = ""),                      threads_params,                      "-u", paste(cuffmerge_output_directory,"/","merged.gtf",sep=""),                      accepted_hits_bam_file_path,                      sep=" ")print(paste(debug_flag, "cuffdiff_cmd: ",cuffdiff_cmd, sep = ""));if(DEBUG_FLAG){    Sys.sleep(time = 20);}## 此處，我們認為使用者已正確修改了accepted_hits_bam_file_path所指代的檔案if(identical(file.info(accepted_hits_bam_file_path),              accepted_hits_bam_file_modify_time)){    print(paste("This file stays unchanged: ", accepted_hits_bam_file_path, sep = ""));    print("cuffdiff is skipped! Please run it manually!!");    print(paste("The cmd for running cuffdiff is: ",cuffdiff_cmd, sep = ""));} else {       print("cuffdiff_cmd is running");    if(!DEBUG_FLAG){        system(command = cuffdiff_cmd, wait = TRUE,);    }}print("Done!");sink();print("Done!");if(!DEBUG_FLAG){    quit(save = "no");}

　　

完成！

用於自動處理高通量測序（RNA-seq）資料的R指令碼