Erreur récurrente avec samtools view

Bonjour,

Je rencontre des erreurs lors de la conversion avec samtools de mes fichiers .sam en fichier .bam et ne comprends pas d'où elle viennent.

Les md5sum de mes fichiers de séquençage sont ok. J'ai également fait le fastqc de mes fichiers et mes banques sont de bonnes qualités.

Concernant l'alignement, j'ai exécuté le script suivant :

#!/bin/bash
#Alignement

#SBATCH --job-name=alignement_AM-36
#SBATCH --mem=40GB
#SBATCH --nodes=2
#SBATCH --cpus-per-task=8

# le script va s'arrêter
# - à la première erreur
# - si une variable n'est pas définie
# - si une erreur est recontrée dans un pipe
set -euo pipefail


#Function of the script: aligning sequences obtained from paired-end sequencing

##Loading of modules##
module load bowtie2/2.5.1


#Alignement AM_19
##Paths##
REF_GENOME="/shared/projects/organoid/genome/hg38/bowtie2/hg38" #Path of the reference genome (hg38)
FASTQC_FILES="/shared/projects/organoid/data/organoid/AM_36/banque" #Path of the directory containing all the data to analyze
ALIGN_DIR="/shared/projects/organoid/data/organoid/AM_36/alignements_36" #Path of the directory where you want to store your results

#Creation of a function that will do the alignment for a given lane and give an output corresponding to the results of the alignment

echo "=============================================================="
echo "Aligner les reads sur le génome de référence : échantillon 36 S1"
echo "=============================================================="

srun bowtie2 -p "${SLURM_CPUS_PER_TASK}" -x ${REF_GENOME} -1 ${FASTQC_FILES}/AM-36-NT_S1_R1_cat.fastq.gz -2 ${FASTQC_FILES}/AM-36-NT_S1_R2_cat.fastq.gz -S ${ALIGN_DIR}/alignment_AM_36_S1.sam 2> ${ALIGN_DIR}/outputInput_AM_36_S1.txt


echo "=============================================================="
echo "Alignement fini"
echo "=============================================================="

L'alignement se déroule a priori bien et j'obtiens un bon taux d'alignement (autour de 90%).

Pour la conversion de mes fichiers .sam en fichiers .bam j'utilise le script suivant :

#!/bin/bash

#SBATCH --job-name=conversion
#SBATCH --mem=40GB
#SBATCH --cpus-per-task=8

# le script va s'arrêter
# - à la première erreur
# - si une variable n'est pas définie
# - si une erreur est recontrée dans un pipe
set -euo pipefail

##Loading of modules##


module load samtools/1.18 #version 1.18 utilisée (la plus récente dispo)

# En se plaçant dans le répertoire /shared/projects/organoid/data/organoid/AM_36/
cd /shared/projects/organoid/data/organoid/AM_36/

##Conversion##
srun samtools view -f 0x2 -b -@ 5 alignements_36/alignment_AM_36_S1.sam > alignements_36_bam/alignment_AM_36_S1.bam
echo "Conversion AM_36 : done"

##Sorting##
srun samtools sort -@ 5 alignements_36_bam/alignment_AM_36_S1.bam -o alignements_36_bam/alignment_AM_36_S1_sorted.bam
echo "Sorting AM_36 : done"

##Indexing##
srun samtools index -@ 5 alignements_36_bam/alignment_AM_36_S1_sorted.bam
echo "Index AM_36 : done"

Cependant, quelques secondes après le lancement du script, celui-ci se stoppe et j'obtiens cette erreur :

[E::sam_parse1] CIGAR and query sequence are of different length
samtools view: error reading file "alignements_36/alignment_AM_36_S1.sam": Input/output error
samtools view: error closing "alignements_36/alignment_AM_36_S1.sam": -5
srun: error: cpu-node-61: task 0: Exited with exit code 1

