New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

-E option not working #75

Closed
yang-yangfeng opened this Issue May 12, 2017 · 4 comments

Comments

Projects
None yet
3 participants
@yang-yangfeng
Contributor

yang-yangfeng commented May 12, 2017

I pulled the latest master. It looks like there's a problem with the -E option. Without -E, it works fine.

yafeng$ regtools cis-splice-effects identify -E hcc1395/clinseq_7/hcc1395_filtered.vcf.gz hcc1395/clinseq_7/tumor_rna_alignments.bam GRCh37.fa /gscmnt/gc2602/griffithlab/regtools/GRCh37.87.exons.gtf

Program: regtools
Version: 0.4.0
Variant file: hcc1395/clinseq_7/hcc1395_filtered.vcf.gz
Alignment file: hcc1395/clinseq_7/tumor_rna_alignments.bam
Reference fasta file: GRCh37.fa
Annotation file: /gscmnt/gc2602/griffithlab/regtools/GRCh37.87.exons.gtf
Warning: The index file is older than the data file: hcc1395/clinseq_7/hcc1395_filtered.vcf.gz.tbi

Variant 1 985359 985360 -1
Variant region is 1:4294967295-0
Unable to iterate to region within BAM.yafeng$

@gatoravi

This comment has been minimized.

Show comment
Hide comment
@gatoravi

gatoravi May 18, 2017

Contributor

Hi @yanyangfeng could you try this when you get a chance, hopefully this fixes it..

Contributor

gatoravi commented May 18, 2017

Hi @yanyangfeng could you try this when you get a chance, hopefully this fixes it..

@yang-yangfeng

This comment has been minimized.

Show comment
Hide comment
@yang-yangfeng

yang-yangfeng May 19, 2017

Contributor

@gatoravi @malachig so I can at least run it now and superficially it seems to be working. I ran it on hcc1395 with the filtered list of variants and got 30535 junctions. Compare that to 291080 junctions from junctions extract and 221 junctions from cse identify without the -E option. My filtered input vcf had 35072 variants - of these, 1851 unique variants were associated with at least one reported junction (with the E option).

