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
ADD REPLY
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}
ADD REPLY
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.

ADD REPLY
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

ADD REPLY

Login before adding your answer.

Traffic: 1115 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6