pRegion for all #172

lldelisle · 2020-01-09T16:55:14Z

Hi,
There was an argument in HiCMatrix to be able to load only part of the cool matrix. Unfortunately, through the different changes, this pRegion argument were not used anymore.
I put it back and I thought it could be usefull to add it to all other Tracks, for example, bed/bedgraph/links etc...
This should speed the use of pygenometracks on real data.
If it is merged and #162 also, the message in #162 should be adapted.

LeilyR

If I got you correctly, you try to load only certain part of matrix instead of all of it and then plot some regions in that part, right? As you have mentioned yourself this is only possible if .cool. I am not sure if from your changes I could see how you could handle that. With .h5 matrices the entire matrix needs to be loaded anyway. Or are pRegions basically the same regions which are given via --bed?

LeilyR · 2020-01-10T08:54:29Z

pygenometracks/tracks/HiCMatrixTrack.py

        try:
            self.hic_ma = HiCMatrix.hiCMatrix(self.properties['file'], pChrnameList=region)
        except Exception:
            region = [str(self.properties['region'][0]) + ':' + str(start) + '-' + str(self.properties['region'][2])]
-            self.hic_ma = HiCMatrix.hiCMatrix(self.properties['file'], pChrnameList=region)
+            try:


Is this necessary?

I thought that if you would give chrX instead of X it would not work but it worked, so you are right, I will remove this.

What I experienced is that if you give a list which has a length of 2 this exit.

From the code, it looks like it could fail:
https://github.com/deeptools/HiCMatrix/blob/cd54a2e7982bc0880e17536a44b09459620cb6db/hicmatrix/lib/cool.py#L106-L120
but I did not manage to make it fail maybe this is linked to the version of cooler.

I see how you mean. I think it means that it can only load one region at the time, this can be confirmed by @joachimwolff . Maybe if you have more than one region you need to do it recursively ?

I managed:
If the chromosome in the cool is UCSC style (chr1, chr2...) and you want to plot --region Y:2500000-2600000 you get:
Wrong chromosome format. Please check UCSC / ensembl notation. but as it is an exit you cannot except it...

So, I will remove the try/catch and check it is UCSC format as UCSC format works on ensembl format but not the contrary.
@joachimwolff, in a next release, could you transform all the exit in exception, so we could handle them?

I think it was not a bad idea to fail or at least warn the user about format. just adding chr at the beginning of every single chr can make troubles

I think we really need to be able to catch the errors, @joachimwolff do you want me to do a PR?
I realized that if the region to plot is the end of the chromosome, then it will exit and we have no way to deal with it (and this but in independent of this PR this is a current bug).
If you want to deal with matrices in both format (ucsc and ensembl) it would be good to be able to plot both...

My prefered solution would be:

update HiCMatrix to be able to catch the error and even better, check that the pChromList does not go over the chromosome size.

use try/Catch to deal with ucsc/ensembl and if it is not implemented in HiCMatrix, to deal with region to get matrix above the chr size.

have a good conversion UCSC/Ensembl because we will still have issues with contigs.

LeilyR · 2020-01-10T08:57:16Z

If I got you correctly, you try to load only certain part of matrix instead of all of it and then plot some regions in that part, right? As you have mentioned yourself this is only possible if .cool. I am not sure if from your changes I could see how you could handle that. With .h5 matrices the entire matrix needs to be loaded anyway. Or are pRegions basically the same regions which are given via --bed?

I went through your changes and I got the answer to my question, so you can ignore it.

lldelisle · 2020-01-10T08:59:26Z

Hi,

For the matrix, the restriction is already implemented in HiCMatrix with the argument pChrnameList:
https://github.com/deeptools/HiCMatrix/blob/cd54a2e7982bc0880e17536a44b09459620cb6db/hicmatrix/HiCMatrix.py#L30
That we were using:

pyGenomeTracks/pygenometracks/tracks/HiCMatrixTrack.py

Lines 99 to 105 in 8bbb22e

    
           # try to open with end region + depth to avoid triangle effect in the plot 
        
           # if it fails open it with given end region. 
        
           try: 
        
               self.hic_ma = HiCMatrix.hiCMatrix(self.properties['file'], pChrnameList=region) 
        
           except Exception: 
        
               region = [str(self.properties['region'][0]) + ':' + str(start) + '-' + str(self.properties['region'][2])] 
        
               self.hic_ma = HiCMatrix.hiCMatrix(self.properties['file'], pChrnameList=region)

