Strange mapping on miRNAs from RNA-seq library

Question

[manekineko]

Hi,
I'm observing strange mapping on the miRNA precursors from RNA-seq PE reads library. I have checked the library and almost 100% of the reads are 100nt. So probably bowtie/tophat make some trimming and mapping here (I made a mapping only with bowtie (by Galaxy - option end-to-end fast) the result is picture below. Strangely I observed this only on miRNAs. The mapped fragments are 22-24nt. (PS. This is an general TrueSeq RNA-seq library 100 PE). I'm confused what it can be?

Answer 1

[andrew.j.skelton73]

There seems to be higher coverage to only a fraction of the hairpin sequence, which might imply that what's aligning is a 3p or 5p mature sequence possibly?

[manekineko]

Yes, but the question is how, as there is no small RNAs in the sample library....and no <30nt reads in the file....

[Brian Bushnell]

I assume these are clipped alignments. If you hover over one of these reads it will display the cigar string... I suggest you post that here. But more importantly, I recommend you adapter-trim your reads prior to mapping, so they don't need to be clipped.

[manekineko]

They are adapter-trimmed:

BBDuk version 35.14
maskMiddle was disabled because useShortKmers=true
Initial:
Memory: max=33585m, free=32358m, used=1227m

Added 250901 kmers; time:     0.723 seconds.
Memory: max=33585m, free=30431m, used=3154m

Input is being processed as paired
Started output streams:    0.100 seconds.
Processing time:           407.442 seconds.

Input:                      91156618 reads         9206818418 bases.
KTrimmed:                   717213 reads (0.79%)     39491028 bases (0.43%)
Result:                     91043132 reads (99.88%)     9167327390 bases (99.57%)

Time:               408.319 seconds.
Reads Processed:      91156k     223.25k reads/sec
Bases Processed:       9206m     22.55m bases/sec

[andrew.j.skelton73]

According to that <1% of reads had adapters, doesn't that seem odd?

[Brian Bushnell]

It doesn't, to me... as long as the average insert size is substantially greater than read length, and particularly if the reads are size-selected to a longer insert size, you can get pretty close to zero percent of reads having adapters. To validate this, an insert-size histogram from mapping or merging is quite useful.

[andrew.j.skelton73]

valid point - I agree that a histogram of insert sizes would be useful

[manekineko]

RNA was selected 200-400nt, general TrueSeq protocol.

I used Picard to estimate the insert size mean - ~193, so the distance between the pairs is around zero.

[Brian Bushnell]

It looks to me like these are proper small RNAs, and only a small fraction (<1%) of the reads map to them. Why are you doubting this result? Did you do some chemical process that you expected to remove all of the small RNAs?

[manekineko]

I'm doubting because there is no smallRNA reads after adapter removal according histogram (100% of reads fall into 100nt group):

0    0    0.000%    45578309    100.000%    0    0.000%    4603409209    100.000%
10    0    0.000%    45578309    100.000%    0    0.000%    4603409209    100.000%
20    0    0.000%    45578309    100.000%    0    0.000%    4603409209    100.000%
30    0    0.000%    45578309    100.000%    0    0.000%    4603409209    100.000%
40    0    0.000%    45578309    100.000%    0    0.000%    4603409209    100.000%
50    0    0.000%    45578309    100.000%    0    0.000%    4603409209    100.000%
60    0    0.000%    45578309    100.000%    0    0.000%    4603409209    100.000%
70    0    0.000%    45578309    100.000%    0    0.000%    4603409209    100.000%
80    0    0.000%    45578309    100.000%    0    0.000%    4603409209    100.000%
90    0    0.000%    45578309    100.000%    0    0.000%    4603409209    100.000%
100    45578309    100.000%    45578309    100.000%    4603409209    100.000%    4603409209    100.000%
110    0    0.000%    0    0.000%    0    0.000%    0    0.000%

[Brian Bushnell]

There's something wrong here... you started with 101bp reads. Adapter removal trimmed 0.79% of them and discarded 0.12%. Therefore, you should still have 0.67% of reads remaining that are shorter than 101bp, but the histogram shows zero. Are you sure you ran it on the trimmed reads?

Also, the cigar string of the read mapped is 23M, indicating that there is no clipping and it is 23bp in length.

[manekineko]

Oh by bad I have mistaken the filnames so after re-run of the trimmed reads from R1 pair:

So it seems that these are actual smallRNAs probably contaminations in RNA-seq? Do you think I can use it to see come diff.expressed miRs or the mapped reads are very few, As I have some values from cuffdiff on miR genes from that?

0    0    0.000%    45521566    100.000%    0    0.000%    4584970307    100.000%
10    35053    0.077%    45521566    100.000%    526628    0.011%    4584970307    100.000%
20    43342    0.095%    45486513    99.923%    1037073    0.023%    4584443679    99.989%
30    25747    0.057%    45443171    99.828%    883980    0.019%    4583406606    99.966%
40    18330    0.040%    45417424    99.771%    814247    0.018%    4582522626    99.947%
50    17896    0.039%    45399094    99.731%    977696    0.021%    4581708379    99.929%
60    19892    0.044%    45381198    99.692%    1284650    0.028%    4580730683    99.908%
70    29270    0.064%    45361306    99.648%    2191778    0.048%    4579446033    99.880%
80    83122    0.183%    45332036    99.584%    7113941    0.155%    4577254255    99.832%
90    0    0.000%    45248914    99.401%    0    0.000%    4570140314    99.677%
100    45248914    99.401%    45248914    99.401%    4570140314    99.677%    4570140314    99.677%
110    0    0.000%    0    0.000%    0    0.000%    0    0.000%

[Brian Bushnell]

The differential expression program will have a statistical model in it that determines whether the number of reads is sufficient to be significant. However, if you did some selection process that was supposed to remove small RNAs, then you should ignore them because that will bias the results.