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
•
written
by
storm1907 •
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
Thanks! This time I am trying to align ss RNA virus, so maybe bwa could be suitable
Log in to answer this question.
