SeqMonk: Visualisation and analysis for large mapped data sets

[simonandrews]

Originally Posted by Neuromancer

Hey all,

how does Seqmonk count paired-end reads? Is each pair only counted once (or once per gene?) or is each read counted individually? In any case: is there any way how to switch between these two modes?

How paired end reads show up in the program is determined when you first import the data. If you choose to import your data as paired end then the reads will appear as a single location which spans the region inferred to be covered by the pair. The direction of this single region will be taken from the first read in the pair.

If you import as single end then both ends of the pair will be shown as separate reads but there will be no connection between them in the internal data model so you can't switch between the two views within the same set of imported data.

One of the main trade offs which seqmonk makes in order to allow it to handle large datasets quickly is that it doesn't maintain links between alignment segments, either for paired reads, or for splice segments in spliced reads. For our internal quantitation of spliced data we use the relative length of each aligned segment to infer how many reads we should count when we're summing up the contribution of different spliced segments.

In your case as long as you're correcting for total read counts then it shouldn't matter too much that you have a mix of single and paired end data. In terms of counts the paired data will be somewhat similar to simply doubling a single end sample, and the global correction will normalise this away. If you want to more explicity correct for this then you could apply a manual correction to halve the counts for your paired end data (this is one of the quantitation options).

[Neuromancer]

Hey Simon,