Scanning through the results, a few variant associations did make me scratch my head a bit. I used the default options besides for -E, so the window to look for associated junctions should be the default, which is one exon away. If you look at this set of junctions, all of them are associated with variant A (1:2286742-2286743) and some are associated with both A and B (1:2303900-2303901). (IGV session here: https://drive.google.com/open?id=0BzhxOtH9JL8vVk1GWlZjYTNmQTQ)

1 2286903 2288871 JUNC00000110 10 - GT-AG 1 0 1 DA 1 1 1 MORN1 ENST00000378529,ENST00000378531,ENST00000606372 1:2286742-2286743
1 2289037 2290031 JUNC00000111 9 - GT-AG 1 0 0 DA 1 1 1 MORN1 ENST00000378529,ENST00000378531,ENST00000606372 1:2286742-2286743
1 2290154 2303920 JUNC00000112 14 - GT-AG 2 1 3 DA 1 1 1 MORN1 ENST00000378529,ENST00000378531,ENST00000419785,ENST00000469374,ENST00000606372,ENST00000607342 1:2286742-2286743
1 2303344 2303920 JUNC00000113 44 - GT-AG 1 0 1 DA 1 1 1 MORN1 ENST00000378529,ENST00000378531,ENST00000419785,ENST00000469374,ENST00000606372,ENST00000607342 1:2286742-2286743,1:2303900-2303901
1 2304030 2305900 JUNC00000114 111 - GT-AG 1 0 0 DA 1 1 1 MORN1 ENST00000378529,ENST00000378531,ENST00000419785,ENST00000469374,ENST00000606372,ENST00000607342 1:2286742-2286743,1:2303900-2303901
1 2305996 2316417 JUNC00000115 128 - GT-AG 4 1 5 DA 1 1 1 MORN1 ENST00000378525,ENST00000378529,ENST00000378531,ENST00000419785,ENST00000449373,ENST00000469374,ENST00000606372,ENST00000607342 1:2286742-2286743,1:2303900-2303901
1 2316148 2316249 JUNC00000116 2 + GT-AG 0 0 0 N 0 0 0 NA NA 1:2286742-2286743,1:2303900-2303901
1 2316504 2317246 JUNC00000117 138 - GT-AG 3 0 1 DA 1 1 1 MORN1 ENST00000378525,ENST00000378529,ENST00000378531,ENST00000419785,ENST00000449373,ENST00000469374,ENST00000606372 1:2286742-2286743,1:2303900-2303901
1 2317336 2318735 JUNC00000118 6 - GT-AG 1 0 1 A 0 1 0 MORN1 ENST00000378525,ENST00000378529,ENST00000378531,ENST00000449373,ENST00000469374,ENST00000475812,ENST00000606372 1:2286742-2286743,1:2303900-2303901
1 2317336 2318825 JUNC00000119 6 - GT-AG 1 0 1 A 0 1 0 MORN1 ENST00000378525,ENST00000378529,ENST00000378531,ENST00000449373,ENST00000469374,ENST00000475812,ENST00000606372 1:2286742-2286743,1:2303900-2303901
1 2317336 2318858 JUNC00000120 187 - GT-AG 1 0 1 DA 1 1 1 MORN1 ENST00000378525,ENST00000378529,ENST00000378531,ENST00000449373,ENST00000469374,ENST00000475812,ENST00000606372 1:2286742-2286743,1:2303900-2303901
1 2318968 2319678 JUNC00000121 190 - GT-AG 1 0 1 DA 1 1 1 MORN1 ENST00000378525,ENST00000378529,ENST00000378531,ENST00000449373,ENST00000469374,ENST00000475812,ENST00000494279,ENST00000606372 1:2286742-2286743
1 2319776 2321364 JUNC00000122 191 - GT-AG 1 0 0 DA 1 1 1 MORN1 ENST00000378525,ENST00000378529,ENST00000378531,ENST00000449373,ENST00000494279,ENST00000606372 1:2286742-2286743
1 2319776 2322897 JUNC00000123 26 - GT-AG 2 1 2 DA 1 1 1 MORN1 ENST00000378525,ENST00000378529,ENST00000378531,ENST00000449373,ENST00000494279,ENST00000606372 1:2286742-2286743
1 2321435 2322897 JUNC00000124 151 - GT-AG 1 0 1 DA 1 1 1 MORN1 ENST00000378525,ENST00000378529,ENST00000378531,ENST00000494279,ENST00000606372 1:2286742-2286743

It doesn't quite make sense to me how regtools is associating the junctions with variants in some of these examples. For instance, both blocks of JUNC113 are clearly one exon away from variant B in a number of transcripts, but there is no transcript for which either block of JUNC113 is one exon away from variant A. It seems possible that since one of the transcripts associated with JUNC113 ends with the 5' exon of this junction, that regtools thinks that this junction is one variant away since there are no intervening exons. But that wouldn't explain why JUNC124 is only associated with variant A. Anyway, perhaps we could take a look at this together. I don't immediately see why adding the -E option would have caused this issue, so perhaps it was an existing issue.

Contributor

yang-yangfeng commented May 19, 2017

@gatoravi @malachig so I can at least run it now and superficially it seems to be working. I ran it on hcc1395 with the filtered list of variants and got 30535 junctions. Compare that to 291080 junctions from junctions extract and 221 junctions from cse identify without the -E option. My filtered input vcf had 35072 variants - of these, 1851 unique variants were associated with at least one reported junction (with the E option).

Scanning through the results, a few variant associations did make me scratch my head a bit. I used the default options besides for -E, so the window to look for associated junctions should be the default, which is one exon away. If you look at this set of junctions, all of them are associated with variant A (1:2286742-2286743) and some are associated with both A and B (1:2303900-2303901). (IGV session here: https://drive.google.com/open?id=0BzhxOtH9JL8vVk1GWlZjYTNmQTQ)

1 2286903 2288871 JUNC00000110 10 - GT-AG 1 0 1 DA 1 1 1 MORN1 ENST00000378529,ENST00000378531,ENST00000606372 1:2286742-2286743
1 2289037 2290031 JUNC00000111 9 - GT-AG 1 0 0 DA 1 1 1 MORN1 ENST00000378529,ENST00000378531,ENST00000606372 1:2286742-2286743
1 2290154 2303920 JUNC00000112 14 - GT-AG 2 1 3 DA 1 1 1 MORN1 ENST00000378529,ENST00000378531,ENST00000419785,ENST00000469374,ENST00000606372,ENST00000607342 1:2286742-2286743
1 2303344 2303920 JUNC00000113 44 - GT-AG 1 0 1 DA 1 1 1 MORN1 ENST00000378529,ENST00000378531,ENST00000419785,ENST00000469374,ENST00000606372,ENST00000607342 1:2286742-2286743,1:2303900-2303901
1 2304030 2305900 JUNC00000114 111 - GT-AG 1 0 0 DA 1 1 1 MORN1 ENST00000378529,ENST00000378531,ENST00000419785,ENST00000469374,ENST00000606372,ENST00000607342 1:2286742-2286743,1:2303900-2303901
1 2305996 2316417 JUNC00000115 128 - GT-AG 4 1 5 DA 1 1 1 MORN1 ENST00000378525,ENST00000378529,ENST00000378531,ENST00000419785,ENST00000449373,ENST00000469374,ENST00000606372,ENST00000607342 1:2286742-2286743,1:2303900-2303901
1 2316148 2316249 JUNC00000116 2 + GT-AG 0 0 0 N 0 0 0 NA NA 1:2286742-2286743,1:2303900-2303901
1 2316504 2317246 JUNC00000117 138 - GT-AG 3 0 1 DA 1 1 1 MORN1 ENST00000378525,ENST00000378529,ENST00000378531,ENST00000419785,ENST00000449373,ENST00000469374,ENST00000606372 1:2286742-2286743,1:2303900-2303901
1 2317336 2318735 JUNC00000118 6 - GT-AG 1 0 1 A 0 1 0 MORN1 ENST00000378525,ENST00000378529,ENST00000378531,ENST00000449373,ENST00000469374,ENST00000475812,ENST00000606372 1:2286742-2286743,1:2303900-2303901
1 2317336 2318825 JUNC00000119 6 - GT-AG 1 0 1 A 0 1 0 MORN1 ENST00000378525,ENST00000378529,ENST00000378531,ENST00000449373,ENST00000469374,ENST00000475812,ENST00000606372 1:2286742-2286743,1:2303900-2303901
1 2317336 2318858 JUNC00000120 187 - GT-AG 1 0 1 DA 1 1 1 MORN1 ENST00000378525,ENST00000378529,ENST00000378531,ENST00000449373,ENST00000469374,ENST00000475812,ENST00000606372 1:2286742-2286743,1:2303900-2303901
1 2318968 2319678 JUNC00000121 190 - GT-AG 1 0 1 DA 1 1 1 MORN1 ENST00000378525,ENST00000378529,ENST00000378531,ENST00000449373,ENST00000469374,ENST00000475812,ENST00000494279,ENST00000606372 1:2286742-2286743
1 2319776 2321364 JUNC00000122 191 - GT-AG 1 0 0 DA 1 1 1 MORN1 ENST00000378525,ENST00000378529,ENST00000378531,ENST00000449373,ENST00000494279,ENST00000606372 1:2286742-2286743
1 2319776 2322897 JUNC00000123 26 - GT-AG 2 1 2 DA 1 1 1 MORN1 ENST00000378525,ENST00000378529,ENST00000378531,ENST00000449373,ENST00000494279,ENST00000606372 1:2286742-2286743
1 2321435 2322897 JUNC00000124 151 - GT-AG 1 0 1 DA 1 1 1 MORN1 ENST00000378525,ENST00000378529,ENST00000378531,ENST00000494279,ENST00000606372 1:2286742-2286743

It doesn't quite make sense to me how regtools is associating the junctions with variants in some of these examples. For instance, both blocks of JUNC113 are clearly one exon away from variant B in a number of transcripts, but there is no transcript for which either block of JUNC113 is one exon away from variant A. It seems possible that since one of the transcripts associated with JUNC113 ends with the 5' exon of this junction, that regtools thinks that this junction is one variant away since there are no intervening exons. But that wouldn't explain why JUNC124 is only associated with variant A. Anyway, perhaps we could take a look at this together. I don't immediately see why adding the -E option would have caused this issue, so perhaps it was an existing issue.

@yang-yangfeng

This comment has been minimized.

Show comment
Hide comment
@yang-yangfeng

yang-yangfeng Aug 24, 2017

Contributor

Seems like there's definitely something weird going on with the -E option. Below are the results for two runs of cse-identify using a gtf file that only contains transcript ENST00000378529. 1) "-e 500" to isolate variant at 1:2286742-2286743 and defaults for window size and everything else. 2) "-E" to isolate variant at 1:2286742-2286743 (the other variant in this area does not overlap with any of the exons in ENST00000378529, so it's not being detected) and again, defaults for window size and everything else. Theoretically, these jobs should give identical results. However, as you can see, the junctions reported are vastly different and the variant regions, which I take to mean the window around the variant to look for junctions in, are also different. Seems like the -e 500 job is doing something more reasonable, so I think it's the -E option that is screwing things up. I'll do some work to clarify exactly what regtools is expected to do/doing and I'll dig into the code to see if I can see what's up with the -E option.

~$ cd /gscmnt/gc2602/griffithlab/regtools/yafeng
yafeng$ grep -E 'ENST00000378529' ../GRCh37.87.exons.chr1.gtf > test.gtf
yafeng$ /gscmnt/gc2502/griffithlab/yafeng/regtools/build/regtools cis-splice-effects identify -e 500 hcc1395/clinseq_7/hcc1395_filtered.vcf.gz hcc1395/clinseq_7/tumor_rna_alignments.bam GRCh37.fa /gscmnt/gc2602/griffithlab/regtools/yafeng/test.gtf

Program:	regtools
Version:	0.4.0
Variant file: hcc1395/clinseq_7/hcc1395_filtered.vcf.gz
Alignment file: hcc1395/clinseq_7/tumor_rna_alignments.bam
Reference fasta file: GRCh37.fa
Annotation file: /gscmnt/gc2602/griffithlab/regtools/yafeng/test.gtf
Warning: The index file is older than the data file: hcc1395/clinseq_7/hcc1395_filtered.vcf.gz.tbi


Variant 1	2286742	2286743	127		
Variant region is 1:2286616-2289037

chrom	start	end	name	score	strand	splice_site	acceptors_skipped	exons_skipped	donors_skipped	anchor	known_donor	known_acceptor	known_junction	genes	transcripts	variant_info
1	2286903	2288871	JUNC00000001	10	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
yafeng$ /gscmnt/gc2502/griffithlab/yafeng/regtools/build/regtools cis-splice-effects identify -E hcc1395/clinseq_7/hcc1395_filtered.vcf.gz hcc1395/clinseq_7/tumor_rna_alignments.bam GRCh37.fa /gscmnt/gc2602/griffithlab/regtools/yafeng/test.gtf

Program:	regtools
Version:	0.4.0
Variant file: hcc1395/clinseq_7/hcc1395_filtered.vcf.gz
Alignment file: hcc1395/clinseq_7/tumor_rna_alignments.bam
Reference fasta file: GRCh37.fa
Annotation file: /gscmnt/gc2602/griffithlab/regtools/yafeng/test.gtf
Warning: The index file is older than the data file: hcc1395/clinseq_7/hcc1395_filtered.vcf.gz.tbi


Variant 1	2286742	2286743	127		
Variant region is 1:2286616-2321435

chrom	start	end	name	score	strand	splice_site	acceptors_skipped	exons_skipped	donors_skipped	anchor	known_donor	known_acceptor	known_junction	genes	transcripts	variant_info
1	2286903	2288871	JUNC00000001	10	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
1	2289037	2290031	JUNC00000002	9	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
1	2290154	2303920	JUNC00000003	14	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
1	2303344	2303920	JUNC00000004	44	-	GT-AG	0	0	1	D	1	0	0	MORN1	ENST00000378529	1:2286742-2286743
1	2304030	2305900	JUNC00000005	111	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
1	2305996	2316417	JUNC00000006	128	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
1	2316148	2316249	JUNC00000007	2	+	GT-AG	0	0	0	N	0	0	0	NA	NA	1:2286742-2286743
1	2316504	2317246	JUNC00000008	138	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
1	2317336	2318735	JUNC00000009	6	-	GT-AG	1	0	0	A	0	1	0	MORN1	ENST00000378529	1:2286742-2286743
1	2317336	2318825	JUNC00000010	6	-	GT-AG	1	0	0	A	0	1	0	MORN1	ENST00000378529	1:2286742-2286743
1	2317336	2318858	JUNC00000011	187	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
1	2318968	2319678	JUNC00000012	190	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
1	2319776	2321364	JUNC00000013	191	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
Contributor

yang-yangfeng commented Aug 24, 2017

Seems like there's definitely something weird going on with the -E option. Below are the results for two runs of cse-identify using a gtf file that only contains transcript ENST00000378529. 1) "-e 500" to isolate variant at 1:2286742-2286743 and defaults for window size and everything else. 2) "-E" to isolate variant at 1:2286742-2286743 (the other variant in this area does not overlap with any of the exons in ENST00000378529, so it's not being detected) and again, defaults for window size and everything else. Theoretically, these jobs should give identical results. However, as you can see, the junctions reported are vastly different and the variant regions, which I take to mean the window around the variant to look for junctions in, are also different. Seems like the -e 500 job is doing something more reasonable, so I think it's the -E option that is screwing things up. I'll do some work to clarify exactly what regtools is expected to do/doing and I'll dig into the code to see if I can see what's up with the -E option.

