BQ of generated reads do not match model specification

[cntrier]

Hi,

I've created a read model with the following script:

mitty create-read-model synth-illumina
100.pkl
--read-length 100
--mean-template-length 250
--std-template-length 20
--bq0 30
--k 200
--sigma 5

And when I check the read model in Mitty with the following:
mitty describe-read-model 100.pkl 100.png

It looks as expected:

But when I generate reads using the model with the following code:
k=HG00632
i=100

mitty -v4 generate-reads GRCh38.p12.fa
./final_vcfs/${k}all.vcf.gz
${k} all_merged_sorted.bed
${i}.pkl
40
7
${k}${i}reads-test.1.fq
${k}${i}-lq.txt
--fastq2 ${k}_${i}reads-test.2.fq
2> vcf-${i}${k}.log

The generated reads have a flat BQ of 9 when I check them with FastQC:

And when I run the god-aligner to create a bam file, I can see in IGV that the reads are a mess. I've tried running different individuals, different read lengths but get the same pattern.

Have I misunderstood something with the read model generation?

Thank you very much for any help you can provide on the matter.

[ghost]

@cntrier thank you for this comprehensive report! Will make it easy to replicate and hopefully isolate and fix.

[ghost]

@cntrier if you use the parameters given in the examples in the Readme, can you replicate the results shown in the Readme? This information is useful to me as a means of understanding if the bug has always been there, or is a recently introduced bug. Thanks!

[cntrier]

Thank you for getting back to me.

I have now downloaded the demo data and ran the same parameters as in the README and encountered the same problem.

I filtered variants and generated a synthetic read model with the code from the README and I have the same output as the README with mitty describe-read-model function:

When I generated reads using the synthetic read model with the parameters from the README, again the reads had a flat BQ of 9.

Here is the code:

FASTA=human_g1k_v37.fa.gz
SAMPLEVCF=1kg.20.22.vcf.gz
SAMPLENAME=HG00119
REGION_BED=region.bed
FILTVCF=HG00119-filt.vcf.gz
READMODEL=1kg-pcr-free.pkl
COVERAGE=30
READ_GEN_SEED=7
FASTQ_PREFIX=HG00119-reads
READ_CORRUPT_SEED=8

docker run --rm -v /Users/cassandra/Desktop/Ous_mitty:/VCF:rw ous_nbs/mitty /opt/conda/envs/mymitty/bin/mitty -v4 filter-variants \
  /VCF/${SAMPLEVCF} \
  ${SAMPLENAME} \
  /VCF/${REGION_BED} \
  -  \
  2> vcf-filter.log | bgzip -c > ${FILTVCF}
 
tabix -p vcf ${FILTVCF}

MODELNAME=synthetic-model.pkl

docker run --rm -v /Users/cassandra/Desktop/Ous_mitty:/VCF:rw ous_nbs/mitty /opt/conda/envs/mymitty/bin/mitty create-read-model synth-illumina \
  /VCF/${MODELNAME} \
  --read-length 121 \
  --mean-template-length 400 \
  --std-template-length 20 \
  --bq0 30 \
  --k 200 \
  --sigma 5 \
  --comment 'Model created for the demo'

docker run --rm -v /Users/cassandra/Desktop/Ous_mitty:/VCF:rw ous_nbs/mitty /opt/conda/envs/mymitty/bin/mitty describe-read-model /VCF/${MODELNAME} /VCF/${MODELNAME}.png
  
docker run --rm -v /Users/cassandra/Desktop/Ous_mitty:/VCF:rw ous_nbs/mitty /opt/conda/envs/mymitty/bin/mitty -v4 generate-reads \
  /VCF/${FASTA} \
  /VCF/${FILTVCF} \
  ${SAMPLENAME} \
  /VCF/${REGION_BED} \
  /VCF/${MODELNAME} \
  ${COVERAGE} \
  ${READ_GEN_SEED} \
  /VCF/${FASTQ_PREFIX}1.fq \
  /VCF/${FASTQ_PREFIX}-lq.txt \
  --fastq2 /VCF/${FASTQ_PREFIX}2.fq \
  --threads 2

Thank you for your time!

[ghost]

@cntrier thank you so much for this report - this is super helpful! I will fix this bug.