Bash script problem

Question

[valastesolutions]

Hi,

I am writing a bash script to call variants in RNA-seq data. At some point, I want to select only ready reads corresponding to a specific genomic region using bedtools intersect and redirect the output to a smaller BAM file more amenable for variant calling. I have used following command for this:

bedtools intersect -wa -a $1 -b $2 > KRAS_s$1

with $1 and $2 taking BAM and BED files from standard input respectively. However, this does not seem to work properly as the output file is created but contains no reads. Running this command outside the bash script delivers the desired result, so the problem does not seem to be software related. I have done some reading, and think the problem relates to >, but I am not sure how to fix it. Can anyone help? Thanks

Steven

[Carambakaracho]

Hi Steven, there's some additional info needed to support you:

Is this the full body for the script? How do you run the bash script? Can you echo $1 and echo $2? Does this yield what you expect?

I don't think > should be the problem...

[ATpoint]

Can you post the script?

[valastesolutions]

#!/bin/bash

echo initiating callvar_rna on $1

echo running samtools sort

samtools sort $1 -o s$1

echo running samtools index

samtools index s$1

echo selecting reads from

bedtools intersect -wa -a s$1 -b $2 > KRAS_s$1

[Carambakaracho]

Just to make sure you're not hunting a ghost: the sorted bam in s$1 contains all reads?

I don't see any obvious show stopper, so maybe your input interferes with the call. You might want to try parameter substitution

bedtools intersect -wa -a s${1} -b ${2} > KRAS_s${1}

[valastesolutions]

Thanks, I will try the parameter substitution. I am absolutely sure that the BAM file in s$1 contains all reads. I will let you know if this works. Thanks all for your time.

[valastesolutions]

Hi,

in fact I build the script first echoing each line to see if the commands are correct, and in fact they are. I have posted the complete script in replay to the remark below. Any help would be welcome.

Steven