Bash script problem
0
0
Entering edit mode
18 months ago

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 Bash Bedtools • 488 views ADD COMMENT 1 Entering edit mode 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... ADD REPLY 0 Entering edit mode Can you post the script? ADD REPLY 0 Entering edit mode #!/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

0
Entering edit mode

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}

0
Entering edit mode

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.

0
Entering edit mode

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