Question: How to pass arguments in a bash script to perform read alignment to the reference?
0
gravatar for Karma
9 months ago by
Karma300
India
Karma300 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 • 352 views
ADD COMMENTlink modified 9 months ago • written 9 months ago by Karma300
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 9 months ago by Pierre Lindenbaum131k

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 9 months ago by ATpoint42k

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

ADD REPLYlink written 9 months ago by Karma300
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 9 months ago • written 9 months ago by _r_am31k

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 9 months ago by lieven.sterck9.1k
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 9 months ago • written 9 months ago by lieven.sterck9.1k

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 9 months ago • written 9 months ago by Karma300

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 9 months ago by ATpoint42k
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 9 months ago by Joe18k

I managed to run this by running

sudo chmod 777 exome_pipeline.sh

and executed the scrpit sh exome_pipeline.sh

ADD REPLYlink modified 9 months ago • written 9 months ago by Karma300
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 9 months ago • written 9 months ago by Joe18k
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 9 months ago by _r_am31k
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: 1104 users visited in the last hour