How to load sam.gz files on IGV?

Question

[marongiu.luigi]

I have mapped reads to a reference genome and generated both a sam.gz and vcf files. I would like to visualize the mapping with the Integrative Genomics Viewer (IGV) and I loaded the reference genome in the viewer. When I load the sam.gz files, though, I get the error:

I can see from the Web that the standard file to use is BAM rather than SAM. I am now using Samtools to make the conversion. THe first step is to create an index, but I get the error:

$ samtools index -b V25_chr1_TrmDedMap.sam.gz 
[E::sam_index_build3] SAM file "HoVa25.sam.gz" not BGZF compressed
samtools index: failed to create index for "V25_chr1_TrmDedMap.sam.gz"

I followed this post to try and create a bam with index, but again I got an error:

even though I have these fules in the same folder:

How can I load them properly into IGV?

Thank you

[Arup Ghosh]

Share the output of file V25_chr1_TrmDedMap.sam.gz.

If it's a text file with gzip compression, use the following commands to convert it into a coordinate-sorted and indexed BAM file.

gunzip -c V25_chr1_TrmDedMap.sam.gz | samtools sort -@ 4 -O BAM -o V25_chr1_TrmDedMap.bam # using 4 threads
samtools index V25_chr1_TrmDedMap.bam

Answer 1

[Pierre Lindenbaum]

file must be compressed with BGZF compression as a BAM file, not just gzip.

assuming your file is gzipped AND sorted on coordinates

gunzip -c  Ori.sam.gz  | samtools view -O BAM -o Ori.bam &&  samtools index Ori.bam

[marongiu.luigi]

Thank you, I tried but I got the error:

$ gunzip -c V25_chr1_TrmDedMap.sam.gz | samtools view -O BAM -o HoVa25_chr1.bam && samtools index HoVa25_chr1.bam
[E::hts_idx_push] NO_COOR reads not in a single block at the end 0 -1
[E::sam_index] Read 'LH00201:89:22VF3NLT3:3:2271:18309:7024 1:N:0:CGTTAGCTGAAC+TGACTTCAGCAG' with ref_name='NZ_CP098033.1 Vibrio alginolyticus strain E110 chromosome 1, complete sequence', ref_length=1819821, flags=99, pos=1582024 cannot be indexed
samtools index: failed to create index for "HoVa25_chr1.bam"

[ATpoint]

Then the file is somehow corrupted. There is little point debuggung this here since the file format sam.gz itself is unstandard, suggesting custom fiddling or other oddities, so upstream processing was at least "not optimal". The person who created this bam should simply arrive in this millenium and produce a coordinate-sorted BAM file (aka aligner + samtools sort, the simplest operation in the world basically) as everyone else does, rather than this unstandard thingy stuff which causes exactly these sorts of problems.

[marongiu.luigi]

Thanks for the kind feedback. The file was created with BBmap from deduplicated files: bbmap.sh in=Trimmed/V2025_1_TrmDed.fq.gz in2=Trimmed/V2025_2_TrmDed.fq.gz out=Mapped/V25_chr1_TrmDedMap.sam.gz ref=Ref/V_chr1_ref.fa. Is there a better command to map these files with BBmap or shall I go for a standard bwa-mem?

[michael.ante]

If your fasta headers have descriptions , I'd add the option trimreaddescriptions=t to only have "NZ_CP098033.1" in the sam-file. I'd also just output the plain sam file and use Pierre's snippet with a cat instead of a gunzip -c.

[GenoMax]

You should have used

out=Mapped/V25_chr1_TrmDedMap.bam

instead of

out=Mapped/V25_chr1_TrmDedMap.sam.gz

BBMap is capable of writing BAM files when samtools or sambamba are available in $PATH.

You could also send the output to STDOUT and then sort/index the output file in one step

bbmap.sh in=Trimmed/V2025_1_TrmDed.fq.gz in2=Trimmed/V2025_2_TrmDed.fq.gz out=stdout.sam ref=Ref/V_chr1_ref.fa | samtools sort --write-index -o Mapped/V25_chr1_TrmDedMap.bam  -