Classical alignment with HISAT2

Question

This question was also asked on Biostars

I have this script

#!/bin/bash

# Defining the number of threads to be used
threads=20

# Alignment of sequences to the reference genome using 'HISAT2'
echo -e "\nNow Running Classical Alignment-based with HISAT2...\n"

mkdir -p hisat2

# Download the human reference genome
wget -c ftp://ftp.ensembl.org/pub/release-100/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz -O hisat2/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
gunzip hisat2/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz

# Download the human annotation file
wget -c ftp://ftp.ensembl.org/pub/release-100/gtf/homo_sapiens/Homo_sapiens.GRCh38.100.gtf.gz -O hisat2/Homo_sapiens.GRCh38.100.gtf.gz
gunzip hisat2/Homo_sapiens.GRCh38.100.gtf.gz

# HISAT2 indexing to build the reference genome index
hisat2-build -p $threads hisat2/Homo_sapiens.GRCh38.dna.primary_assembly.fa hisat2/Homo_sapiens.GRCh38v3_hisat2.idx

# HISAT2 Alignment
for sample in $(cat sample_ids.txt)
do
    hisat2 -p $threads -x hisat2/Homo_sapiens.GRCh38v3_hisat2.idx \
        -1 trimmed_data/${sample}_1_val_1.fq -2 trimmed_data/${sample}_2_val_2.fq \
        -S hisat2/${sample}_hisat2.sam

    # Compress the SAM files to binary format, sort them, and index them.
    samtools view -Sb hisat2/${sample}_hisat2.sam | samtools sort > hisat2/${sample}_hisat2_sorted.bam
    samtools index hisat2/${sample}_hisat2_sorted.bam
rm hisat2/${sample}_hisat2.sam  # Remove the .sam files for storage purposes
done

# Quantifying Aligned reads using the Subread package's script 'featureCounts'
for sample in $(cat sample_ids.txt)
do
featureCounts -T $threads -a hisat2/Homo_sapiens.GRCh38.100.gtf \
        -o hisat2/${sample}_hisat2_counts.txt \
    hisat2/${sample}_hisat2_sorted.bam
done

# Generating analysis report with MultiQC
echo "Running MultiQC on HISAT2 alignment results ........."
multiqc -o multiQCreports hisat2
echo "MultiQC on HISAT2 alignment results completed successfully."

# Display message to show that the script is complete
echo "Script completed!"

However, I have receieved an error below

Load annotation file Homo_sapiens.GRCh38.100.gtf ...                       ||
||    Features : 1378275                                                      ||
||    Meta-features : 60683                                                   ||
||    Chromosomes/contigs : 47                                                ||
||                                                                            ||
|| Process BAM file SRR19170514_hisat2_sorted.bam...                          ||
ERROR: Paired-end reads were detected in single-end read library : hisat2/SRR19170514_hisat2_sorted.bam
ERROR: Paired-end reads were detected in single-end read library : hisat2/SRR19170514_hisat2_sorted.bam

Question: What can I do to solve this?

Answer 1

I believe the error is from featureCounts and the solution would be to add the -p flag to the featureCounts call, stating that you are providing paired-end data.

EDIT: You might want to check this r/bioinformatics (reddit) thread out.