如何通过commandArgs将twoBitToFa输出从Bash传递给R？

当前Bash脚本中，$P（BED文件路径）、$G（基因组版本）、$O（输出格式）三个参数已正常定义。当$O为fasta时，通过UCSC的twoBitToFa工具能成功生成.fa文件，代码如下：

if [ "$O" = "fasta" ]
    then
        awk '{print $0" "$1":"$2"-"$3}' "$P" | 
        twoBitToFa -bed=stdin -udcDir=. "$twobit" stdout > "${P%.bed}".fa
    fi

注：此处stdout与>确实冗余，直接移除stdout即可，twoBitToFa默认输出到标准输出，重定向依然生效。

当$O为bed时，需要将$P、$G、$O以及twoBitToFa的输出内容传递给R脚本GetSeq_R.r，跳过生成.fa文件的步骤，但两次尝试均失败：

失败尝试1：将输出存入变量传递 #

尝试把twoBitToFa的输出存入myvar再传给R，但R脚本仅能收到$P、$G、$O三个参数，无法获取myvar内容：

if [ "$O" = "bed" ]
    then
        awk '{print $0" "$1":"$2"-"$3}' "$P" | 
        myvar=$(twoBitToFa -bed=stdin -udcDir=. "$twobit" stdout) \
        Rscript \
        GetSeq_R.r \
        $P \
        $G \
        $O \
        $myvar
    fi

R脚本开头代码：

args = commandArgs(trailingOnly=TRUE); print(args)

输出结果：

[1] "TAD-1_template.bed" "hg38"               "bed"  

失败尝试2：使用tee传递 #

使用tee后，twoBitToFa的输出能在控制台显示，但依然无法传递给R：

if [ "$O" = "bed" ]
    then
        awk '{print $0" "$1":"$2"-"$3}' "$P" | 
        twoBitToFa -bed=stdin -udcDir=. "$twobit" stdout | tee \
        Rscript \
        GetSeq_R.r \
        $P \
        $G \
        $O
    fi

解决方案 #

方法1：通过标准输入（stdin）传递序列数据 #

twoBitToFa的输出是多行FASTA格式内容，不适合作为命令行参数传递（会被空格/换行分割为多个参数，且存在长度限制）。正确做法是将序列数据通过标准输入传给R脚本，同时把$P、$G、$O作为命令行参数传递：

修改Bash脚本中的bed分支代码：

if [ "$O" = "bed" ]
    then
        awk '{print $0" "$1":"$2"-"$3}' "$P" | 
        twoBitToFa -bed=stdin -udcDir=. "$twobit" stdout |
        Rscript GetSeq_R.r "$P" "$G" "$O"
    fi

修改R脚本GetSeq_R.r，读取标准输入内容：

args = commandArgs(trailingOnly=TRUE)
# 获取命令行参数
bed_path <- args[1]
genome <- args[2]
output_format <- args[3]

# 读取标准输入的FASTA内容
fasta_content <- readLines(file("stdin"))
print("Received FASTA content:")
print(fasta_content)

方法2：使用临时文件传递（适合超大序列数据） #

若twoBitToFa输出内容极大，标准输入效率可能不足，可使用临时文件中转：

if [ "$O" = "bed" ]
    then
        # 创建临时文件
        temp_fasta=$(mktemp)
        awk '{print $0" "$1":"$2"-"$3}' "$P" | 
        twoBitToFa -bed=stdin -udcDir=. "$twobit" stdout > "$temp_fasta"
        # 传递临时文件路径给R
        Rscript GetSeq_R.r "$P" "$G" "$O" "$temp_fasta"
        # 删除临时文件
        rm "$temp_fasta"
    fi

对应R脚本修改：

args = commandArgs(trailingOnly=TRUE)
bed_path <- args[1]
genome <- args[2]
output_format <- args[3]
temp_fasta_path <- args[4]

fasta_content <- readLines(temp_fasta_path)
print(fasta_content)

完整修改后的Bash脚本 #

for arg in "$@"; do
  shift
  case "$arg" in
      "--path")         set -- "$@" "-P" ;; 
      "--genome")       set -- "$@" "-G" ;; 
      "--output")       set -- "$@" "-O" ;; 
      "--help")         set -- "$@" "-h" ;; 
       *)               set -- "$@" "$arg"
  esac
done

while getopts ":P:G:O:h" OPT
do
    case $OPT in
  P) P=$OPTARG;;
  G) G=$OPTARG;;
  O) O=$OPTARG;;
  h) help ;; 
  \?)
      echo "Invalid option: -$OPTARG" >&2
      usage
      exit 1
      ;; 
  :)
      echo "Option -$OPTARG requires an argument." >&2
      usage
      exit 1
      ;; 
    esac
done


num_col=$(awk "{print NF; exit}" "$P")  # 优化：无需cat，直接awk读取文件

if [ "$num_col" = 3 ]
then
echo -e "\n\n3 column bed file detected; no directional considerations for sequences \n\n"
    if [ "$G" = "hg38" ]
    then
        twobit="https://hgdownload.cse.ucsc.edu/goldenpath/hg38/bigZips/hg38.2bit"
    fi    
    if [ "$G" = "hg19" ]
    then
        twobit="https://hgdownload.cse.ucsc.edu/goldenpath/hg19/bigZips/hg19.2bit"
    fi
    if [ "$O" = "fasta" ]
    then
        # 去掉冗余的stdout
        awk '{print $0" "$1":"$2"-"$3}' "$P" | 
        twoBitToFa -bed=stdin -udcDir=. "$twobit" > "${P%.bed}".fa
    fi
    if [ "$O" = "bed" ]
    then
        awk '{print $0" "$1":"$2"-"$3}' "$P" | 
        twoBitToFa -bed=stdin -udcDir=. "$twobit" stdout |
        Rscript GetSeq_R.r "$P" "$G" "$O"
    fi
fi

内容的提问来源于stack exchange，提问作者Doda