FPKM转TPM脚本分享：当Shell脚本嵌入R代码

当Shell脚本中嵌入R代码，这到底是R脚本，还是shell脚本呢？

以下内容均为使用中的个人感受或看法，大家谨慎使用。

Shell脚本中嵌入R代码，听起来有点花里胡哨的。我个人的看法是，还是有一点用处的。

对于R脚本，一般需要安装一些R包，如optparse、getopt、argparse和argparser才能使用长短选项设定参数。这使得脚本的可移植性可能变差，因为有些运行环境是不可以随便更改或者不允许更改的。当然，我之前写过一个R脚本模板，可以只依赖基础函数也允许使用长选项设定参数，也可以设定缺省值，详见《【R脚本模板3】谁说arg <- commandArgs(T)脚本不能设置选项参数》。但是写完这个模板，经过一些简单测试之后，我也还没用到生产或者自己流程脚本中。

而shell编程可以使用getopts/getopt等工具来设定选项和缺省值，这些工具在常用系统环境中是默认配置，因此使得脚本更加灵活【下面脚本并未这样写，大家可以自己改写】。

本次是以FPKM表达矩阵转TPM矩阵为例封装的脚本，下面带大家看下使用方法。

脚本使用

输入文件：FPKM表达矩阵，行名为基因，列名为样本名称，值为FPKM。文件以tab制表符分割，文件名为*_FPKM.xls格式，因为是特定情境下的脚本，所以卡的比较严格。

结果文件：与FPKM表达矩阵格式相同，文件格式为*_TPM.xls

使用：脚本名称为run_all_FPKM2TPM_R.sh

$ ./run_all_FPKM2TPM_R.sh result_genes_FPKM.xls 
==> Start to convert 'result_genes_FPKM.xls' <==                                            
Trying to read result_genes_FPKM.xls...                                                     
Convertting FPKM to TPM...                                                                  
Result: 'result_genes_TPM.xls'                                                              
==> Succeed to convert 'result_genes_FPKM.xls' <==

$ ll result_genes*
-rwxrwxrwx 1 zheng zheng 8.8M Jan 10 20:31 result_genes_FPKM.xls*                           
-rwxrwxrwx 1 zheng zheng 8.7M Jan 25 02:39 result_genes_TPM.xls*

脚本代码

• 需要注意 shell中参数传到R中参数的方法
• 注意R中data.frame 中索引列的写法，不要用$ 了，防止歧义。
• 为了防止R读入数据时被修改列名，这里用readLines()读取第一行后转为列名。
• 对于TPM矩阵结果，第一列的列名与FPKM中第一列列名相同。
• 本次R代码中仅用了基础函数

#!/bin/bash

input="$1" # ./*_FPKM.xls  # 结果与其相同文件夹

Rscript - <<EOF
    options(stringsAsFactors = F)

    # input <- "./ALL_FPKM.xls"  
    input <- "$input" # 必须加引号。原因？
    if(!grepl("_FPKM.xls", input))stop("The input must be '{Sample}_FPKM.xls'")
    output <- gsub(pattern = "_FPKM.xls",replacement = "_TPM.xls", input)

    # chcek input 
    if(!file.exists(input)) stop("FileNotExist: '",input,"'")
    message("==> Start to convert '",input,"' <==")

    # read permission of input # 可读时为0
    if(file.access(input, mode = 4) != 0) stop("File '",input,"' does not have read permission")

    message("Trying to read ",input,"...")
    fpkm <- read.table(input, header = TRUE, sep = "\t",comment.char = "")
    header <- strsplit(readLines(input, n = 1L), "\t")[[1]]

    colnames(fpkm) <- header
    
    geneid <- colnames(fpkm)[1] # "GeneID"
    #if(!"GeneID" %in% colnames(fpkm[,2:ncol)) stop("Not find 'GeneID' column in '",input,"'")
    # if(!geneid %in% colnames(fpkm[,2:ncol)) stop("Not find '"geneid"' column in '",input,"'")
    if( sum(duplicated(fpkm[,geneid])) > 0 ) stop( "There are some genes is repeated.\n", paste0(unique(fpkm[,geneid][duplicated(fpkm[,geneid])]), collapse = ", ") )


    message("Convertting FPKM to TPM...")

    tpm <- apply(fpkm[,2:ncol(fpkm)], 2, function(t){
      exp(log(t)-log(sum(t))+log(1e6))
    })

    result_tpm <- data.frame(GeneID = fpkm[,geneid], # GeneID仅作为临时基因ID名称
                             tpm
                             )
    colnames(result_tpm)[1] <- geneid # 基因ID名称转为原文件中的定义

    write.table(result_tpm, file = output, sep = "\t", quote = F, na = "", col.names = TRUE, row.names = FALSE)
    message("Result: '",output,"'")
    message("==> Succeed to convert '",input,"' <==\n")
EOF

小结

shell中写R代码的一些细节，了解的还不够详细，需要多写多练。
辑器