Question: How to pass arguments in a bash script to perform read alignment to the reference?
0
gravatar for Karma
6 weeks ago by
Karma270
India
Karma270 wrote:

My plan to is to make automated pipeline for processing exome data to derive the vcf file. I coded the alignment step in the following way

 #!usr/bin/bash
    reference="/media/ngs/ref/Homo_sapiens_assembly19/Homo_sapiens_assembly19.fasta"
    outputDir=$1
    sampleID=$2
    trimmedR1=$3
    trimmedR2=$4
    exec &> $outputDir/$sampleID.StdOp.Err.txt; #Redirect the stdout and stderror to a text file

    echo "Starting Run ==" `date +%d/%m/%Y\ %H:%M:%S` 2>> $outputDir/$sampleID.StdOp.Err.txt

    # saving the readGroup to the variable
    readGroup="@RG\tID:${sampleID}\tSM:${sampleID}\tLB:${sampleID}\tPL:Illumina"

    echo "BWA mem runtime start" `date +%d/%m/%Y\ %H:%M:%S` 2>> $outputDir/$sampleID.StdOp.Err.txt
    # bwa mem algorithm is used for alignment
    bwa mem -M -t 16 -R "$readGroup" $reference -o $outputDir/$sampleID/${sampleID}_raw.sam $trimmedR1 $trimmedR2 2>> $outputDir/$sampleID.StdOp.Err.txt
    echo "BWA mem runtime end" `date +%d/%m/%Y\ %H:%M:%S` 2>> $outputDir/$sampleID.StdOp.Err.txt

and I ran the code as follows:

sh exome_pipeline.sh outputDir AB132 AB132_R1.fq.gz AB132_R2.fq.gz

It gives me an error

[main_mem] fail to open file '/AB132 /AB132_raw.sam' : No such file or directory

How can fix this error?

bash bwa exome data analysis • 163 views
ADD COMMENTlink modified 6 weeks ago • written 6 weeks ago by Karma270
2

to make automated pipeline for processing exome data to derive the vcf file

using plain bash ? you'd better learn how to manage a workflow with nextflow / snakemake...

ADD REPLYlink written 6 weeks ago by Pierre Lindenbaum127k

Agreed. If you are getting started with writing pipelines it is a bit of an extra effort to work yourself into a workflow manager but the effort is worth it. You also would need to invest effort into writing the bash script which should have some basic testing and exiting functions in case things go wrong so that you can easily debug. Use a workflow manger.

ADD REPLYlink written 6 weeks ago by ATpoint32k

Is it not possible to create pipeline with plain bash scripting?

ADD REPLYlink written 6 weeks ago by Karma270
1

I built a bash-based pipeline for exome sequencing a couple of years ago. From your script above, I gather that your shell experience is quite basic and nowhere close enough to be able to build a robust pipeline that accounts for various edge cases. Please invest time in learning a workflow manager, because you're going to need to invest effort in learning something anyway.

ADD REPLYlink modified 6 weeks ago • written 6 weeks ago by RamRS26k

yes it is, but as Pierre Lindenbaum and ATpoint pointed out, it's more future-proof if you go for one of the pipeline manager tools

ADD REPLYlink written 6 weeks ago by lieven.sterck7.2k
1

looks like it does not get the $outputDir correctly. and as such tries to write to the root folder / to which you might not have access

can you post all output from your script (and perhaps add an echo $outputDir in it)

ADD REPLYlink modified 6 weeks ago • written 6 weeks ago by lieven.sterck7.2k

I cross checked outputDir is available as directory. I just manually created a directory 'AB132' inside outputDir and executed the script again. I don't know it is working!!!. The directories are not in the root folder if it is related to any permission issues.

ADD REPLYlink modified 6 weeks ago • written 6 weeks ago by Karma270

That is not what lieven.sterck is saying. Check the error message.

fail to open file '/AB132 /AB132_raw.sam'

First, /AB132 is in fact referring to a root directory. So probably outputDir is wrong, so not correctly pared. Also there is a whitespace in that path which is probably messing up things.

ADD REPLYlink written 6 weeks ago by ATpoint32k
1

If that error is verbatim, you have an unescaped space character after the first /AB123. Unix will not play nicely with this.

ADD REPLYlink written 6 weeks ago by Joe16k

I managed to run this by running

sudo chmod 777 exome_pipeline.sh

and executed the scrpit sh exome_pipeline.sh

ADD REPLYlink modified 6 weeks ago • written 6 weeks ago by Karma270
1

Why are you using sudo here?

You've given that script admin privileges and that is not a safe thing to do.

You may have treated the symptom of your script failing, but not the cause.

ADD REPLYlink modified 6 weeks ago • written 6 weeks ago by Joe16k
1

You've given that script admin privileges

No, OP has changed permissions using a sudo command, which is meaningless in this context if OP owns the script file.

You may have treated the symptom of your script failing, but not the cause.

OP has treated neither. OP is not running sudo sh exome_pipeline.sh, just sudo chmod followed by sh exome_pipeline.sh.

ADD REPLYlink written 6 weeks ago by RamRS26k
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: 820 users visited in the last hour