However, the region variable was always equal to None because we were generating the tracks before looking at the region, thus the pRegion was not set:

pyGenomeTracks/pygenometracks/plotTracks.py

Line 272 in 8bbb22e

    
           trp = PlotTracks(args.tracks.name, args.width, fig_height=args.height, fontsize=args.fontSize, dpi=args.dpi, track_label_width=args.trackLabelFraction)

pyGenomeTracks/pygenometracks/tracksClass.py

Lines 64 to 69 in 8bbb22e

    
           class PlotTracks(object): 
        
               def __init__(self, tracks_file, fig_width=DEFAULT_FIGURE_WIDTH, 
        
                            fig_height=None, fontsize=None, dpi=None, 
        
                            track_label_width=None, 
        
                            pRegion=None):

So for the matrix, I just put back the pRegion which is the plotted region +/- the depth to avoid the triangle effect.

For the other format (bed, bedgraph and co, links):
Nothing was implemented but following advice from the gffutils guy (Is it possible to create a db only for one chromosome? daler/gffutils#143) I used pybedtools. I chose arbitrary to use the plotted region +/- 100kb

joachimwolff · 2020-01-10T09:50:18Z

If I got you correctly, you try to load only certain part of matrix instead of all of it and then plot some regions in that part, right? As you have mentioned yourself this is only possible if .cool. I am not sure if from your changes I could see how you could handle that. With .h5 matrices the entire matrix needs to be loaded anyway. Or are pRegions basically the same regions which are given via --bed?

HiCMatrix is implemented in a way that you simply pass the region you want to load and HiCMatrix takes care of it, independent if it is a cool file (and partial load is supported) or h5 (and HiCMatrix loads all and trims down to the requested region).

joachimwolff · 2020-01-10T09:57:32Z

@lldelisle I see you remove the UCSC/ensemble check. Do you cover UCSC/ensemble checking somewhere else?

LeilyR · 2020-01-10T09:57:36Z

pygenometracks/tracks/HiCMatrixTrack.py

+            # The chromosome name will be UCSC format because
+            # cooler v0.8.5 can fetch UCSC format in Ensembl-like format
+            # but not the contrary
+            if not chrom.startswith('chr'):


what if it is mitochondrial dna or some contigs?

I agree it would be better to be able to catch the exception from HiCMatrix but for the moment we cannot as it is exit.

I will try to see how it behave.

it is a very bad practice to add chr to just everything, that is why cooler doesn't do that too, i suppose. We already have tests to check for ensemble vs ucsc and i think that failure message in hicmatrix makes sense. This is my opinion. Maybe others have better idea.

What do you mean by

We already have tests to check for ensemble vs ucsc

lldelisle · 2020-01-10T10:12:47Z

@lldelisle I see you remove the UCSC/ensemble check. Do you cover UCSC/ensemble checking somewhere else?

I did not removed it is still here but latter on:

pyGenomeTracks/pygenometracks/tracks/HiCMatrixTrack.py

Lines 174 to 182 in 8bbb22e

    
           if chrom_region not in chrom_sizes: 
        
               chrom_region_before = chrom_region 
        
               chrom_region = self.change_chrom_names(chrom_region) 
        
               if chrom_region not in chrom_sizes: 
        
                   self.log.warning("*Warning*\nNeither " + chrom_region_before 
        
                                    + " nor " + chrom_region + " existss as a " 
        
                                    "chromosome name on the matrix. " 
        
                                    "This will generate an empty track!!\n") 
        
                   return

And this is still useful.

lldelisle · 2020-01-10T11:08:15Z

I was wrong, cooler cannot deal with UCSC vs Ensembl in none of the way but it raises exception with good messages:

import cooler

cooler_file = cooler.Cooler('./pygenometracks/tests/test_data/Li_et_al_2015.cool')
matrix = cooler_file.matrix(balance=False, sparse=True).fetch('chrX:0-1000').tocsr()
# ValueError: Unknown sequence label: chrX
matrix = cooler_file.matrix(balance=False, sparse=True).fetch('X:0-10000000000').tocsr()
# ValueError: Genomic region out of bounds: [0, 10000000000)

I will remove the lines I added, however, I still think the best solution is:

update HiCMatrix to be able to catch the error if possible keeping the message (or even better, check that the pChromList does not go over the chromosome size).
use try/Catch to deal with ucsc/ensembl and if it is not implemented in HiCMatrix, to deal with region to get matrix above the chr size.
have a good conversion UCSC/Ensembl because we will still have issues with contigs for this there is a repo https://github.com/dpryan79/ChromosomeMappings, I will try to see how to implement it.

lldelisle · 2020-01-10T13:37:30Z

I solved 3. in #174

bgruening · 2020-05-05T10:32:13Z

This one needs an extensive rebase :(

lldelisle · 2020-05-05T12:50:47Z

Yes but this one needs to work on HiCMatrix first.. So this is in standby...

…into pRegionForAll

lldelisle · 2020-06-03T20:51:47Z

Dear all,
Thanks to the update of HiCMatrix this PR is back on the track. I would like that @LeilyR, @joachimwolff and @bgruening (if you have time) accept before we merge it. There is no emergency but this will be a massive gain of time for all plots except the one using h5 but these files can be converted to cool.
What is still missing:

advice users to use cool instead of h5
advice users to do multiple plots instead of using --BED if they are not using h5.

lldelisle · 2020-06-04T05:50:03Z

Sorry, finally, I am still working on it and I will call you when it is ready.

…cept in HiC.

lldelisle · 2020-06-04T10:04:16Z

@LeilyR @joachimwolff @bgruening I am ready for review.
Here is a summary:
When the tracks are created, the behaviour is changed:

for bed, gtf: the file is first intersected with all plotted regions +/- 100kb (to get gene extremities) before being loaded.
for bedgraph and all children: narrowPeak, bedGraphMatrix, Epilogos, if it is not a tabix file, then the file is first intersected with all plotted regions sharp (no need to do +/- anything).
for links: the file is first intersected with the whole chromosome of all plotted regions (to be sure we are not missing anything as there are chrom1 start1 end1 and chrom2 start2 end2).
for matrix: if all plotted regions are on the same chromosome, a single region from the minimum of all start to the maximum of all ends extended of the depth is used in HiCMatrix. If it was a cool file, this will increase the speed like hell...

I think this is a massive improvement for people using bedgraph/bed or gtf/cool files.

bgruening · 2020-06-04T17:56:53Z

@lldelisle I will only have time over the weekend, sorry. You are way too fast for me :)

lldelisle · 2020-06-04T18:42:06Z

But weekend is more than fine.

…into pRegionForAll

lldelisle · 2020-06-05T11:12:19Z

Just to let you know these changes implies that most color scales where min_value and/or max_value were not specified will change because they will be evaluated on the plotted regions instead of the whole genome. But I don't think it is an issue.

bgruening

I'm not super happy with disabling the log messages for some lines of code, but I also don't have a better idea :(

Looks great @lldelisle ... feel free to merge.

lldelisle · 2020-07-05T19:39:23Z

Thanks for the review. Yes, for pybedtools, I don't like it neither... but the only alternative I see is to create only one temporary file and write into it all logs of bedtools and give the path at the beginning... Do you think it is better?

…into pRegionForAll

bgruening · 2020-07-05T19:52:17Z

I think it ok'ish like it is. Thanks @lldelisle!

lldelisle added 4 commits January 9, 2020 12:39

use pRegion for HiC tracks to speed up

404c264

added a test with a bed with multiple chr

28d5653

added pybedtools

7e38dfa

added pRegion everywhere there is an intervaltree

9d7a3a4

lldelisle mentioned this pull request Jan 9, 2020

Big files #162

Closed

LeilyR reviewed Jan 10, 2020

View reviewed changes

remove try/catch but use UCSC format

1cea057

LeilyR reviewed Jan 10, 2020

View reviewed changes

joachimwolff self-requested a review January 10, 2020 10:24

remove the chr added to the region

34cb06a

lldelisle added 5 commits January 20, 2020 17:23

use one cool in the test of hic rasterize

3c0a900

deal with cases where the matrix is empty on the region

3b051fa

change error to warning and trycatch the vmax as percentile

de5b0fd

linting

329b153

new test

20c929f

lldelisle changed the title ~~pRegion for all~~ [WIP] pRegion for all May 5, 2020

joachimwolff requested a review from LeilyR May 6, 2020 07:09

lldelisle added 3 commits June 3, 2020 16:02

Merge branch 'develop' of https://github.com/deeptools/pyGenomeTracks …

0a3a3da

…into pRegionForAll

remove ununsed matplotlib

bd88273

put back HUGE_NUMBER

815c8ea

lldelisle added 8 commits June 3, 2020 18:47

add test for bed vs region

16c76b1

remove y_axis when no data

9e9c062

adapt pChrnameList to exceptions

3c39a02

update and add new tests

40bc04b

remove debug prints

f68a2b7

require hicmatrix>=13

4f559b7

update doc

4144cc7

linting

0823bef

lldelisle changed the title ~~[WIP] pRegion for all~~ pRegion for all Jun 3, 2020

lldelisle added 4 commits June 3, 2020 22:58

remove redundant test

2318fa8

mv self.img=None on top

a140022

remove useless comments

fbf9eb5

linting

5357f7d

lldelisle added 4 commits June 4, 2020 10:28

change pRegion to plot_regions and intersect with all plot_regions ex…

c9a32f7

…cept in HiC.

update tests/doc to use imbricated bed

96fbd6b

linting

254c911

make prettier the coding

e46f72f

Merge branch 'develop' of https://github.com/deeptools/pyGenomeTracks …

a68dbbe

…into pRegionForAll

bgruening approved these changes Jul 4, 2020

View reviewed changes

Merge branch 'develop' of https://github.com/deeptools/pyGenomeTracks …

63e103d

…into pRegionForAll

bgruening merged commit b6cb3c6 into deeptools:develop Jul 5, 2020

lldelisle deleted the pRegionForAll branch July 10, 2020 14:28

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

pRegion for all #172

pRegion for all #172

lldelisle commented Jan 9, 2020

LeilyR left a comment •

edited

Loading

LeilyR Jan 10, 2020 •

edited

Loading

lldelisle Jan 10, 2020

lldelisle Jan 10, 2020

lldelisle Jan 10, 2020

LeilyR Jan 10, 2020 •

edited

Loading

lldelisle Jan 10, 2020

lldelisle Jan 10, 2020

LeilyR Jan 10, 2020 •

edited

Loading

lldelisle Jan 10, 2020

LeilyR commented Jan 10, 2020

lldelisle commented Jan 10, 2020

joachimwolff commented Jan 10, 2020

joachimwolff commented Jan 10, 2020

LeilyR Jan 10, 2020

lldelisle Jan 10, 2020

lldelisle Jan 10, 2020

LeilyR Jan 10, 2020

lldelisle Jan 10, 2020

lldelisle commented Jan 10, 2020

lldelisle commented Jan 10, 2020

lldelisle commented Jan 10, 2020

bgruening commented May 5, 2020

lldelisle commented May 5, 2020

lldelisle commented Jun 3, 2020

lldelisle commented Jun 4, 2020

lldelisle commented Jun 4, 2020

bgruening commented Jun 4, 2020

lldelisle commented Jun 4, 2020

lldelisle commented Jun 5, 2020

bgruening left a comment

lldelisle commented Jul 5, 2020

bgruening commented Jul 5, 2020

pRegion for all #172

pRegion for all #172

Conversation

lldelisle commented Jan 9, 2020

LeilyR left a comment • edited Loading

Choose a reason for hiding this comment

LeilyR Jan 10, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

LeilyR Jan 10, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

LeilyR Jan 10, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

LeilyR commented Jan 10, 2020

lldelisle commented Jan 10, 2020

joachimwolff commented Jan 10, 2020

joachimwolff commented Jan 10, 2020

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

lldelisle commented Jan 10, 2020

lldelisle commented Jan 10, 2020

lldelisle commented Jan 10, 2020

bgruening commented May 5, 2020

lldelisle commented May 5, 2020

lldelisle commented Jun 3, 2020

lldelisle commented Jun 4, 2020

lldelisle commented Jun 4, 2020

bgruening commented Jun 4, 2020

lldelisle commented Jun 4, 2020

lldelisle commented Jun 5, 2020

bgruening left a comment

Choose a reason for hiding this comment

lldelisle commented Jul 5, 2020

bgruening commented Jul 5, 2020

LeilyR left a comment •

edited

Loading

LeilyR Jan 10, 2020 •

edited

Loading

LeilyR Jan 10, 2020 •

edited

Loading

LeilyR Jan 10, 2020 •

edited

Loading