Closed ymcki closed 4 months ago
ymcki
commented
1 year ago [2023-08-11 11:54:45.957] [info] > Simplex reads basecalled: 1305230 [2023-08-11 11:54:45.957] [info] > Simplex reads filtered: 532 [2023-08-11 11:54:45.957] [info] > Duplex reads basecalled: 5242 [2023-08-11 11:54:45.957] [info] > Duplex reads filtered: 4275 [2023-08-11 11:54:45.957] [info] > Duplex rate: 0.078480266% [2023-08-11 11:54:45.957] [info] > Basecalled @ Bases/s: 4.441025e+06
My numbers are pretty bad :(
tijyojwad
commented
1 year ago Hi @ymcki - This yield looks extremely low. I've seen this kind from non duplex flow cell runs where we hardly expect duplex pairs to be found (and those which are found are likely false positives).
Also, is this short read data by any chance? dorado duplex isn't setup to handle short reads right now
ymcki
commented
1 year ago How do I know it is short read data?
This is the output from pod5 inspect. So it is 4kHz data with RAD114 kit. Flowcell is PRO114M. Is this the duplex flowcell?
pod5 inspect debug ./pod5/split_by_channel/channel-1935.pod5
Contains 1 reads, in 1 batches: [1]
Reads span from sample 22594207 to 22606826
12619 samples, 10429 bytes: 41.3 % signal compression ratio
Run info 0:
acquisition_id: cfd7091fd46b02e226c8b66fb4fcc31dbaa11d8b
acquisition_start_time: 2023-02-08 20:52:38.378000+00:00
adc_max: 2047
adc_min: 0
context_tags: {'barcoding_enabled': '0', 'experiment_duration_set': '5760', 'experiment_type': 'genomic_dna', 'local_basecalling': '0', 'package': 'bream4', 'package_version': '7.3.5', 'sample_frequency': '4000', 'selected_speed_bases_per_second': '400', 'sequencing_kit': 'sqk-rad114'}
experiment_name:
flow_cell_id: PAK61167
flow_cell_product_code: FLO-PRO114M
protocol_name: krill_scripts/sequencing_PRO114_DNA_e8_2_400K:FLO-PRO114M:SQK-RAD114:400
protocol_run_id: b141cebc-a924-4f36-a3a3-8a95c24fd762
protocol_start_time: 2023-02-08 20:51:18.261000+00:00
sample_id: ABC
sample_rate: 4000
sequencing_kit: sqk-rad114
sequencer_position: 3D
sequencer_position_type: promethion
software: python-pod5-converter
system_name: PC48B200
system_type: PRO-PRC048
tracking_id: {'asic_id': '0004A30B001D5E6F', 'asic_id_eeprom': '0004A30B001D5E6F', 'asic_temp': '34.006222', 'asic_version': 'Unknown', 'auto_update': '0', 'auto_update_source': 'https://mirror.oxfordnanoportal.com/software/MinKNOW/', 'bream_is_standard': '0', 'configuration_version': '5.3.8', 'device_id': '3D', 'device_type': 'promethion', 'distribution_status': 'stable', 'distribution_version': '22.10.7', 'exp_script_name': 'krill_scripts/sequencing_PRO114_DNA_e8_2_400K:FLO-PRO114M:SQK-RAD114:400', 'exp_script_purpose': 'sequencing_run', 'exp_start_time': '2023-02-08T12:52:38.378283-08:00', 'flow_cell_id': 'PAK61167', 'flow_cell_product_code': 'FLO-PRO114M', 'guppy_version': '6.3.9+799d5bd', 'heatsink_temp': '34.179996', 'host_product_code': 'PRO-PRC048', 'host_product_serial_number': 'PC48B200', 'hostname': 'PC48B200', 'hublett_board_id': '0133c0d5a87e0e5c', 'hublett_firmware_version': '2.1.10', 'installation_type': 'nc', 'ip_address': '', 'local_firmware_file': '1', 'mac_address': '', 'operating_system': 'ubuntu 20.04', 'protocol_group_id': 'ABC', 'protocol_run_id': 'b141cebc-a924-4f36-a3a3-8a95c24fd762', 'protocol_start_time': '2023-02-08T12:51:18.261306-08:00', 'protocols_version': '7.3.5', 'run_id': 'cfd7091fd46b02e226c8b66fb4fcc31dbaa11d8b', 'sample_id': 'ABC', 'satellite_board_id': '0139c958aa701b78', 'satellite_firmware_version': '2.1.9', 'sequencer_product_code': 'PRO-SEQ048', 'sequencer_serial_number': '', 'usb_config': 'fx3_0.0.0#fpga_0.0.0#unknown#unknown', 'version': '5.3.1'}I noticed that even the source of this data could only basecall 104.8MB fq.gz for duplex but they were getting 25.1GB fq.gz for simplex using . So actually the number of duplex reads are not much at all for this data. Probably that means R10.4.1 ultra long read duplex were not ready for prime time back in Feb 2023?
Anyway, I noticed that the fastqs generated by the source is having much higher qualities for duplex reads that are expected for duplex reads. So I think I still need to figure out how they did that.
I will try older version of dorado (0.1.1 and/or 0.2.1) and see if I can replicate their base qualities.
ymcki
commented
1 year ago I used 0.2.1 and is now able to get similar number of duplex reads (actually a bit more than the source as they used 0.1.1) for orig and splitduplex and also similar base quality. The input I use is the split_by_channel pod5 I fed to 0.3.4.
So probably something was wrong with 0.3.4? Or I missed some parameters? I noticed 0.3.4 has quite many <500bp duplex reads while 0.2.1 has none but 0.2.1 has quite many <5bp duplex reads that 0.3.4 doesn't have.
Here is a list of commands for me to generate the duplex_orig and duplex_splitduplex with dorado-0.2.1 ./dorado basecaller dna_r10.4.1_e8.2_400bps_sup@v4.1.0 ./pod5/split_by_channel/ --emit-moves > unmapped_reads_with_moves.bam duplex_tools pair --threads 96 --output_dir pairs_from_bam unmapped_reads_with_moves.bam ./dorado duplex -t 96 dna_r10.4.1_e8.2_400bps_sup@v4.1.0 ./pod5/split_by_channel/ --pairs pairs_from_bam/pair_ids_filtered.txt > duplex_orig.bam duplex_tools split_pairs --threads 96 unmapped_reads_with_moves.bam ./pod5/split_by_channel/ pod5s_splitduplex/ ./dorado duplex -t 96 dna_r10.4.1_e8.2_400bps_sup@v4.1.0 pod5s_splitduplex/ --pairs split_duplex_pair_ids.txt > duplex_splitduplex.bam
ymcki
commented
1 year ago I find that unmapped_reads_with_moves.bam contains reads from duplex_splitduplex.bam and duplex_orig.bam.
If I want to get a bam for simplex only, I should remove those reads, correct?
tijyojwad
commented
1 year ago Hi @ymcki - 0.2.1 is quite an old version and we've made many improvements since then. Also 0.2.1 didn't have end to end duplex, so it needed the duplex_tools toolchain to generate pairs files and call duplex. That flow is now deprecated. For proper duplex runs though, I would expect dorado 0.3.4 to perform similarly (if not slightly better).
I checked with our team internally and the RAD kit will not give good duplex yield. Maybe single digit % yield at best.
It sounds like you're attempting to reproduce a previous result. Is the dataset public? Could you share it with us?
ymcki
commented
1 year ago Thanks for letting me know the RAD114 duplex only have single digit % duplex yield.
Out of curiosity, I also went 0.3.4 with duplex_tools workflow. It is worse than 0.2.1 in every measure.
Sample 0.2.1_duplex_origg 0.2.1_duplex_splitduplex 0.2.1_duplex 0.3.4_duplex_origg 0.3.4_duplex_splitduplex 0.3.4_duplex throughput 131,755,956 81,728,072 213,484,028 4,313,238 3,183,788 7,497,026 reads# 2,807 1,716 4,523 1,601 1,199 2,800 Avg Read Length 46,938.35 47,627.08 47,199.65 2,694.09 2,655.37 2,677.51 Longest Read 307,535 360,343 360,343 14,246 29,869 29,869 N50 72,947 90,820 79,017 3,046 4,656 3,339 L50 595 287 874 444 216 629 >=50K bases 92,511,776 61,319,963 153,831,739 0 0 0 >=100K bases 40,961,253 36,750,510 77,711,763 0 0 0 >=900K 0.0000% 0.0000% 0.0000% 0.0000% 0.0000% 0.0000% >=500K 0.0000% 0.0000% 0.0000% 0.0000% 0.0000% 0.0000% >=400K 0.0000% 0.0000% 0.0000% 0.0000% 0.0000% 0.0000% >=300K 0.0356% 0.2331% 0.1105% 0.0000% 0.0000% 0.0000% >=200K 0.6056% 1.8648% 1.0834% 0.0000% 0.0000% 0.0000% >=100K 10.7232% 14.1608% 12.0274% 0.0000% 0.0000% 0.0000% Mean Reported Base Quality 51.73867375 49.74932494 50.97709152 12.17244121 12.17771723 12.1746818
I am not sure if my data is public or not as I found it by browsing a public ftp site. I didn't see it being publicized at all. How about you give me an ftp site to upload the split_by_channel pod5 directory for you to debug?
ymcki
commented
1 year ago I also tried LSK114 duplex data from the same source but it still doesn't work for 0.3.4 but it worked for 0.2.1.
Anyway, I aligned the duplex reads and simplex reads from the RAD114 duplex data in previous reply to HG002 genome and got these stats.
Sample duplex_orig duplex_splitduplex duplex simplex Aligned (qs>=10) 2,733 1,580 4,313 1,011,180 Unmapped (qs>=10) 58 86 144 86,321 Aligned (qs<10) 10 24 34 137,736 Unmapped (qs<10) 6 26 32 62,538 Bases Aligned (qs>=10) 131,062,596 80,896,726 211,959,322 27,075,362,319 %Bases Aligned (qs>=10) 99.4738% 98.9828% 99.2858% 85.0999% %mismatch (qs>=10) 1.0857% 1.3421% 1.1796% 2.2476% %mismatch (qs<10) 7.1154% 10.4193% 9.4476% 16.0698%
We can see that the improvement in mismatch rate is from 2.25% to 1.09%, ie Q16.5 to Q19.6. I was expecting at a jump of 10 q-values as indicated by the jump of the reported quality. From my analysis of another 8 runs of non-duplex RAD114 run, the average mismatch% is 1.77%. So unless you want absolute best quality, there doesn't seem to make sense to run duplex. The conclusion seems to be while duplex can improve quality but the jump is not as dramatic as claimed. Of course, this can be due to the inherent defects in the RAD114 duplex kit.
Is it expected that the duplex_splitduplex has worse quality than duplex_orig as it is a chimeric read consists of two reads from different strands? https://github.com/nanoporetech/duplex-tools/blob/master/fillet.md
I will try the same exercise on the LSK114 duplex data and see what I observe.
tijyojwad
commented
1 year ago Hi @ymcki
I found it by browsing a public ftp site
If it's a public ftp site, then we should be a bee to access sit too. Could you post that? I think we can take a look and let you know if it's usable as a duplex data source or not.
We do typically see the same increase in alignment accuracy (and that is what is generally reported when we talk about duplex improvements, not just the mean base q-scores). I would be careful drawing a lot of conclusions from this dataset since we don't know yet if this is a proper duplex run. We can check once you share the link to the data.
ymcki
commented
1 year ago Thanks for your reply. But I am still downloading the data. I think I should at least finish downloading them before I show you the link. How about you find R10.4.1 duplex data that were generated internally before Feb 08, 2023 and see if you can observe the same problem with 0.3.x.
In the last London calling, Clive mentioned the R10.4HD flowcell. Is it a flowcell that came after Feb 08, 2023 such that the flowcell now is a different beast than the one back then?
tijyojwad
commented
1 year ago For questions about flow cells, you can reach out to support@nanoporetech.com
Once you are able to share the data, we'd be happy to take a look. You can also reach out to me over email at joyjit dot daw at nanoporetech dot com to share the data if you prefer.
ymcki
commented
1 year ago Thanks for your reply. It probably takes me a few more days to download the data. Then I will email you the link.
I am also done with analysis of the LSK duplex data dated Jan 2023 using 0.2.1. The yield is about 26% vs 0.67% for the ULK duplex data. Error profile is similar. I remember in one of Clive's R10.4HD slide, he was talking about 100Gbp duplex and 50Gbp simplex from one flowcell. I suppose that would mean 50% yield (100Gbp duplex => 50Gbp basecalled). So 26% is a bit disappointing but of course compare to the UL one, it is more useful. Perhaps R10.4HD is a different beast?
QC stats for unmapped bam
Sample R10.4.1_UL_Duplex_1_orig R10.4.1_UL_Duplex_1_splitduplex R10.4.1_UL_Duplex_1_simplex R10.4.1_Duplex_1_orig R10.4.1_Duplex_1_splitduplex R10.4.1_Duplex_1_simplex throughput 131,755,956 81,728,072 31,815,975,055 11,105,053,330 7,707,946,491 53,141,467,609 reads# 2,807 1,716 1,297,775 433,122 322,472 3,167,812 Avg Read Length 46,938.35 47,627.08 24,515.79 25,639.55 23,902.68 16,775.45 Longest Read 307,535 360,343 995,521 166,765 159,803 643,857 N50 72,947 90,820 97,894 38,115 37,572 46,256 L50 595 287 96,964 111,131 78,030 378,724 >=50K bases 92,511,776 61,319,963 23,160,312,438 2,634,956,045 1,763,844,587 23,885,540,683 >=100K bases 40,961,253 36,750,510 15,620,198,375 26,903,891 15,527,597 5,109,480,386 >=900K 0.000000% 0.000000% 0.000077% 0.000000% 0.000000% 0.000000% >=500K 0.000000% 0.000000% 0.014717% 0.000000% 0.000000% 0.000095% >=400K 0.000000% 0.000000% 0.074281% 0.000000% 0.000000% 0.000221% >=300K 0.035625% 0.233100% 0.342047% 0.000000% 0.000000% 0.000631% >=200K 0.605629% 1.864802% 1.601664% 0.000000% 0.000000% 0.031820% >=100K 10.723192% 14.160839% 7.247327% 0.054950% 0.043105% 1.311505% Mean Reported Base Quality 51.73867375 49.74932494 34.93595437 49.79552859 50.22268231 36.20575354
QC stats after aligning to HG002v0.7
Sample R10.4.1_UL_Duplex_1_orig R10.4.1_UL_Duplex_1_splitduplex R10.4.1_UL_Duplex_1_simplex R10.4.1_Duplex_1_orig R10.4.1_Duplex_1_splitduplex R10.4.1_Duplex_1_simplex Aligned (qs>=10) 2,733 1,580 1,011,180 387,932 285,742 2,478,713 Unmapped (qs>=10) 58 86 86,321 24,372 21,967 295,753 Aligned (qs<10) 10 24 137,736 7,151 2,792 265,472 Unmapped (qs<10) 6 26 62,538 13,667 11,971 127,874 Bases Aligned (qs>=10) 131,062,596 80,896,726 27,075,362,319 10,874,613,257 7,615,017,322 46,320,768,485 %Bases Aligned (qs>=10) 99.4738% 98.9828% 85.0999% 97.9249% 98.7944% 87.1650% %mismatch (qs>=10) 1.0857% 1.3421% 2.2476% 1.1852% 1.2818% 2.1841% %mismatch (qs<10) 7.1154% 10.4193% 16.0698% 9.9108% 9.5205% 15.5533%
vellamike
commented
4 months ago Closing due to inactivity. Please request a reopen if you still require help.
I am getting average Q5 reads from my 4kHz dorado-0.3.4-rc2. I had similar problem with 0.3.2
The command line
When I grep the duplex reads at the bam file, here are the first three lines. As you can see quite many bases with quality " which means 1 and will be flagged as fail.
My data was in fast5 format. Then I converted to pod5 and then split them by channel as recommended by README.
The output from my simple hdf5 script on one of the fast5.
What went wrong?