RNA-seq QC and alignment error in script

Question

I am analyzing bulk RNA-seq data for Paired-End. I have separate scripts for fastqc, STAR & qualimap but want to run them in a single script which looks like this

USAGE: sh rna-qc.sh <path/to/sample_file/sample_R1.fastq>

The script takes as input the R1 fastq file ONLY

#!/bin/bash/

# initialize a variable to store the name of the input fastq file for Read1
fq1=$1

# grab the path information to use for loading the Read2 fastq file
path=`temp=$( realpath "$fq1" ) && dirname "$temp"`

# get base of filename 
samplename=`basename $fq _1.fastq.gz`
echo "Starting analysis of sample $samplename"

# create a variable to store the read 2 file
fq2=${path}/${samplename}_2.fastq.gz

# directory with the genome index files + name of the gene annotation file
genome=~/star/genomeIndex/
gtf=~/star/

# create the whole path
mkdir -p results/fastqc/
mkdir -p results/STAR/
mkdir -p results/qualimap/

# set up output filenames and locations
fastqc_out=results/fastqc/
align_out=results/STAR/${samplename}
align_out_bam=results/STAR/${samplename}_Aligned.sortedByCoord.out.bam
qualimap_out=results/qualimap/${samplename}.qualimap

echo "Processing file $fq"

echo "Starting QC for $samplename"

# Run FastQC
fastqc -t 8 -o $fastqc_out $fq1 $fq2

# Run STAR
STAR --runThreadN 8 --genomeDir $genome --readFilesIn $fq1 $fq2 --outFileNamePrefix $align_out --outSAMtype BAM SortedByCoordinate --outSAMunmapped Within --outSAMattributes Standard

# Run Qualimap
../qualimap_v2.3/qualimap rnaseq \
-outdir $qualimap_out \
-a proportional \
-bam $align_out_bam \
-gtf $gtf \
-pe
--java-mem-size=16G

This script named rna-qc.sh is then run using following loop

for fq in ~/rna_seq/raw_data/*_1.fastq.gz
do
rna-qc.sh $fq
done

then following is shown:

Starting analysis of sample _1.fastq.gz
Processing file
Starting QC for _1.fastq.gz
Skipping '/home/saad/rna_seq/raw_data/_1.fastq.gz_2.fastq.gz' which didn't exist, or couldn't be read
Started analysis of ERR4231010_1.fastq.gz
Approx 5% complete for ERR4231010_1.fastq.gz
.......

so on and only fastQC analysis completes on files. Please guide where could be issue in naming?

Also I want to know what argument should be used for qualimap of PE samples? Thanks!

Answer 1

#a bash script

#!/bin/bash


# Set your file paths and variables
raw_data_dir="~/rna_seq/raw_data"
output_dir="~/rna_seq/results"
genome="path/to/genome_index"
qualimap_dir="../qualimap_v2.3"  # Adjust as needed

# Create output directories
mkdir -p "$output_dir/fastqc"
mkdir -p "$output_dir/star"
mkdir -p "$output_dir/qualimap"

# Loop through all paired-end FASTQ files
for fq1 in "$raw_data_dir"/*_1.fastq.gz; do
    fq2="${fq1/_1.fastq.gz/_2.fastq.gz}"  # Construct the second read file name

    # Run FastQC
    fastqc -t 8 -o "$output_dir/fastqc" "$fq1" "$fq2"

    # Run STAR
    sample_name=$(basename "$fq1" _1.fastq.gz)
    align_out="$output_dir/star/${sample_name}_"
    STAR --runThreadN 8 --genomeDir "$genome" --readFilesIn "$fq1" "$fq2" --outFileNamePrefix "$align_out" --outSAMtype BAM SortedByCoordinate --outSAMunmapped Within --outSAMattributes Standard

    # Run Qualimap
    qualimap_out="$output_dir/qualimap/${sample_name}"
    "$qualimap_dir/qualimap" rnaseq -outdir "$qualimap_out" -a proportional -bam "${align_out}Aligned.sortedByCoord.out.bam" -gtf "$gtf" --java-mem-size=16G
done