BWA alignement for RNA Seq

Question

[storm1907]

I have to write script for bwa alignement with trimmed reads (not with samtools at the moment yet). I have a poor understanding about bwa code. 1) Bash does not recognize command "bwa", although I changed the tool directory before calling the command (export to PATH in the directory did not help either) 2) Does indexing for viral RNA can look like that: bwa index $reference? 3) what module should I load at the start of script? i.e., for Trimmomatics it was module load java

Thank You :)

#!/bin/bash -x
#PBS -N bwa samtools
#PBS -q batch

INPATH=/home/groups/dir/usr/subdir
OUTPATH=/home/groups/dir/usr/subdir
reference=/dir/subdir/refseq.fasta

cd /mnt/home/usr/tools/bwa-0.7.17
for file in $INPATH/*R1*.fastq.gz
do
        bname=$(basename $file 'R1_001.unpaired.fastq.gz')
        echo "file: "$file
        echo $bname
        input1=$INPATH/$bname"R1_001.unpaired.fastq.gz"
        output1=$OUTPATH/$bname"R1_001.unpaired.sam"
        input2=$INPATH/$bname"R2_001.unpaired.fastq.gz"
        output2=$OUTPATH/$bname"R2_001.unpaired.sam"
        outfile=$OUTPATH$bname"pe.sam"
        echo $input1 $output1 $input2 $output2 $outfile
        bwa index -p refseq $reference
        bwa mem aln -t 4 $reference $input1 > $output1
        bwa mem aln -t 4 $reference $input2 > $output2
        bwa sampe $reference $output1 $output2 $input1  $input2 > $outfile
done

Answer 1

[Nicolas Rosewick]

you don't need to reindex the reference transcriptome each time. You can put bwa index -p refseq $reference outside the for loop. Also why do you use bwa to map RNA-Seq against transcriptome. A better way would be to map the RNA-Seq reads against the reference genome using a splice-aware aligner (STAR, HiSAT2, .. ).

An other option if you do not have a reference genome for your organism of interest but a transcriptome file, you can try a pseudoalignment method e.g. kallisto or salmon

[storm1907]

Thanks! This time I am trying to align ss RNA virus, so maybe bwa could be suitable