~$ cd /gscmnt/gc2602/griffithlab/regtools/yafeng
yafeng$ grep -E 'ENST00000378529' ../GRCh37.87.exons.chr1.gtf > test.gtf
yafeng$ /gscmnt/gc2502/griffithlab/yafeng/regtools/build/regtools cis-splice-effects identify -e 500 hcc1395/clinseq_7/hcc1395_filtered.vcf.gz hcc1395/clinseq_7/tumor_rna_alignments.bam GRCh37.fa /gscmnt/gc2602/griffithlab/regtools/yafeng/test.gtf

Program:	regtools
Version:	0.4.0
Variant file: hcc1395/clinseq_7/hcc1395_filtered.vcf.gz
Alignment file: hcc1395/clinseq_7/tumor_rna_alignments.bam
Reference fasta file: GRCh37.fa
Annotation file: /gscmnt/gc2602/griffithlab/regtools/yafeng/test.gtf
Warning: The index file is older than the data file: hcc1395/clinseq_7/hcc1395_filtered.vcf.gz.tbi


Variant 1	2286742	2286743	127		
Variant region is 1:2286616-2289037

chrom	start	end	name	score	strand	splice_site	acceptors_skipped	exons_skipped	donors_skipped	anchor	known_donor	known_acceptor	known_junction	genes	transcripts	variant_info
1	2286903	2288871	JUNC00000001	10	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
yafeng$ /gscmnt/gc2502/griffithlab/yafeng/regtools/build/regtools cis-splice-effects identify -E hcc1395/clinseq_7/hcc1395_filtered.vcf.gz hcc1395/clinseq_7/tumor_rna_alignments.bam GRCh37.fa /gscmnt/gc2602/griffithlab/regtools/yafeng/test.gtf