after upgrading to the new version Seqmonk doesn't recognize the the chromosome names of the bowtie2 standard sam-output file anymore (still using the same bowtie2-provided index for mm10) Error message says, it cant make a name of chrM... However, in a sorted bam file, it reads all the lines until it reaches chrM and then quits with the same error, leaving no reads to look at...
Is that because of the new version or did I do something crazy here?
If not, can I somehow exclude importing unreadable chr names (as I'm not interested in chM anyways) and still importing all others?

Thanks a lot for your help so far and (hopefully) in the future

r

edit:
rolling back to 0.22 solves the problem. the same sam/bam/sortedbam are read without complaints! Well, complaints are coming about chrM in the end, but the program just reports that and doesn't break in between...

[simonandrews]

This is a bug in the latest version caused by the way we handle chromosome name matches internally. I've just put a fix into the development version and writing the release notes for a new release right now.

There should be an update out in a couple of hours which will fix this. Sorry for any trouble this has caused.

[simonandrews]

There is a new release of seqmonk (v0.23.1) which should fix the import problems found in v0.23.0. It's up on the project page now and you should be able to get it.

Please let me know if any problems persist with this new version.

[Neuromancer]

Originally Posted by simonandrews

There is a new release of seqmonk (v0.23.1) which should fix the import problems found in v0.23.0. It's up on the project page now and you should be able to get it.

Please let me know if any problems persist with this new version.

Yepp - up and running!
Thanks for that REALLY fast update!
It know works as the version beofre, chrM could, however still not be extracted, but the workaround for that can be found here...

[mathew]

I ahve a question about using scatter plot in seqmonk. When you scatter plot two read counts/ expression there is cloud of red dots in in the center what that mean? Alternatively Can someone point me what various colors in such scatter plot mean

[simonandrews]

The colours in the scatterplot represent the density of points which are overlaid at that point in the plot. There are normally way too many points to be able to show each one, so we use the colours to show where the plot high densities of points are being found. The colour scheme is the standard cold to hot colours used elsewhere in the program.

[mathew]

Thanks Simon for quick reply So from the attached image - red area represents a perfect correlation (or close to 1) and as we move away from line it decrease.
Is it reasonable or I misunderstood something.

Thanks

[simonandrews]

The red area simply shows regions of the plot where many probes are packed on top of each other. There are never enough pixels on the screen to show every probe independently so the colours simply relate to how many probes are overlaid at a particular position. In your plot it shows that the largest number of probes are in a region of the plot showing very little change between the two conditions you've plotted out.

[mjp]

Q1: What is the shortest way to obtain a list of probes that have 'FWD only', 'RVR only' and 'NONE reads' covering them across multiple data stores independently. The idea is to avoid going through the same steps for each datasets.

Currently I'm following this workflow, however I can't get it to work on multiple datastores:
1. Defining my probes.
2. Quantintating FWD reads only
3. Filtering on value -> Individual probes -> value above 1, for exactly 1 of the one selected 1
4. Quantitating RVR reads only
5. Filtering on value -> Individual probes -> value above 1, for exactly 1 of the one selected 1
6. Filtering by combining existing list:
6.1. 'RVR value above 1' BUTNOT 'FWD value above 1'
6.2. 'FWD value above 1' BUTNOT 'RVR value above 1'

Q2: By following this sequence of steps, will 6.2 produce probes with some value for FWD only and not RVR (given that last quantitation was made for RVR reads only?

I believe that should work for single data store given the answer to the Q2 is yes. When I tried (in step 3) to go through multiple datasets I'm a little confused which option to choose.

Would you be able to give me a hint here?

Thanks in advance.

[simonandrews]

I'm not exactly clear what you want to know - do you want a quick way to determine if a given probe is forward or reverse only in all of a set of stores, or are you looking for a quick way to make separate lists for several stores where you have them?

To do the analysis across several stores you could basically repeat the process you outlined but selecting all of the stores, and making your values filters use 'at least 1' rather than 'exactly 1' to pull out probes which had a read in that direction in any of your stores. You could also put all of your data into a single data group and then treat it as a single dataset.

You might also want to try using a difference quantitation, where you could do forward reads as percentage of all reads, and then filter for either 100% forward or 0% forward which might be easier than going through lots of values filters. The only hiccup with this would be that empty features would also show up with 0% so you'd need to have done a read count quantitaiton first and created lists for each of your datasets of probes containing no reads and then use the combine filter to subtract these from the reverse only set (0%) to get the true reverse only count.

[simonandrews]

I have acted on one of my New Year's resolutions and have finally got round to producing some more tutorial videos showing the use of SeqMonk to analyse a number of different datasets covering RNA-Seq (both simple and complex experimental setups), ChIP-Seq, BS-Seq and Hi-C data.

All of the videos can be found at our YouTube channel at:

https://www.youtube.com/user/babrahambioinf

If you have any suggestions for other tutorials which might be useful then please let me know and I'll have a go at putting them together.

[honey]

Hi Simon,

Great job! One suggestion if you can add demo to plot density graphs across TSS or more precisely, how one can plot density graph across TSS as discussed in http://seqanswers.com/forums/showthr...ht=tss+density
will be great. I know close to this can be plotted in Seqmonk, however if we can select specific probes and plot the graphics as suggested will be very helpful.

[simonandrews]

Originally Posted by honey

Hi Simon,

Great job! One suggestion if you can add demo to plot density graphs across TSS or more precisely, how one can plot density graph across TSS as discussed in http://seqanswers.com/forums/showthr...ht=tss+density
will be great. I know close to this can be plotted in Seqmonk, however if we can select specific probes and plot the graphics as suggested will be very helpful.

The ChIP-Seq tutorial constructs one of these types of graph after doing the peak detection part (TSS probes is 6:20 and the plot is about 8:20) so hopefully that covers much of what you wanted. The plot will display whichever probes you have selected through the filters.

The only thing which looks different in the post you linked is that in seqmonk the probes are ordered in the plot by the number of reads covering each probe, whereas the plot you linked to looks like they did something else to order the probes to get some of the patterns you saw. In some older versions we used to cluster the probes but this often produced a messier result than you'd hope. I'm happy to hear suggestions for other ways we could order these plots if there's anything we could do better.

[mjp]

Originally Posted by simonandrews

do you want a quick way to determine if a given probe is forward or reverse only in all of a set of stores, or are you looking for a quick way to make separate lists for several stores where you have them?

Sorry for this small delay and for not being specific enough. I see that the title of my post was wrong.
I wanted to have a list of probes that have just strand specific reads covering them. So the list would contain probes that have only fwd, only rvrs, none, and both type of reads covering them. Ideally I would like to have such a list for each of my stores independently. Is it possible to create that in single SeqMonk pipeline or does it have to be repeated for each store separately?

Originally Posted by simonandrews

To do the analysis across several stores you could basically repeat the process you outlined but selecting all of the stores, and making your values filters use 'at least 1' rather than 'exactly 1' to pull out probes which had a read in that direction in any of your stores.

I'm not quite sure if that does what I want - this is where the confusion starts. I created 2 sample stores for yeast chr1 and run the workflow that I outlined previously but this time selecting 'at least 1' for all selected stores (please see the attached image). Probes generated using running window, size 1 step 1.
What I see is that the first visible cluster of probes for the bottom dataset (on the screenshot) is not spread entirely over the read. Instead it covers the section of the read that met similar criteria for the top dataset.

I was thinking that this would produce probes with value 1 for the top dataset as it is currently seen. However for the bottom dataset I would have probes of value 1 for entire read, wider than that of top dataset.

Originally Posted by simonandrews

You could also put all of your data into a single data group and then treat it as a single dataset.

These are independent samples. So I would like to avoid doing that.

I hope I didn't make it more complicated.

[simonandrews]

Originally Posted by mjp

Sorry for this small delay and for not being specific enough. I see that the title of my post was wrong.
I wanted to have a list of probes that have just strand specific reads covering them. So the list would contain probes that have only fwd, only rvrs, none, and both type of reads covering them. Ideally I would like to have such a list for each of my stores independently. Is it possible to create that in single SeqMonk pipeline or does it have to be repeated for each store separately?

I don't think you can avoid having to do some filtering per-store since I think I understand that you want to end up with a separate set of lists for each store? You'll also need to do two quantitations. I reckon the quickest way to get these lists would be:

1) Quantitate using the difference quantitation using the option to quantitate forward reads as percentage of all reads.

2) Use the values filter to select probes with a value of 0 or 100 for each store. This should be pretty quick since you can set the parameters and then just select each store in turn in the filter and re-run it without having to reopen the dialog.

