alexyermanos / Platypus

R package for the analysis of single-cell immune repertoires
GNU General Public License v3.0
37 stars 16 forks source link

Error GSEA function #14

Closed beginner984 closed 2 years ago

beginner984 commented 2 years ago

Hello

I am getting a permanent error in two of the functions

>   gsea_EAE <- GEX_GSEA(GEX.cluster.genes.output = gene_expression_cluster[[5]], MT.Rb.filter = T, path.to.pathways = "/user/data/data1/")
  Error in path.to.pathways[[2]] : subscript out of bounds
  In addition: Warning message:
    In GEX_GSEA(GEX.cluster.genes.output = gene_expression_cluster[[5]],  :
                  No genes matching filter found
                >   gsea_EAE <- GEX_GSEA(GEX.cluster.genes.output = gene_expression_cluster[[4]], MT.Rb.filter = T, path.to.pathways = "/user/data/data1/")
                Error in path.to.pathways[[2]] : subscript out of bounds
                In addition: Warning message:
                  In GEX_GSEA(GEX.cluster.genes.output = gene_expression_cluster[[4]],  :
                                No genes matching filter found
                              > 

I am sure c1.all.v7.4.symbols.gmt is in the path

alexyermanos commented 2 years ago

can you please post the first few lines of gene_expression_cluster[[5]] ?

can you also try setting the Mt.Rb.filter to F?

also if you have another findmarkers output can you test if this is working ? (e.g. gene_expression_cluster[[1]])

beginner984 commented 2 years ago

Thank you very much

I get this

> head(gene_expression_cluster[[5]])
p_val avg_logFC pct.1 pct.2 p_val_adj SYMBOL cluster
CD8A      0  1.630774 0.465 0.063         0   CD8A       4
GZMK      0  2.792085 0.734 0.034         0   GZMK       4
GZMA      0  1.461000 0.649 0.103         0   GZMA       4
CTSW      0  1.350435 0.752 0.185         0   CTSW       4
KLRB1     0  2.361302 0.686 0.144         0  KLRB1       4
CCL5      0  2.178712 0.939 0.175         0   CCL5       4
> 

> gsea_EAE <- GEX_GSEA(GEX.cluster.genes.output = gene_expression_cluster[[5]], MT.Rb.filter = FALSE, path.to.pathways = "/user/data/data1/")
Error in path.to.pathways[[2]] : subscript out of bounds
> gsea_EAE <- GEX_GSEA(GEX.cluster.genes.output = gene_expression_cluster[[0]], MT.Rb.filter = FALSE, path.to.pathways ="/user/data/data1/")
Error in gene_expression_cluster[[0]] : 
  attempt to select less than one element in get1index <real>
  > gsea_EAE <- GEX_GSEA(GEX.cluster.genes.output = gene_expression_cluster[[1]], MT.Rb.filter = FALSE, path.to.pathways ="/user/data/data1/")
Error in path.to.pathways[[2]] : subscript out of bounds
> gsea_EAE <- GEX_GSEA(GEX.cluster.genes.output = gene_expression_cluster[[2]], MT.Rb.filter = FALSE, path.to.pathways = "/user/data/data1/")
Error in path.to.pathways[[2]] : subscript out of bounds
>
vickreiner commented 2 years ago

Hi! I was able to replicate this issue on my machine. For reading the .gmt file, we are using fgsea::gmtPathways function. This one needs a direct path to the .gmt object. In your case this would probably look like path.to.pathways = "/user/data/data1/c1.all.v7.4.symbols.gmt"

I hope this solves the issue! Thanks!

beginner984 commented 2 years ago

Thank you so much Another error :(

> gsea_EAE <- GEX_GSEA(GEX.cluster.genes.output = gene_expression_cluster[[2]], MT.Rb.filter = F, path.to.pathways = "/user/data/data1/c1.all.v7.4.symbols.gmt")
Error in `dplyr::arrange()`:
  ! Problem with the implicit `transmute()` step.
x Problem while computing `..1 = IRanges::desc(NES)`.
Caused by error:
  ! unable to find an inherited method for function ‘desc’ for signature ‘"numeric"’
Run `rlang::last_error()` to see where the error occurred.
vickreiner commented 2 years ago

Hi! I assume this is a version incompatibility, but to make sure we can still check the raw Fgsea output. If you run the function with verbose = TRUE you should get raw results printed to the console. Could you share theses? Thanks a lot!

raphkuhn commented 2 years ago

Thank you for reporting this second error! A wrong package was assigned to the dplyr::desc() function. We fixed the function now and it can be downloaded from our github: https://github.com/alexyermanos/Platypus/blob/master/R/GEX_GSEA.R

beginner984 commented 2 years ago

Thank you very much

Sorry, might be irrelevant to the main post, but, I only have TCR + GEX

But I am wondering how I get B cells here

clonal_out <- VDJ_clonal_expansion(VDJ = vgm[[1]],celltype = "Bcells",clones = "30", group.by = "sample_id", color.by = "isotype", isotypes.to.plot = "all", treat.incomplete.clones = "exclude", treat.incomplete.cells = "proportional")
  #group by specifies how many separate plots should be generated. If vgm contains global clonotype information this can be set to "none"
  print(clonal_out[[1]])

Rplot35

And when I try for T cells I get nothing

clonal_out <- VDJ_clonal_expansion(VDJ = vgm[[1]],celltype = "Tcells",clones = "30", group.by = "sample_id", color.by = "isotype", isotypes.to.plot = "all", treat.incomplete.clones = "exclude", treat.incomplete.cells = "proportional")
  #group by specifies how many separate plots should be generated. If vgm contains global clonotype information this can be set to "none"
  print(clonal_out[[1]])

Rplot36

alexyermanos commented 2 years ago

Thanks for posting this - we will look into this but given that the distribution is identical for both arguments and there is no isotype information, the clonal distribution is most likely only describing the T cell expansion in both plots. this argument was more designed in the case that B and T cells are present in the same VGM

beginner984 commented 2 years ago

Sorry maybe this is a non related question

But what does the color scale bar here means?

Red means we have more fold of a given gene or more cells express this given gene?

Picture 1

alexyermanos commented 2 years ago

We would need to know the parameters used to call this function - it would either be the fraction of clones or fraction of cells using a given TRBV-TRVA pairing in this sample's repertoire

beginner984 commented 1 year ago

Excuse me,

What does this error mean please?


Integrating VDJ and GEX 
Error in VDJ_GEX_matrix(VDJ.out.directory.list = c(VDJ.out.directory.list_T,  : 
  object 'batches' not found

Here is the full thing


>                       vgm <- VDJ_GEX_matrix(VDJ.out.directory.list = c(VDJ.out.directory.list_T, VDJ.out.directory.list_B),
                                              +                                             GEX.out.directory.list = GEX.out.directory.list,
                                              +                                             Seurat.in = pbmc,
                                              +                                             GEX.integrate = T,
                                              +                                             VDJ.combine = T,
                                              +                                             integrate.GEX.to.VDJ = T,
                                              +                                             integrate.VDJ.to.GEX = T, #This will adjunct the VDJ information as metadata to the GEX object
                                              +                                             exclude.GEX.not.in.VDJ = F,
                                              +                                             filter.overlapping.barcodes.GEX = T,
                                              +                                             filter.overlapping.barcodes.VDJ = T,
                                              +                                             #exclude.on.cell.state.markers = c("CD3E"), #Exclude T cells from this analysis
                                                +                                             get.VDJ.stats = T,
                                              +                                             parallel.processing = "mclapply", #see note at the end of this chunk
                                              +                                             trim.and.align = T, #Do not align BCR sequences to reference 
                                              +                                             group.id = c(1,1,2,2),
                                              +                                             n.feature.rna=2000,
                                              +                                             n.count.rna.min=200,
                                              +                                             n.count.rna.max=2500,
                                              +                                             mito.filter = 5,
                                                                                             integration.method = "harmony") 
17:25:03 
Loading in data
17:25:45 Loaded VDJ data

Getting VDJ GEX stats
Starting with 1 of 12
Starting with 2 of 12
Starting with 3 of 12
Starting with 4 of 12
Starting with 5 of 12
Starting with 6 of 12
Starting with 7 of 12
Starting with 8 of 12
Starting with 9 of 12
Starting with 10 of 12
Starting with 11 of 12
Starting with 12 of 12
Getting 10x stats
17:25:53 Got VDJ GEX stats

For input Seurat object: 14261 cells assigned barcodes in GEX
For input Seurat object GEX and VDJ barcode overlap is: 9141

For sample 1: 2293 cells assigned with high confidence barcodes in VDJ

For sample 2: 2074 cells assigned with high confidence barcodes in VDJ

For sample 3: 2393 cells assigned with high confidence barcodes in VDJ

For sample 4: 100 cells assigned with high confidence barcodes in VDJ

For sample 5: 2077 cells assigned with high confidence barcodes in VDJ

For sample 6: 606 cells assigned with high confidence barcodes in VDJ

For sample 7: 1038 cells assigned with high confidence barcodes in VDJ

For sample 8: 341 cells assigned with high confidence barcodes in VDJ

For sample 9: 314 cells assigned with high confidence barcodes in VDJ

For sample 10: 28 cells assigned with high confidence barcodes in VDJ

For sample 11: 678 cells assigned with high confidence barcodes in VDJ

For sample 12: 106 cells assigned with high confidence barcodes in VDJ

Removed a total of 184 cells with non unique barcodes in VDJ
17:25:53  Starting VDJ barcode iteration 1 of 12...
Started mcapply cluster with 48 cores
17:30:19 Done with 1 of 12
17:30:19  Starting VDJ barcode iteration 2 of 12...
Started mcapply cluster with 48 cores
17:30:50 Done with 2 of 12
17:30:50  Starting VDJ barcode iteration 3 of 12...
Started mcapply cluster with 48 cores
17:35:59 Done with 3 of 12
17:35:59  Starting VDJ barcode iteration 4 of 12...
Started mcapply cluster with 48 cores
17:36:12 Done with 4 of 12
17:36:12  Starting VDJ barcode iteration 5 of 12...
Started mcapply cluster with 48 cores
17:40:18 Done with 5 of 12
17:40:18  Starting VDJ barcode iteration 6 of 12...
Started mcapply cluster with 48 cores
17:41:37 Done with 6 of 12
17:41:37  Starting VDJ barcode iteration 7 of 12...
Started mcapply cluster with 48 cores
17:43:48 Done with 7 of 12
17:43:48  Starting VDJ barcode iteration 8 of 12...
Started mcapply cluster with 48 cores
17:44:31 Done with 8 of 12
17:44:31  Starting VDJ barcode iteration 9 of 12...
Started mcapply cluster with 48 cores
17:45:10 Done with 9 of 12
17:45:10  Starting VDJ barcode iteration 10 of 12...
Started mcapply cluster with 48 cores
17:45:15 Done with 10 of 12
17:45:15  Starting VDJ barcode iteration 11 of 12...
Started mcapply cluster with 48 cores
17:46:42 Done with 11 of 12
17:46:42  Starting VDJ barcode iteration 12 of 12...
Started mcapply cluster with 48 cores
17:46:43 Done with 12 of 12

Preparing Seurat.in object
17:46:48 Done with Seurat.in

Integrating VDJ and GEX 
Error in VDJ_GEX_matrix(VDJ.out.directory.list = c(VDJ.out.directory.list_T,  : 
                                                     object 'batches' not found
                                                   In addition: Warning messages:
                                                     1: In data.frame(..., check.names = FALSE) :
                                                     row names were found from a short variable and have been discarded
                                                   2: In VDJ_GEX_matrix(VDJ.out.directory.list = c(VDJ.out.directory.list_T,  :
                                                                                                     Filtering of overlapping barcodes in GEX is not performed for Seurat.in input
                                                                                                   3: In SeuratObject::FetchData(GEX.proc, vars = c("orig.ident", "seurat_clusters",  :
                                                                                                                                                      The following requested variables were not found: tSNE_1, tSNE_2