max_coverage issue for running IsoQuant

[skkjy]

Hi, Andrey

I am graduate student majoring in bioinformatics. I previously questioned regarding the use of IsoQuant during 2 samples (normal, disease paired) test. I am reaching out again while scaling up the analysis to a larger cohort. (I am using ONT long-read sequencing data)

In this current run, I attempted to process 21 samples simultaneously using a merged GTF file and 21 corresponding BAM files. Due to the high depth of the merged data, the coverage at certain loci exceeded tens of millions of reads (in many loci), leading to a drastic increase in computational time. I think this process seems infeasible, I re-attempted the analysis by setting max_coverage_normal_chr to 5,000,000, while keeping the chrM coverage cutoff at its default value. Despite these adjustments, the process still took 17 days to complete.

I suspect that the primary cause is the amplification of read counts at specific chromosomal positions when multiple samples are processed together. However, I am concerned that imposing a coverage cutoff might result in the loss of quantitative information, potentially underestimating the expression levels of highly expressed genes or transcripts in the final output matrix.

To address these issues, I would like to ask for your expert opinion on the following:

My opinion (question): Would it be a feasible alternative to run IsoQuant on each sample individually, and then merge the result matrix files in Seurat for downstream processing and batch correction?

Related to Cutoff: Do you think the results obtained using the max_coverage_normal_chr 5000000 cutoff are reliable for downstream analysis, or would this significantly bias the quantification?

Alternative Solutions: Are there any other recommended strategies to optimize performance for a 21 samples cohort without compromising data integrity?

Thanks for reading my question.

[andrewprzh]

Dear @skkjy

I re-attempted the analysis by setting max_coverage_normal_chr to 5,000,000, while keeping the chrM coverage cutoff at its default value. Despite these adjustments, the process still took 17 days to complete.

Could you send me the log file please? I wonder which part took the longest time.

My opinion (question): Would it be a feasible alternative to run IsoQuant on each sample individually, and then merge the result matrix files in Seurat for downstream processing and batch correction?

If you use the same GTF file as an input, then merging counts should not be the issue.

Related to Cutoff: Do you think the results obtained using the max_coverage_normal_chr 5000000 cutoff are reliable for downstream analysis, or would this significantly bias the quantification?

Depends on your goal. If you compare expression between samples (more or less common use case), then the read count proportions of the same gene will remain the same, even if the cutoff is applied (since reads are skipped randomly and uniformly in all samples, if samples are processed together). Same goes to isoform counts and computing differential isoform expression within a gene.
However, proportions between different genes may become biased and so can overall quantification.

Note that processing samples separately and using coverage cutoff a the same time may lead to unusable and biased results.

Alternative Solutions: Are there any other recommended strategies to optimize performance for a 21 samples cohort without compromising data integrity?

You can also try to equally downsample all files (e.g. take 20% from each file), run the analysis, and later verify individual genes using full data if needed.
This might have a negative effect on detecting low-expressed isoforms, but should not affect the whole picture in general.

Best
Andrey

[skkjy]

Dear Andrey,

Thank you for your detailed explanation. Regarding the log file, I am currently working within a high-security internal network, which makes it difficult to export files externally. I am sorry about writing directly an example I was encountering.

WARNING - Genomic region 48962049-48994816 on chromosome chr19 (region legnth 32767, chromosome length 58617616) has coverage 13020347, which exceed coverage cutoff 5000000

As this log indicates, the cumulative coverage from the 21 pooled samples frequently exceeds the 5 million threshold.

Since my analysis aims to compare gene expression within a single sample (comparing gene expression between different cell types) as well as between condition groups (Normal vs Disease). So it is important to maintaining the integrity of inter-gene quantification.

Note that processing samples separately and using coverage cutoff a the same time may lead to unusable and biased results.

I am considering setting both --max_coverage_small_chr and --max_coverage_normal_chr to 10,000,000 instead of their default values. Since I will be processing each sample individually, I expect the coverage to stay well below this 10 million threshold. Do you think this would be an appropriate way to ensure data integrity and prevent any potential bias issues?

I have decided to adopt your suggestion of processing each sample individually.

• 2 Steps:

1. Run IsoQuant for each sample separately using the merged-GTF(21 samples') as a consistent reference and the individual sample BAM files. (setting both --max_coverage_small_chr and --max_coverage_normal_chr to 10,000,000)
   --> I expect to stay well below the coverage threshold for each run.
2. I will merge the resulting count matrices in Seurat for downstream analysis and batch correction.

I will try this approach first and, if performance issues persist, I will then consider the 20% downsampling strategy you suggested.

Thank you again for your guidance.

[andrewprzh]

Dear @skkjy

Do you think this would be an appropriate way to ensure data integrity and prevent any potential bias issues?

Yes, it makes sense. But if you expect the coverage to be lower that this threshold, you can set both parameters to -1 (unlimited) and see if it goes though.

Also, since you already have a GTF, use --no_model_construction, this might save a lot of time. Read assignment and quantification time should proportional to the number of reads, but model construction may have higher complexity in some cases.

Best
Andrey

[skkjy]

Dear @andrewprzh

Thank you for your kind answers to my repetitive questions. My apologies for the late reply as there was a holiday period here.

I’ve noted your suggestion about --no_model_construction option. As I intend to use the merged GTF for my upcoming steps, this may be suitable to my workflow.

I truly appreciate your dedication to maintaining the tool and providing helpful support.

Thanks Andrey