Empty seqz.gz output from bam2seqz

[PRAKHAR SONIK]

HI all,

I am trying to get the seqz.gz output from bam2seqz feature of Sequenza however, my output file comes out to be empty totally. 

I checked the order of the file it seems to be same as well and no progress.

this is the code I use:

sequenza-utils bam2seqz \
  -n Normal_Aligned.out.markdup.recal.bam \
  -t Tumor_Aligned.out.markdup.recal.bam \
  -gc hg38.gc50Base.ordered.wig.gz \
  -F hg38_sorted.fa \
  -o out.seqz.gz \
  -T usr/bin/tabix

[sming liu]

I have tried version 2.1.9999b0 and 3.0, both obtained empty results. I finally found an alternative solution: call mutation based on the paired BAMs, then use snp2seqz to obtain the sample seqz file.

[favero.f...@gmail.com]

Yes, that’s the “new” way I think sequenza files should be generated.

I’ve added new features in the R package and python package that are soon to be out.

I use strelka germline calls (generated using both tumor and normal bams) in ~2 minutes (vs days) in wgs samples

from iPhone

Il giorno 24 ott 2023, alle ore 05:31, sming liu <smli...@gmail.com> ha scritto:

I have tried version 2.1.9999b0 and 3.0, both obtained empty results. I finally found an alternative solution: call mutation based on the paired BAMs, then use snp2seqz to obtain the sample seqz file.

--
You received this message because you are subscribed to the Google Groups "Sequenza User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sequenza-user-g...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sequenza-user-group/b194cc96-2826-4159-bd1a-1bf723ee0f30n%40googlegroups.com.

[sming liu]

Woh, Great! 

I'm now using somatic mutations generated by Mutect2, should I turn to strelka for better performance? I'm also curious about why bam2seqz generates an empty file but the vcf based on the same BAMs not.

[Francesco Favero]

In theory whatever caller it would do. It’s supposed to have an interface to define which tag to parse in a vcf (but practically I am not sure it covers 100% of the vcf flavours).  

We use strelka because it just works well and it’s stable, but to be honest I’m concerned that’s not really open source and it’s not being updated in years.

Illumina is moving toward Dragen and I doubt they’ll update strelka any time soon.

So something like gatk haplotype caller or freebayes should also work.

I unfortunately have limited time to spend on developing this, but the concept is to use multiple  vcf files to form a seqz: 

    process the germline vcf calls, keeping the het positions only (I’ve added a parameter in sequenza-utils), obtaining a ~30-50Mb seqz

    process the somatic vcf calls keeping the pass variants (indels and SNVs) from the variant caller obtaining a tiny seqz file.

        Use sequenza utils to merge the two files

Other efforts are going to add the subclonality calls, I have various scripts that do this since almost 10 years now, and I never managed to clean the process and bundle it in the package.

So basically, the new sequenza will be faster to run, generating cleaner calls (also a colleague of mine is working on a new segmentation algorithm that I want to bundle in, replacing the deprecated copy number package) and annotating mutations and copy number with its relative cancel cell fraction (clonal/subclonal) 

I hope to be able to deliver this soon :D

To view this discussion on the web visit https://groups.google.com/d/msgid/sequenza-user-group/267fed74-dff4-42db-85c3-889091611993n%40googlegroups.com.

[Francesco Favero]

I have never experience the empty seqz file situation, I mean I did, but it was always for some human error :D

When the bam2seqz is run on the “full” bams, it relies on the order of the various sequences in the 2 bams, and in the GC wiggle file.

If for some reason the chromosome orders are so different that the 3 files go completely out of sync, it could generate an empty file.

For this I run bam2seqz in parallel for each chromosome, and merge the results later.

It could also be the “chr1” vs “1” situation. If one of the 3 file uses a different one, the sequenza utils is not checking for this, and it will give no results.

Mutect calls usually are made from a very standardised pipeline, eg using gatk best practices for alignment, variant calling, and using reference file from broad directly.

On 24 Oct 2023, at 10.19, sming liu <smli...@gmail.com> wrote:

To view this discussion on the web visit https://groups.google.com/d/msgid/sequenza-user-group/267fed74-dff4-42db-85c3-889091611993n%40googlegroups.com.

[sming liu]

Thanks for your kind reply.

The BAMs should be in the same order since they were generated using the standardized GTAK Somatic mutation pipeline. BAMs and GC wiggle file chromosomes are all in 'chr' format. I also tried to to run bam2seqz using --parallel and -C params, still obtained empty files for given chromosomes, unfortunately. Something must be turned wrong, but I don't know :(

[sming liu]

Only using  mutect2 mutations to generate seqz file may not be a good choice, this will generate a small file and may lead to an error while load data using sequenza.extract:

Processing chr1: Error in data.frame(base.count = as.integer(n.base.mut), maj.base.freq = as.numeric(max.freqs[,  :
  arguments imply differing number of rows: 196, 188

This can be soloved by using samtools mpileup (produce pielup file from BAM) + sequenza-utils bam2seqz (produce seqz file based on pielup file) + sequenza-utils seqz_binning (binning seqz), the final binned seqz file can be loaded sucessfully and produce usable results.

[Francesco Favero]

You can have a look at how we generate the seqz from variants here

https://bitbucket.org/weischenfeldt/pype_modules/src/master/variants/snippets/variants_vcf2seqz.md

It’s still not perfect (it could have issues on chrX in males) the key element is to generate germline calls using the tumor and the normal bams together (we use strelka)

Best

Francesco