Skip to content
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

Missed call - not clear why #491

Open
gevro opened this issue Aug 29, 2023 · 11 comments
Open

Missed call - not clear why #491

gevro opened this issue Aug 29, 2023 · 11 comments

Comments

@gevro
Copy link

gevro commented Aug 29, 2023

Hi,
We did RNA-seq on a sample on which we know there is a very clear and strong splicing defect, as a positive control for the DROP pipeline.

The exon skipping is very obvious:
Screenshot 2023-08-28 at 9 37 43 PM

However, DROP did not detect it.

Here are our config parameters:


aberrantSplicing:
    run: true
    groups:
        - GTEX100
    recount: false
    longRead: false
    keepNonStandardChrs: false
    filter: true
    minExpressionInOneSample: 20
    minDeltaPsi: 0.05
    implementation: PCA
    padjCutoff: 0.1
    maxTestedDimensionProportion: 6
    genesToTest: null
    FRASER_version: "FRASER2"
    deltaPsiCutoff : 0.1
    quantileForFiltering: 0.75

Here is the call in results_gene_all.tsv, but it was not in the final FRASER output:

seqnames	start	end	width	strand	sampleID	hgncSymbol	type	pValue	psiValue	deltaPsi	counts	totalCounts	meanCounts	meanTotalCounts	nonsplitCounts	nonsplitProportion	nonsplitProportion_99quantile	annotatedJunction	pValueGene	padjustGene
**	**	**	3650	+	Sample_RNA	**	jaccard	0.032264	0.31	-0.14	146	464	123.63	128.96	59	0.13	NA	both	0.90339	1

Can you help us understand why this did not get into the final call set and how we can fix this?

Also would appreciate to know if there is a fix, how we can do that without rerunning the whole pipeline, which takes a long time.

Thanks

@vyepez88
Copy link
Collaborator

Hi,
Thanks for using DROP and reporting this.
It would be interesting to see the results for all 3 of the junctions involved in the exon skipping event you are focusing on here. Can you run the following command to get the results table of all the junctions of your sample of interest?

results(fds, sampleIDs = {your sample of interest},  aggregate = FALSE,  all = TRUE)

Also, from the one junction that you provided the results, it's hard for us to understand what might be happening. Can you provide the output of the plotExpectedVsObservedPsi function which can be a starting point for understanding why FRASER reports a relatively small, non-significant delta jaccard value in this case?

@gevro
Copy link
Author

gevro commented Sep 16, 2023

Sorry, none of these work. I'm getting these errors:

> results(fds,sampleIDs = "Sample_RNA",  aggregate = FALSE,  all = TRUE)
Sat Sep 16 16:06:31 2023: Collecting results for: psi3
Error: BiocParallel errors
  element index: 1, 2, 3
  first error: error in evaluating the argument 'x' in selecting a method for function 'rowMeans': the supplied 'dimnames' must have one list element per dimension
> plotExpectedVsObservedPsi(fds,type="psi5")
Error in normarg_dimnames(dimnames, seed_dim) : 
  the supplied 'dimnames' must have one list element per dimension

@AtaJadidAhari
Copy link
Collaborator

Hi @gevro, do you have any updates on this issue? You could maybe try again with DROP 1.4.0. I suspect this would fix the issue with rowMeans that you encountered when running results function and we could see how the results for the missed junctions look like.

@gevro
Copy link
Author

gevro commented Oct 4, 2024

We reran with 1.4.0 and we will check to see if it is now detected.

@gevro
Copy link
Author

gevro commented Oct 7, 2024

We checked and this splice event is still missed in drop 1.4.0.

@AtaJadidAhari
Copy link
Collaborator

AtaJadidAhari commented Oct 7, 2024

Then let's try to redo the results section and see what's going on. First, read in the fds object:

# read in the snakemake from step 08_extract_results_FRASER.R 
# tmp_dir is probably {path_to_where_you_run_drop}/.drop/tmp/ 

snakemake <- readRDS("{tmp_dir}/AS/{dataset}--{annotation}/08_results.Rds")
annotation    <- snakemake@wildcards$annotation
dataset    <- snakemake@wildcards$dataset
workingDir <- snakemake@params$workingDir

fds <- loadFraserDataSet(dir=workingDir, name=paste(dataset, annotation, sep = '--'))

# check the results for the sample
results(fds, sampleIDs = {your sample of interest},  aggregate = FALSE,  all = TRUE)

And check the metrics for the junctions that are missed. This time probably it will run through since you updated DROP as opposed to the previous try one year ago. Let us how how the results looks like!

@gevro
Copy link
Author

gevro commented Oct 7, 2024

What R libraries do I need to load to run this?

@vyepez88
Copy link
Collaborator

vyepez88 commented Oct 7, 2024

only FRASER

@gevro
Copy link
Author

gevro commented Oct 7, 2024

Which RDS object?

./processed_results/aberrant_splicing/datasets/savedObjects/GTEX100--v32/fds-object.RDS
./processed_data/aberrant_splicing/datasets/savedObjects/GTEX100/fds-object.RDS
./processed_data/aberrant_splicing/datasets/savedObjects/raw-local-GTEX100/fds-object.RDS

@gevro
Copy link
Author

gevro commented Oct 7, 2024

Getting errors

snakemake@wildcards
Error: no slot of name "wildcards" for this object of class "FraserDataSet"

@gevro
Copy link
Author

gevro commented Oct 7, 2024

Ok, figured out how to pull out the data. Here is the data for that junction. Why wasn't this detected? It is very obvious in the raw reads. Note, I redacted the location of the splice junction and sampleID and hgncSymbol. I can e-mail separately if helpful.

seqnames	start	end	width	strand	sampleID	hgncSymbol	type	pValue	padjust	psiValue	deltaPsi	counts	totalCounts	meanCounts	meanTotalCounts	nonsplitCounts	nonsplitProportion	nonsplitProportion_99quantile	annotatedJunction
**	**	**	7559	+	**	**	jaccard	0.92387	1	0.47	0.02	259	556	1.8	301.55	59	0.11	0.14	start

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants