ShapeMapper QC Fail

[camdenbair]

I am having an issue similar to "QC reliability, #42". I am processing a ShapeMap data set where I used the amplicon workflow for a ~400 nt RNA. I pass the first three quality control checks but fail on the "Number highly reactive check".
I have received this error with with 3 different sources of 1M7 (commercial and two different batches which were synthesized by a collaborator). We have also replaced several of the reagents required for the workflow (DMSO, MnCl2, Superscript II, nTPs, etc).
Any help narrowing down the problem would be very helpful.

Thank you in advance!

ShapeMapper Histograms.pdf

ShapeMapper Profiles.pdf

[Psirving]

I think many of the points I made #42 are valid here as well. Based on your profile, I think the most likely cause is "extremely highly structured RNA", and that your data mostly look okay, the experiment worked. The main thing I notice is that the quality near the 3' end is low, correlated with the low read depths. I would mask the region around 355 to the end. For a region of this size, often the small RNA workflow is better, because you can sequence the entire region with a 500/600 cycle kit. This avoids the issue of losing sequencing depth at the ends of the RNA. Also, you should consider one of the newer, more reactive SHAPE reagents, but this doesn't appear to causing an issue in your case.

[camdenbair]

@Psirving Thank you very much for your insight.
The current data set does reflect a highly structured RNA. I considered using the small RNA workflow, but that option was not available at the time of sequencing. I'm currently determining if that is an option now, considering the drop in coverage; thank you for the suggestion.

During the same sequencing run, I also processed a 16s rRNA sample (in vitro using 1M7) using the randomer approach and we also received the same error. This data set looks very different from the ShapeMapper2 16s example data set.
I thought we might have a systemic problem since we received the same error in multiple RNAs.

From my limited experience, it appears the 1M7 is not as reactive as it should be (we're using it at 10mM final).
Any additional thoughts would be greatly appreciated!

16S_rRNA_singlestranded_histograms.pdf

16S_rRNA_singlestranded_profiles.pdf

[Psirving]

It looks like your "denatured" control much more reactive than the "modified" sample, suggesting your 1M7 is at least somewhat reactive and that structure is playing a role. I also see individual positions with decent reactivity in the modified sample. It might be worth analyzing without the denatured control, which sometimes causes more trouble than its worth. My hunch is that your 1M7 is okay, maybe not 100%.

Can you get your hands on that other 16S data set? If so, I would want to see linear regressions comparing the raw mutation rates of each of the samples:

• modified vs modified
• untreated vs untreated
• denatured vs denatured
• modified-untreated vs modified-untreated

Another thing that I suspect can happen is that structured RNAs can hold on to cellular Mg++ (depending on purification method) and affect the RT reaction (which needs Mn++ for MaP to work). It could be that the denatured sample has "let go" of the cellular Mg++ and the MaP read out is affected. This is just speculation.

[camdenbair]

I did a Pearson correlation (using the Modified_rate, Untreated_rate, and Denatured_rate values) between the two data sets. The sample data set was downloaded form https://weekslab.org/software/ under the "SHAPE-MaP" heading, "Recommended sample data in a ready to use format".

The R2 values for each comparison are as follows:
-modified vs modified: 0.4553
-untreated vs untreated: 0.4508
-denatured vs denatured: 0.4298
-modified-untreated vs modified-untreated: 0.0001973

I also reran my dataset removing the denatured control which resulted in profiles that were largely similar the profiles with the denatured control. The attached image show the SHAPE reactivity profiles of my dataset ran with and without denatured control as well as the sample "Weeks" dataset.

[Psirving]

I still think your Zika result looks like decent 1M7 data (except after ~350). There's obviously something going on with the 16S data. I'm mostly convinced 1M7 is not the issue, otherwise the denatured mutation rate would not be elevated. Are there any significant differences between the probing of these two RNAs, or between our probing and your probing of 16S?

[camdenbair]

I am also mostly convinced that the ZIKV results is accurate (compared with previously published SHAPE profiles). I haven't, however, been able to get my hands on raw data sets for similar RNAs to know if failing the 4th QC ("Number highly reactive check") is normal or not.
As far as the 16S, I am not sure if there were differences in treating our 16S vs what was done with the Weeks 16S data set.
On our end, after purification of the 30S they're heated to 42c for 5 minutes, then treated with 1M7 following the Smola 2015 Nature Protocol. I've never been sure how how the example data set RNA was purified and prepared prior to 1M7 treatment. Certainly, differences in preparation/treatment could affect SHAPE reactivity.
With all of that being said, we're going to try a different shape reagent (2A3) as you suggested in Issue #42 to determine if we can get a stronger signal.

[camdenbair]

Thank you @Psirving for your help with this issue.
After some consideration, we repeated the SHAPE-MaP experiment with the more reactive reagent 2A3 (Marinus et al., NAR 2021) and successfully passed the quality control checkpoints. After mapping the SHAPE reactivities to the predicted secondary structure of my RNA, I found the 1M7 and 2A3 experiments were very similar (although the 2A3 experiment had higher reactivities overall).

@Psirving I have taken some time to review your initial comments and suggestions including using the small RNA workflow:
"For a region of this size, often the small RNA workflow is better, because you can sequence the entire region with a 500/600 cycle kit. This avoids the issue of losing sequencing depth at the ends of the RNA."

Just to confirm, are you suggesting to use a 500 cycle kit (paired-end 250 bp reads) to cover a 400-500 nucleotide RNA? I want to confirm that ShapeMapper2 does not expect the R1 and R2 to have 100% overlap. In other words, using a 500 cycle kit (PE-250) for a 400 nucleotide RNA, the R1 and R2 will have some region of homology/overlap in the center of the RNA, but the 5' ends of the R1 read will not overlap with the R2 read, and visa versa.

Thank you!

[Psirving]

Yes, ShapeMapper2 will correctly handle paired reads, no matter the degree of overlap. Overlapping reads are combined pre-alignment using BBMerge. Non-overlapping reads are combined after alignment in the mutation counting step.