J'avais déjà réalisé ces mêmes instructions sur d'autres données de séquençage et n'avais jamais rencontré le moindre problème.

Il me reste 1T d'espace, je ne pense donc pas que ce soit lié à un problème d'espace.

En outre, j'ai une erreur quelle que soit la version de samtools que j'utilise.

J'ai également relancé plusieurs fois mes alignements avec bowtie2 et obtiens à chaque fois un message d'erreur à cette étape.

J'ai également lancé la commande sur le terminal :

samtools view -f 0x2 -b -@ 5 alignements_36/alignment_AM_36_S1.sam > alignements_36_bam/alignment_AM_36_S1.bam

et j'obtiens toujours une erreur, qui peut varier (bien que la commande reste la même) :

[E::sam_parse1] CIGAR and query sequence are of different length
[E::parse_cigar] CIGAR length invalid at position 1 (TTTCTTTCTTTTTTTTTTTCCCCCGAGATGAAGTCTCGCTCTGTCACCCAGGCTAGAGTGCAGTGGTGTGATCTCA AAAAAEEEEEEEEEEEEEEEAAEAEE/EEEEEEEAEEEEEEEEEAEEEAAAE<EEEEEEAAEEEEEEEEAEEEEA< AS:i:-11 XN:i:0 XM:i:0 XO:i:1 XG:i:2 NM:i:2 MD:Z:74 YS:i:-5 YT:Z:CP)
[E::aux_parse] unrecognized type 'E'
samtools view: error reading file "alignements_36/alignment_AM_36_S1.sam": Input/output error
samtools view: error closing "alignements_36/alignment_AM_36_S1.sam": -5

ou encore :

[E::sam_parse1] CIGAR and query sequence are of different length
samtools view: error reading file "alignements_36/alignment_AM_36_S1.sam": Input/output error
samtools view: error closing "alignements_36/alignment_AM_36_S1.sam": -5

ou encore :

[E::sam_parse1] CIGAR and query sequence are of different length
[E::aux_parse] unrecognized type 'E'
[E::sam_parse1] CIGAR and query sequence are of different length
[E::sam_parse1] SEQ and QUAL are of different length
[E::sam_parse1] SEQ and QUAL are of different length
[E::sam_parse1] SEQ and QUAL are of different length
[E::sam_parse1] [E::parse_cigar] CIGAR length invalid at position 1 (TTTCTTTCTTTTTTTTTTTCCCCCGAGATGAAGTCTCGCTCTGTCACCCAGGCTAGAGTGCAGTGGTGTGATCTCA AAAAAEEEEEEEEEEEEEEEAAEAEE/EEEEEEEAEEEEEEEEEAEEEAAAE<EEEEEEAAEEEEEEEEAEEEEA< AS:i:-11XN:i:0 XM:i:0 XO:i:1 XG:i:2 NM:i:2 MD:Z:74 YS:i:-5 YT:Z:CP)
SEQ and QUAL are of different length
[E::sam_parse1] CIGAR and query sequence are of different length
samtools view: error reading file "alignements_36/alignment_AM_36_S1.sam": Input/output error
samtools view: error closing "alignements_36/alignment_AM_36_S1.sam": -5

Également, ces mêmes commandes ont été lancées, et la conversion s'est faite sans problème avec samtools view (samtools/1.18) par une autre personne travaillant sur un autre cluster que celui de l'IFB.

Je voulais alors savoir si ce problème avait déjà été rencontré et à quoi pouvait-il être dû.

Je vous remercie d'avance pour votre aide,

Aurélie

Bonjour @amasson, d'après certains sujets sur biostars il semblerait que ce problème viendrait plutôt de l'alignement en lui-même et pas de samtools. Est-ce que vous pouvez essayer en alignant avec bwa (module bwa-mem2) ?

Ici CIGAR and query sequence are of different length when trying to convert from sam to bam? une personne suggère d'utiliser picard tools pour valider l'alignement, ca peut peut être aider également