Entering edit mode
8 months ago
Emmanuel
•
0
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?
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
Oh! I'm really sorry for that and I promise I won't repeat the same
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.Oh! I appreciate you for that. Well the contents of the sample_ids.txt are; SRR19170508, SRR19170510, SRR19170512, SRR19170514