Problem while working with sequenza - Chromosomes out of order

Question

[Camilo Andres]

Hi, I'm trying to work with sequenza in order to calculate HRD score of a sample using WES data. When I run sequenza, I get a message saying that "chromosomes are out of order", and I don't know how to solve it, so if you could help me I would really appreciate it.

This is the code I'm using:

sequenza-utils bam2seqz -gc hg38/hg38.wig.gz --fasta hg38/hg38.fa -n MCF10A/SRR3090727_recall_reads.bam -t MDA-MB-231/SRR1021654_recall_reads.bam -C chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22 chr23 chr24 chrX --output sample.seqz.gz
`

and this is the output I get:

[mpileup] 1 samples in 1 input files
[mpileup] 1 samples in 1 input files
[E::bam_plp_push] [E::bam_plp_push] The input is not sorted (chromosomes out of order)The input is not sorted (chromosomes out of order)

samtools mpileup: error reading from input file
samtools mpileup: error reading from input file

The bam files were sorted using picard's SortSam tool, so I don't really understand what is happening.

Thanks in advance!

[Camilo Andres]

when I run

samtools view -H MDA-MB-231/SRR1021654_recall_reads.bam | head

I get

@HD VN:1.6  SO:coordinate
@SQ SN:chr1 LN:248956422
@SQ SN:chr10    LN:133797422
@SQ SN:chr10_GL383545v1_alt LN:179254
@SQ SN:chr10_GL383546v1_alt LN:309802
@SQ SN:chr10_KI270824v1_alt LN:181496
@SQ SN:chr10_KI270825v1_alt LN:188315
@SQ SN:chr10_KN196480v1_fix LN:277797
@SQ SN:chr10_KN538365v1_fix LN:14347
@SQ SN:chr10_KN538366v1_fix LN:85284

and tail

@SQ SN:chrY LN:57227415
@SQ SN:chrY_KI270740v1_random   LN:37240
@SQ SN:chrY_KN196487v1_fix  LN:101150
@SQ SN:chrY_KZ208923v1_fix  LN:48370
@SQ SN:chrY_KZ208924v1_fix  LN:209722
@SQ SN:chrY_MU273398v1_fix  LN:865743
@RG ID:sample   LB:Paired-end   PL:Illumina SM:sample   PU:Unknown
@PG ID:bwa  PN:bwa  VN:0.7.17-r1188 CL:bwa mem -t 6 -M hg38/hg38.fa MDA-MB-231/SRR1021654_1.fastq.gz MDA-MB-231/SRR1021654_2.fastq.gz
@PG ID:GATK ApplyBQSR   VN:4.4.0.0  CL:ApplyBQSR --output MDA-MB-231/SRR1021654_recall_reads.bam --bqsr-recal-file MDA-MB-231/SRR1021654_recall_data.table --input MDA-MB-231/SRR1021654_sorted_RG.bam --reference hg38/hg38.fa --preserve-qscores-less-than 6 --use-original-qualities false --quantize-quals 0 --round-down-quantized false --emit-original-quals false --global-qscore-prior -1.0 --interval-set-rule UNION --interval-padding 0 --interval-exclusion-padding 0 --interval-merging-rule ALL --read-validation-stringency SILENT --seconds-between-progress-updates 10.0 --disable-sequence-dictionary-validation false --create-output-bam-index true --create-output-bam-md5 false --create-output-variant-index true --create-output-variant-md5 false --max-variants-per-shard 0 --lenient false --add-output-sam-program-record true --add-output-vcf-command-line true --cloud-prefetch-buffer 40 --cloud-index-prefetch-buffer -1 --disable-bam-index-caching false --sites-only-vcf-output false --help false --version false --showHidden false --verbosity INFO --QUIET false --use-jdk-deflater false --use-jdk-inflater false --gcs-max-retries 20 --gcs-project-for-requester-pays  --disable-tool-default-read-filters false    PN:GATK ApplyBQSR
@PG ID:samtools PN:samtools PP:GATK ApplyBQSR   VN:1.17 CL:samtools view -H MDA-MB-231/SRR1021654_recall_reads.bam
`

[Istvan Albert]

This error seems to come from the samtools mpileup, run that command on each bam file to troubleshoot. Maybe resort?

Answer 1

[Camilo Andres]

The error was that the chromosomes were not sequentially ordered, so what I did was to order the fasta file of the reference genome, and then create all the files related (index and others). Once this was done, I did the alignment again and all the posterior steps and it worked.

[Istvan Albert]

thanks for following up with an answer