3) Requantitate the data with a simple read count quantitation (no corrections or transformations).

4) Use a values filter to select probes with some reads in them for each store.

5) Use the combine probes filter to select the subset of the 0% results which actually have some data in them. Again you can do this in a single filter session by changing which lists you're using so it shouldn't be too horrible.

Does this sound OK?

[mjp]

That does sound OK indeed. Thanks!

As another alternative I could do simple read count for all probes and create an annotated probe report #1 (not annotating with anything) for all the stores which would give me a list of probes with '0's for probes not having any reads over them.

Do the difference quantitation of forward as percentage of all, which would give me the list probes with '100' for the probes with only forward reads across all stores. Probe report #2.

Same for the reverse. Probe report #3.

Having these three reports it would be easy to parse it outside of SeqMonk.
If probe in #1 = 0 => no reads.
If probe in #1 > 0 and probe in #2 = 100 => then only forward
If probe in #1 > 0 and probe in #3 = 100 => then only reverse.
If probe in #1 > 0 and probes in #2 and #3 different that 0 and 100 => both reads

I think this way I will get what I need the fastest for all stores.

One way or another, your input about difference quantitation was invaluable.

Thanks again.

[shadow19c]

Hello,
I want to know how can I vizualise teh bedgrap file from bismark after methylation call?

Thanks

[simonandrews]

Originally Posted by shadow19c

Hello,
I want to know how can I vizualise the bedgraph file from bismark after methylation call?

SeqMonk is designed to to the quantitation of your data within the program rather than taking in externally quantitated files. Rather than trying to load the BedGraph file from Bismark you'd instead import the raw data from the methylation extractor and then quantitate this however you wanted inside SeqMonk to be able to visualise the methylation levels.

I put up a tutorial video covering some of the basics for working with bisulphite data on our youtube channel which should give you an idea how to get started with this.

[glados]

Dear Simon.

I sent you a private message a few weeks ago. Perhaps you can take a look at it? It was a question regarding installing a custom genome.