rrwick / Porechop

adapter trimmer for Oxford Nanopore reads
GNU General Public License v3.0
323 stars 124 forks source link

Barcode demultiplexing after albacore #39

Closed olgakarinsky closed 6 years ago

olgakarinsky commented 6 years ago

Dear all, I have encountered something that I can't seem to explain in the results of running Porechop on Albacore demultiplexed data. As I understand, the final results for all barcodes should only contain the reads for which both Porechop and Albacore assigned the same barcode, otherwise the read should be in the 'none.fastq' file.

When looking at the results I receive, I see some cases where the final assignment of porechop is different than the assignment of albacore and the reads are assigned to a certain barcode and not filtered to the 'none' file.

Albacore adds a tag to the read ID in the fastq file with the barcode ID that was chosen. When I look at 'BC01.fastq' file, I expect all headers to contain the value barcode=barcode01 (since disagreements should be filtered out), but instead I see the following histogram of barcodes:

barcode01 | 58075 barcode02 | 13 barcode03 | 9 barcode04 | 7 barcode05 | 3 barcode06 | 1 barcode08 | 2 barcode09 | 9 barcode10 | 5 barcode11 | 1 barcode12 | 2 unclassified | 7409

It is clear that in most cases they indeed agree, but ~7,500 reads were assigned differently.

In order to explore this further, I have ran Porechop using verbosity = 3 and looked at a specific read which was assigned to BC02 by Albacore but to BC01 after running Porechop on it's output and this is the relevant results I see:

ab8c241b-5a89-420a-b723-98b7790d8e44 runid=77c14fc9c92785d9274f3d85723a4061540e1a40 read=1970 ch=34 start_time=2017-12-11T12:44:48Z barcode=barcode02 start: TCATGCTTCGTTCAGTTACGTATTGCTTCAGTTCGTTTACATGAAAGTCGTCTGTGTTTTCGCATTTATCGACAAACGCTTTCGCGTTTCGTGCGCCGCTTCACGCAGGCCTTCATCTTTCTGGCTGATGTACCACATTGCCTGCCGCGA... start alignments: SQK-NSK007, full score=89.285714, partial score=89.285714, read position: 2-27 Rapid, full score=92.0, partial score=92.0, read position: 55-103 1D^2 part 2, full score=81.818182, partial score=81.818182, read position: 5-33 Rapid barcoding 1 (full sequence), full score=76.315789, partial score=76.315789, read position: 2-103 Rapid barcoding 2 (full sequence), full score=85.185185, partial score=85.185185, read position: 2-103 Rapid barcoding 3 (full sequence), full score=75.221239, partial score=75.221239, read position: 2-103 Rapid barcoding 5 (full sequence), full score=78.378378, partial score=78.378378, read position: 2-103 Rapid barcoding 6 (full sequence), full score=79.279279, partial score=79.279279, read position: 2-103 Rapid barcoding 7 (full sequence), full score=79.464286, partial score=79.464286, read position: 2-103 Rapid barcoding 8 (full sequence), full score=77.876106, partial score=77.876106, read position: 2-103 Rapid barcoding 9 (full sequence), full score=78.378378, partial score=78.378378, read position: 2-103 Rapid barcoding 10 (full sequence), full score=78.378378, partial score=78.378378, read position: 2-103 Rapid barcoding 11 (full sequence), full score=75.438596, partial score=75.438596, read position: 2-103 Rapid barcoding 12 (full sequence), full score=79.090909, partial score=79.090909, read position: 2-103 end: ...CTGGACATTAGTACTACTCATTTTACATTAATGCTGGAATGCGTCTGGAGAAACGGGAATATGACGTTACCTCTCCCCAAAGAGAATGGGGAGAAGATTTTCTGCATGTGCATTTTATGCACGGCGATAAAGAATGCTGCTGGGTGGAGT Barcodes: start barcodes: BC01 (75.0%), BC02 (70.0%), BC03 (69.0%), BC04 (65.5%), BC05 (58.6%), BC06 (60.0%), BC07 (55.2%), BC08 (45.8%), BC09 (25.0%), BC10 (8.3%), BC11 (59.3%), BC12 (12.5%) end barcodes: BC01 (8.3%), BC02 (64.3%), BC03 (46.4%), BC04 (55.6%), BC05 (34.6%), BC06 (66.7%), BC07 (4.2%), BC08 (16.7%), BC09 (35.7%), BC10 (12.0%), BC11 (4.2%), BC12 (12.5%) best start barcode: BC01 (75.0%) best end barcode: BC06 (66.7%) final barcode call: BC01

I will add that I see around the same percentage of reads with different assignments in all the barcodes (not only BC01).

Am I missing something? What is the possible explanation for this?

Thank you in advance, Olga Karinsku

rrwick commented 6 years ago

Hi Olga,

This was an interesting one!

Porechop can only determine Albacore's barcodes if it is given a whole Albacore directory for its input (-i). Then for each fastq it finds in the directory, it looks for the barcode name in the path. For example, /barcode01/ in the path for this fastq: workspace/pass/barcode01/fastq_runid_4409128d0616a0c2bb18e7405515d7fa445a812b_0.fastq Porechop does not use the info in the read header (though that's not a bad idea - I'll consider adding it). So I suspect you ran Porechop only on the one file, not on the whole Albacore directory.

If my assumption is correct, Porechop never knew what the Albacore barcode call was. If you run Porechop like this porechop -i albacore_dir -b output_dir, then I think it will put the read in the 'none' bin, as you expected.

What I found interesting about this case is this: Albacore got the barcode correct (barcode02) and Porechop got it wrong (barcode01).

Here's the alignment of barcode01 to your read's start:

TCATGCTTCGTTCAGTTACGTATTGCTTCAGTTCGTTTACATGAAAGTCGTCTGTGT-TTT-CGCATTTATCGAC
                                        | |||||| ||| |||| |||  |           
----------------------------------------AAGAAAGTTGTCGGTGTCTTTGTG-----------

And here's the alignment of barcode02:

TCATGCTTCGTTCAGTTACGTATTGCTTCAG-TT-CGTTTACATGAAAGTCGTCTGTGTTTTCGCATTTATCGAC
                           || | || ||||    ||  ||||||||||                  
---------------------------TC-GATTCCGTT----TG-TAGTCGTCTGT------------------

Barcode01 indeed looks like a better match (which is why Porechop called it as such), but I can see from the flanking sequences that the barcode01 alignment is in the wrong place and the barcode02 alignment is in the right place! In a rapid kit, the barcode is preceded by ...CGTATTGCT and followed by GTTTTCGCA..., which matches the barcode02 alignment.

It's an interesting case, and it makes me wonder if and how Porechop should be considering those flanking sequences. But I'm not too worried about it - I suspect it's an unusual case and Porechop would 'none' bin the read because of the Albacore disagreement.

I hope that clears it up. I'm going to close this issue now, but if something's still not right or you still have questions, please let me know!

Ryan