direct RNA adapter sequences are available

[phiweger]

hi Ryan,

here

I would need to incorporate them in the adapter.py script, correct?

thank you very much,
Adrian

[rrwick]

Hi Adrian,

Thanks for the link! Yes, adding those adapters to the adapters.py file should be all that's required. I went ahead and did that on Porechop's development branch, so give it a try if you'd like!

To use Porechop's development branch, follow the install from source instructions, but use this clone command: git clone -b development https://github.com/rrwick/Porechop.git

I could not, however test these new adapters. I don't have any direct RNA read sets myself, and I looked for publicly available ones on SRA/ENA but didn't find much. There were some read sets here, but they don't seem to have any sequences matching the direct RNA adapters. Maybe they were trimmed before uploading?

Do you have a direct RNA read set that you could share with me? I would promise to keep the data private and only use it for Porechop development. Even a smaller piece of a read set (100 Mbp or so) would be helpful. Or do you know of anywhere else I could get my hands on one? I've asked Oxford Nanopore if they have internal datasets to share, but I had no luck on that front.

Ryan

[rrwick]

Hi Adrian, I just saw in issue #27 that you said 'I could send you some data'. If by that you mean a read set, then yes please! 😄

Dropbox works for me, I have an account tied to my email: [email protected]. Or else I'm open to other suggestions for moving around large files.

Ryan

[phiweger]

There should be something in your inbox now :)

[rrwick]

Thanks for the dataset, and thanks @nickschurch for that link. I've taken a look at both and come to a tentative conclusion.

Very short version:
I don't think it's going to work.

Short version:
I suspect that the adapter which comes after the RNA is actually DNA, which throws off the basecaller. So while reads do contain sequences which should be trimmed, they are too inconsistent for Porechop to work with.

Full version:
As far as I could tell from these diagrams and other ONT materials, the Reverse transcription adapter (RTA) should be the first thing encountered after the poly-A tail. Interesting side note: direct RNA sequencing seems to go through the pore 'backwards' (3' to 5') so the adapter is at the start of the signal but at the end of the basecalled read. To find the adapter sequence, I looked in the direct RNA datasets for poly-A tails near the end of reads and examined sequences that followed.

What I found was quite inconsistent. It often looked something like UCCCAUCCCAUACACUCCCACAU, but with a lot of variation. Gs were very uncommon and they did not look at all like the RTA sequence ONT provides, which has a lot of Gs. I found a clue on this page. In the section describing how to make custom RTA sequences, it says 'Both oligo A and oligo B are DNA.'

So here's what I think is happening: the sequenced strand is mostly RNA, but the adapter part is actually DNA. Since the basecaller was trained on RNA, it doesn't quite know what to do with the signal in the DNA region, and what results is a weird G-less mess. The signal image on this page shows the poly-A tail in green, and the red part, which has a distinctly different looking signal from the rest, is what I think is the DNA adapter.

Ultimately this means Porechop isn't the right tool for the job. Porechop trims adapters based on sequence identity, and if the basecaller isn't making consistent sequence, it won't have anything good to work with. This sort of trimming should really be done at the signal level, and is therefore probably a job for the basecallers. Hopefully a future version of Albacore/Guppy can do this.

Phew, I think that's everything! 😅 In conclusion, I'm not going to put the direct RNA adapters in Porechop because I don't think it will work. I'll raise a support ticket with ONT, asking them if I got anything wrong or if they have anything else to add. I'll close this issue now because I don't think there's anything more to be done with Porechop, but feel free to keep posting if you still have questions.

Ryan

[wdecoster]

I wonder if those DNA adapter sequences get a (very) low quality score from the RNA basecaller, providing an alternative method to detect/clip them?

[phiweger]

I really appreciate your looking into this. I trimmed the 3' end of some reads and did come to a similar conclusion, i.e. there was pretty much no identity to speak of to the ONT direct RNA adapters. I find your RNA-DNA-basecaller hypothesis convincing. I will probably hard trim the poly(A) + 3' sequence for now. Thanks a lot.

[rrwick]

@wdecoster, you are right! I grabbed reads with a long poly-A (25 bp or more), and trimmed them such that the poly-A ended at position 100. So everything after 100 is what should presumably be trimmed. I then put them through FastQC and got this:

The pre-poly-A quality averages around 9-10. It then spikes for the poly-A itself, and then falls to around 5-6 for the remaining sequence.

Whether this could be used to reliably trim reads, even when there is little poly-A, I'm not sure. I'd be interested if others had luck here, but I think this falls outside of Porechop's domain, so I don't plan on attempting it myself.

Ryan

[rrwick]

Here is the file I analysed if you're interested: polyA_ends.fastq.txt

The reads are from this dataset. And sorry about the '.txt' - GitHub won't attach the file without it.

[rrwick]

Final update:
I heard back from ONT, and my suspicions were correct: the direct RNA adapter is DNA so that's why it's not being basecalled to a consistent sequence. Also, they are planning to include adapter trimming in the basecallers at some point in the future, but no ETA yet for that feature.