我正在尝试创建一个脚本，从 sam 文件的标题中删除读取组。从命令行运行的代码如下。
samtools view -H e2_20.indel.recal.dedup.bam | awk ' BEGIN {FS = "\t"} {split($2,a,":")} {if ($1 != "@RG" || ($1 =="@RG" && a[2] == "e2_20")) print}' | samtools reheader - e2_20.indel.recal.dedup.bam | samtools view -H -
这是输入文件的示例
@HD VN:1.4 GO:none SO:坐标
@SQ SN:chr10 LN:129993255
@SQ SN:chr11 LN:121843856
@RG ID:e2_20 PL:illumina
@RG ID:e2_9 PL:illumina
@PG ID:GATK IndelRealigner
运行上述命令后，输出为
@HD VN:1.4 GO:none SO:坐标
@SQ SN:chr10 LN:129993255
@SQ SN:chr11 LN:121843856
@RG ID:e2_20 PL:illumina
@PG ID:GATK IndelRealigner。
基本上，我只是删除以“@RG”开头且不是 ID e2_20 的行。
问题是，如果我在 bash 脚本中运行此命令，则该命令不起作用。
脚本如下。
#!/bin/bash
#$ -j y
#$ -S /bin/bash
#$ -V
#$ -cwd
source /apps1/modules/init/bash
module load samtools/gnu/1.1
input=$1
output=$2
samp=$3
samtools view -H ${input} | awk ' BEGIN {FS = "\t"} {split($2,a,":")} {if ($1 != "@RG" || ($1 =="@RG" && a[2] == "${samp}")) print}' | samtools reheader - ${input} > ${output}
echo "samtools view -H ${input} | awk ' BEGIN {FS = "\t"} {split($2,a,":")} {if ($1 != "@RG" || ($1 =="@RG" && a[2] == "${samp}")) print}' | samtools reheader - ${input} > ${output}"
该脚本在 SGE 集群中运行，这就是脚本顶部出现奇怪语法的原因。脚本的名称是 reHeadBams.bash。我通过输入以下命令在 shell 上运行脚本：
qsub reHeadBams.bash e2_20.indel.recal.dedup.bam e2_20.prac.indel.recal.dedup.bam e2_20
该命令的参数是输入文件，然后是输出文件，最后是我要查找的样本或读取组。
脚本的输出如下所示：
@HD VN:1.4 GO:none SO:坐标 @SQ SN:chr10 LN:129993255 @SQ SN:chr11 LN:121843856
@PG ID:GATK IndelRealigner。
因此脚本删除了所有读取组，而不是 ID 为 e2_9 的组。
我从脚本中回显了命令，输出是
samtools view -H e2_20.indel.recal.dedup.bam | awk ' BEGIN {FS = '\t'} {split(e2_20.prac.indel.recal.dedup.bam,a,:)} {if (e2_20.indel.recal.dedup.bam != '@RG' || (e2_20.indel.recal.dedup.bam =='@RG' && a[2] == 'e2_20')) print}' | samtools reheader - e2_20.indel.recal.dedup.bam > e2_20.prac.indel.recal.dedup.bam
虽然我可能弄错了，但问题似乎是 awk 使用我的命令行参数而不是输入文件的列作为 $1 和 $2。有谁知道为什么会这样？
很抱歉对我的问题进行了冗长的描述。如果您需要任何澄清，请告诉我。