Question: Bash script problem
0
gravatar for valastesolutions
6 weeks ago by
valastesolutions0 wrote:

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 • 152 views
ADD COMMENTlink modified 6 weeks ago by lieven.sterck8.0k • written 6 weeks ago by valastesolutions0
1

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 REPLYlink written 6 weeks ago by Carambakaracho2.2k

Can you post the script?

ADD REPLYlink written 6 weeks ago by ATpoint36k
#!/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 REPLYlink modified 6 weeks ago by Jean-Karim Heriche22k • written 6 weeks ago by valastesolutions0

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 REPLYlink modified 6 weeks ago • written 6 weeks ago by Carambakaracho2.2k

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 REPLYlink written 6 weeks ago by valastesolutions0

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 REPLYlink written 6 weeks ago by valastesolutions0
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1391 users visited in the last hour