Program:	regtools
Version:	0.4.0
Variant file: hcc1395/clinseq_7/hcc1395_filtered.vcf.gz
Alignment file: hcc1395/clinseq_7/tumor_rna_alignments.bam
Reference fasta file: GRCh37.fa
Annotation file: /gscmnt/gc2602/griffithlab/regtools/yafeng/test.gtf
Warning: The index file is older than the data file: hcc1395/clinseq_7/hcc1395_filtered.vcf.gz.tbi


Variant 1	2286742	2286743	127		
Variant region is 1:2286616-2321435

chrom	start	end	name	score	strand	splice_site	acceptors_skipped	exons_skipped	donors_skipped	anchor	known_donor	known_acceptor	known_junction	genes	transcripts	variant_info
1	2286903	2288871	JUNC00000001	10	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
1	2289037	2290031	JUNC00000002	9	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
1	2290154	2303920	JUNC00000003	14	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
1	2303344	2303920	JUNC00000004	44	-	GT-AG	0	0	1	D	1	0	0	MORN1	ENST00000378529	1:2286742-2286743
1	2304030	2305900	JUNC00000005	111	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
1	2305996	2316417	JUNC00000006	128	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
1	2316148	2316249	JUNC00000007	2	+	GT-AG	0	0	0	N	0	0	0	NA	NA	1:2286742-2286743
1	2316504	2317246	JUNC00000008	138	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
1	2317336	2318735	JUNC00000009	6	-	GT-AG	1	0	0	A	0	1	0	MORN1	ENST00000378529	1:2286742-2286743
1	2317336	2318825	JUNC00000010	6	-	GT-AG	1	0	0	A	0	1	0	MORN1	ENST00000378529	1:2286742-2286743
1	2317336	2318858	JUNC00000011	187	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
1	2318968	2319678	JUNC00000012	190	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
1	2319776	2321364	JUNC00000013	191	-	GT-AG	1	0	0	DA	1	1	1	MORN1	ENST00000378529	1:2286742-2286743
@yang-yangfeng

This comment has been minimized.

Show comment
Hide comment
@yang-yangfeng

yang-yangfeng Sep 7, 2017

Contributor

With the -e option, window (variant region) seems to be what we expect, which is the bounds from the central (casette) exon in which the variant was found to the exon downstream (and potentially the variant upstream).

With the -E option, it looks like it's grabbing the entire transcript as the window.

Contributor

yang-yangfeng commented Sep 7, 2017

With the -e option, window (variant region) seems to be what we expect, which is the bounds from the central (casette) exon in which the variant was found to the exon downstream (and potentially the variant upstream).

With the -E option, it looks like it's grabbing the entire transcript as the window.

gatoravi added a commit to gatoravi/regtools that referenced this issue Sep 11, 2017

yang-yangfeng added a commit that referenced this issue Oct 9, 2017

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment