Low recall and precision using v3.6.0

[LilyWang07]

I ran NanoCaller v3.6.0 on ONT NA24385 data using the default parameters:
NanoCaller --bam ${bam} --ref ${reference} --cpu 15

I also tried specifying the SNP and indel models:
NanoCaller --bam ${bam} --ref ${reference} --seq ont --snp_model ONT-HG002 --indel_model ONT-HG002 --sample HG002

I compared the calling results using hap.py and the hg19 reference. When using both sets of parameters, the evaluated indel recall and precision were ~0.086 and ~.007, respectively. The SNP recall and precision were ~0.44 and ~0.81. Do these low values make sense? If not, please let me know how to fix this. Thanks!

[umahsn]

Hi,

These results are very unexpected. Can you please some more information about the sample and the ground truth used? In particular, what is the coverage and average read length of the sample? What is the ground truth for NA24385 that you are using? Are you also using a high confidence interval for your ground truth?

What is the performance if you evaluate using the method shown here [ONT Case Study.md](https://github.com/WGLab/NanoCaller/blob/master/docs/ONT Case Study.md)? You can replace GRCh38 with GRCh37 in the case study since you are using hg19.

[LilyWang07]

Hi,

Thanks for your quick reply. The sample coverage is ~46x and average read length is 8098. We used the following for the ground truth:

VCF: https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/AshkenazimTrio/HG002_NA24385_son/NISTv4.2.1/GRCh37/HG002_GRCh37_1_22_v4.2.1_benchmark.vcf.gz
High confidence region BED: https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/AshkenazimTrio/HG002_NA24385_son/NISTv4.2.1/GRCh37/HG002_GRCh37_1_22_v4.2.1_benchmark_noinconsistent.bed

I will try the case study method. Please let me know if I can provide any other info.

[umahsn]

Thank you. Is the sequencing data publicly available or would you be able to share? I would love to run NanoCaller on it and see if I can reproduce the issue. Also, can you share the hap.py command you used for evaluation, as well as the alignment commands?

[LilyWang07]

Yes, the sequencing reads are from:

wget https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/AshkenazimTrio/HG002_NA24385_son/UCSC_Ultralong_OxfordNanopore_Promethion/GM24385_1.fastq.gz
wget https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/AshkenazimTrio/HG002_NA24385_son/UCSC_Ultralong_OxfordNanopore_Promethion/GM24385_2.fastq.gz
wget https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/AshkenazimTrio/HG002_NA24385_son/UCSC_Ultralong_OxfordNanopore_Promethion/GM24385_3.fastq.gz

and they are aligned with minimap2:
minimap2 -t 30 --MD -Y -L -a -x map-ont $ref $fastq

for hap.py evaluation:
hap.py ${bench_vcf} ${query_vcf} -f ${conf_bed} -r ${reference} -o ${output_dir}/${prefix}

Let me know if I can provide more information. Thanks!

[umahsn]

Thanks, I am running some tests on this dataset and will let you know soon.

[umahsn]

Hi Lily,

I used the sample you mentioned above and here are my steps and results.

I downloaded the following reference genome for GRCh37 and indexed it:

wget http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz -O -| gunzip -c > GRCh37.fa
samtools faidx GRCh37.fa

# SDF file needed for RTG evaluation later
rtg RTG_MEM=4G format -f fasta GRCh37.fa -o GRCh37.sdf

I downloaded the ground truth VCF and high confidence BED files:

wget https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/AshkenazimTrio/HG002_NA24385_son/NISTv4.2.1/GRCh37/HG002_GRCh37_1_22_v4.2.1_benchmark.vcf.gz
wget https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/AshkenazimTrio/HG002_NA24385_son/NISTv4.2.1/GRCh37/HG002_GRCh37_1_22_v4.2.1_benchmark.vcf.gz.tbi
wget https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/AshkenazimTrio/HG002_NA24385_son/NISTv4.2.1/GRCh37/HG002_GRCh37_1_22_v4.2.1_benchmark_noinconsistent.bed

Then, I aligned the reads to GRCh37 using minimap2 (v2.28-r1209):

 minimap2 -t 28 --MD -Y -L -a -x map-ont GRCh37.fa GM24385_1.fastq.gz GM24385_2.fastq.gz GM24385_3.fastq.gz|samtools view -Shu|samtools sort -O BAM -o HG002.GRCh37.bam -@ 4

Then I used NanoCaller v3.6.0:

python NanoCaller-3.6.0/NanoCaller --bam HG002.GRCh37.bam --ref GRCh37.fa --seq ont --snp_model ONT-HG002 --indel_model ONT-HG002 --sample HG002 --output test_hg19_all --cpu 32
tabix  test_hg19_all/variant_calls.vcf.gz

Then, I used RTGtool's vcfeval function for evaluation:

rtg  RTG_MEM=12G vcfeval -b HG002_GRCh37_1_22_v4.2.1_benchmark.vcf.gz -c test_hg19_all/variant_calls.vcf.gz -t GRCh37.sdf -e HG002_GRCh37_1_22_v4.2.1_benchmark_noinconsistent.bed -Z -o evaluation -f 'QUAL'

I used the following command to display the evaluation results:

for TYPE in {snp,non_snp};do echo "$TYPE performance"; printf 'Quality_Cutoff\tTP_baseline\tFP\tTP_call\tFN\tPrecision\tRecall\tF1\n'; tail -n +8 evaluation/${TYPE}_roc.tsv |awk -v max=0 '{if($8>max){want=$0; max=$8}}END{print want}';tail -n +8 evaluation/${TYPE}_roc.tsv |awk -v max=0 '{if($8>=max){want=$0; max=$8}}END{print want}';tail -1 evaluation/${TYPE}_roc.tsv; echo "";done

and here are the results:

SNP Performance

Quality_Cutoff	TP_baseline	FP	TP_call	FN	Precision	Recall	F1
5.968	3188097.56	45164.00	3187329.00	164581.44	0.9860	0.9509	0.9682
4.688	3197340.52	55019.00	3196567.00	155338.48	0.9831	0.9537	0.9682
3.010	3208804.00	77221.00	3208023.00	143875.00	0.9765	0.9571	0.9667

Indel Performance

Quality_Cutoff	TP_baseline	FP	TP_call	FN	Precision	Recall	F1
22.490	105600.90	285709.00	105889.00	416725.10	0.2704	0.2022	0.2314
22.340	107021.99	296594.00	107313.00	415304.01	0.2657	0.2049	0.2314
0.220	163742.00	13957060.00	163707.00	358584.00	0.0116	0.3135	0.0224

I noticed that when I initially used GRCh37 downloaded from (https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000001405.25/), I had similarly low performance. Specifically, I used https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.25_GRCh37.p13/GCF_000001405.25_GRCh37.p13_genomic.fna.gz. I believe the poor performance was due to presence of alt contigs. Please let me know which version of the GRCh37 reference genome you used.

I was unable to get hap.py to run since it uses a very antiquated version of python, and used rtgtool's vcfval function instead.

[LilyWang07]

Hi,

Thanks so much for your help! I can try using the reference you mentioned.

My reference was downloaded from http://support.10xgenomics.com/genome-exome/software/pipelines/latest/advanced/references and is labeled hg19-2.1.0. Please let me know if I can provide any additional info.