Classical alignment with HISAT2

Question

[Emmanuel]

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 received 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?

[Ram]

You created multiple posts on multiple sites for the same question. That is bad etiquette and only serves to annoy people in both communities. Remember you're asking volunteers in two places to spend their time on you while not telling them that you've asked another set of volunteers as well. You're not getting quotes from competing businesses, you're wasting volunteers' time. Please do not repeat this in the future.

Link to cross-post: https://bioinformatics.stackexchange.com/questions/21987/classical-alignment-with-hisat2

[Emmanuel]

Oh! I'm really sorry for that and I promise I won't repeat the same

[GenoMax]

This error does not make sense since SRR19170514 is a paired-end library.

What did you put in the sample_ids.txt file when you ran this script?

BTW: I formatted your post with code 10101 button that you see in editor window. Please use that in future when posting code.

[Emmanuel]

Oh! I appreciate you for that. Well the contents of the sample_ids.txt are; SRR19170508, SRR19170510, SRR19170512, SRR19170514

Answer 1

[Ram]

This has been answered on bioinformatics SE: https://bioinformatics.stackexchange.com/a/21988/650

[GenoMax]

When using -p one needs to also add the following parameter which is not mentioned in the above answer. This is a critical parameter.

--countReadPairs    If specified, fragments (or templates) will be counted
                      instead of reads. This option is only applicable for
                      paired-end reads.