UMAP using only a Subset of Loci from bed file?

[mfrenkel16]

Hi,

Rather than use reads distributed across the entire genome, I would like to subset the regions of the genome that are input into the dimension reduction before running the UMAP (i.e. I have regions of known high variability in a bed file and would like to subset the remainder of the analyses just to these regions). Is this possible? If I didn't give enough information to answer meaningfully, of course please let me know and I can try to clarify. It's also possible that my question is misplaced and maybe I should figure out how to do the subsetting for each cell barcode to only include those loci that overlap with the regions of the bed file at the level of the .tsv files before they're even input into ArchR?

Thanks!
Max

[rcorces]

create a PeakMatrix based off of your regions of interest using addPeakSet() and addPeakMatrix(). Then run addIterativeLSI() using that PeakMatrix as input. Then run addUMAP() etc on that new reducedDims object.

[mfrenkel16]

Thanks, Ryan! At the step of addInterativeLSI, I get a lengthy error which seems to boil down to the following (log files attached too):

ERROR Found in .LSIPartialMatrix...Error Found Iteration 1 :
[1] "Error in rownames<-(*tmp*, value = new_rownames) : \n invalid rownames length\n"
<simpleError in rownames<-(*tmp*, value = new_rownames): invalid rownames length>

I'm not sure how to resolve this. Are these rownames referring to those of the GRanges object with the peaks? (I notably did not get an error at the addPeakSet or addPeakMatrix steps). I made the GRanges object from a bed file, so I wonder if this means that I should reformat the original bed file (e.g. right now it's in the format of chr1, chr2...chrY but maybe ArchR requires a different alias for the chromosome names)? Maybe I'm entirely wrong about what the problem is though.

Thanks!
ArchR-addIterativeLSI-1466ba91c7a-Date-2022-11-03_Time-17-40-39.log
ArchR-addPeakMatrix-14644890000-Date-2022-11-03_Time-17-09-49.log

[rcorces]

@mfrenkel16 - apologies for the delay in replying. I dont have a good intuition for what is causing your error. If you're able to, providing a traceback() immediately after getting the error might help.

Since I'm unable to recapitulate the issue with the tutorial dataset, its hard for me to know what is happening on your end. My code below.

> projHeme5 <- loadArchRProject(path = "./Save-ProjHeme5/")
Successfully loaded ArchRProject!

                                                   / |
                                                 /    \
            .                                  /      |.
            \\\                              /        |.
              \\\                          /           `|.
                \\\                      /              |.
                  \                    /                |\
                  \\#####\           /                  ||
                ==###########>      /                   ||
                 \\##==......\    /                     ||
            ______ =       =|__ /__                     ||      \\\
        ,--' ,----`-,__ ___/'  --,-`-===================##========>
       \               '        ##_______ _____ ,--,__,=##,__   ///
        ,    __==    ___,-,__,--'#'  ==='      `-'    | ##,-/
        -,____,---'       \\####\\________________,--\\_##,/
           ___      .______        ______  __    __  .______      
          /   \     |   _  \      /      ||  |  |  | |   _  \     
         /  ^  \    |  |_)  |    |  ,----'|  |__|  | |  |_)  |    
        /  /_\  \   |      /     |  |     |   __   | |      /     
       /  _____  \  |  |\  \\___ |  `----.|  |  |  | |  |\  \\___.
      /__/     \__\ | _| `._____| \______||__|  |__| | _| `._____|
    
> ps <- getPeakSet(ArchRProj = projHeme5)
> ps
GRanges object with 143400 ranges and 12 metadata columns:
        seqnames              ranges strand |     score replicateScoreQuantile
           <Rle>           <IRanges>  <Rle> | <numeric>              <numeric>
   Mono     chr1       752502-753002      * |  15.35900                  0.723
      B     chr1       762685-763185      * |  48.88240                  0.933
  CD4.M     chr1       779917-780417      * |   7.02613                  0.335
    GMP     chr1       793283-793783      * |   3.48882                  0.384
      B     chr1       801006-801506      * |  14.66540                  0.697
    ...      ...                 ...    ... .       ...                    ...
   PreB     chrX 154840785-154841285      * |   3.21394                  0.664
  CD4.N     chrX 154841356-154841856      * |   1.35161                  0.282
   PreB     chrX 154842430-154842930      * |   5.32226                  0.866
     NK     chrX 154862053-154862553      * |  12.37540                  0.494
   Mono     chrX 154997960-154998460      * |   8.40370                  0.590
        groupScoreQuantile Reproducibility         GroupReplicate
                 <numeric>       <numeric>            <character>
   Mono              0.509               2  Mono._.scATAC_PBMC_R1
      B              0.845               2     B._.scATAC_PBMC_R1
  CD4.M              0.066               2 CD4.M._.scATAC_PBMC_R1
    GMP              0.098               2 GMP._.scATAC_CD34_BM..
      B              0.421               2     B._.scATAC_PBMC_R1
    ...                ...             ...                    ...
   PreB              0.278               2            PreB._.Rep2
  CD4.N              0.050               2 CD4.N._.scATAC_PBMC_R1
   PreB              0.599               2            PreB._.Rep2
     NK              0.200               2    NK._.scATAC_PBMC_R1
   Mono              0.308               2  Mono._.scATAC_PBMC_R1
        distToGeneStart nearestGene    peakType distToTSS  nearestTSS        GC
              <integer> <character> <character> <integer> <character> <numeric>
   Mono           10149   LINC00115      Distal     10149  uc001aau.3    0.4830
      B              32   LINC00115    Promoter        32  uc001aau.3    0.6926
  CD4.M           17195   LINC01128    Intronic     16988  uc021oeh.1    0.4890
    GMP           18648      FAM41C      Exonic     18648  uc021oei.1    0.3852
      B           10925      FAM41C      Distal     10925  uc021oei.1    0.4371
    ...             ...         ...         ...       ...         ...       ...
   PreB            1586       TMLHE    Intronic      1586  uc004cin.3    0.4850
  CD4.N            1015       TMLHE    Intronic      1015  uc004cin.3    0.3533
   PreB              57       TMLHE    Promoter        57  uc004cin.3    0.5749
     NK           19680       TMLHE      Distal     19680  uc004cin.3    0.4311
   Mono          155587       TMLHE      Distal       758  uc004cin.3    0.4591
              idx
        <integer>
   Mono         1
      B         2
  CD4.M         3
    GMP         4
      B         5
    ...       ...
   PreB      3413
  CD4.N      3414
   PreB      3415
     NK      3416
   Mono      3417
  -------
  seqinfo: 23 sequences from an unspecified genome; no seqlengths
> ps <- ps[sample(x = 1:length(ps), size = 50000, replace = FALSE)]
> ps2 <- GRanges(seqnames = seqnames(ps), ranges = ranges(ps), strand = strand(ps))
> names(ps2) <- NULL
> ps2
GRanges object with 50000 ranges and 0 metadata columns:
          seqnames              ranges strand
             <Rle>           <IRanges>  <Rle>
      [1]    chr11 131056398-131056898      *
      [2]     chr3   53079268-53079768      *
      [3]     chr7 151106748-151107248      *
      [4]    chr12   50101096-50101596      *
      [5]    chr15   78341921-78342421      *
      ...      ...                 ...    ...
  [49996]     chr3 158445458-158445958      *
  [49997]     chr6   26159463-26159963      *
  [49998]    chr12 131836078-131836578      *
  [49999]    chr15   65809825-65810325      *
  [50000]    chr12 118643548-118644048      *
  -------
  seqinfo: 23 sequences from an unspecified genome; no seqlengths
> projHeme5 <- addPeakSet(ArchRProj = projHeme5, peakSet = ps2, force = TRUE)
> projHeme5 <- addPeakMatrix(ArchRProj = projHeme5, force = TRUE)
ArchR logging to : ArchRLogs/ArchR-addPeakMatrix-3377cb71dd9929-Date-2022-11-07_Time-05-49-50.log
If there is an issue, please report to github with logFile!
2022-11-07 05:49:50 : Batch Execution w/ safelapply!, 0 mins elapsed.
Overriding previous entry for ReadsInPeaks
Overriding previous entry for FRIP
ArchR logging successful to : ArchRLogs/ArchR-addPeakMatrix-3377cb71dd9929-Date-2022-11-07_Time-05-49-50.log
> projHeme5 <- addIterativeLSI(ArchRProj = projHeme5, useMatrix = "PeakMatrix", name = "PeakIterativeLSI")
Checking Inputs...
ArchR logging to : ArchRLogs/ArchR-addIterativeLSI-3377cb69695639-Date-2022-11-07_Time-05-51-24.log
If there is an issue, please report to github with logFile!
2022-11-07 05:51:24 : Computing Total Across All Features, 0.004 mins elapsed.
2022-11-07 05:51:25 : Computing Top Features, 0.02 mins elapsed.
Not Enough Non-Zero Features to Filter!
###########
2022-11-07 05:51:27 : Running LSI (1 of 2) on Top Features, 0.052 mins elapsed.
###########
2022-11-07 05:51:27 : Sampling Cells (N = 10002) for Estimated LSI, 0.053 mins elapsed.
2022-11-07 05:51:27 : Creating Sampled Partial Matrix, 0.053 mins elapsed.
2022-11-07 05:51:31 : Computing Estimated LSI (projectAll = FALSE), 0.114 mins elapsed.
2022-11-07 05:51:50 : Identifying Clusters, 0.444 mins elapsed.
2022-11-07 05:52:05 : Identified 6 Clusters, 0.688 mins elapsed.
2022-11-07 05:52:05 : Saving LSI Iteration, 0.688 mins elapsed.
2022-11-07 05:52:19 : Creating Cluster Matrix on the total Group Features, 0.917 mins elapsed.
2022-11-07 05:52:26 : Computing Variable Features, 1.035 mins elapsed.
###########
2022-11-07 05:52:26 : Running LSI (2 of 2) on Variable Features, 1.035 mins elapsed.
###########
2022-11-07 05:52:26 : Creating Partial Matrix, 1.035 mins elapsed.
2022-11-07 05:52:30 : Computing LSI, 1.097 mins elapsed.
2022-11-07 05:52:51 : Finished Running IterativeLSI, 1.451 mins elapsed.

[mfrenkel16]

Hi Ryan,

Thanks for the code. I'm pretty sure I'm doing the same thing as you are. I've provided the traceback() immediately after the error in case that helps you diagnose what's happening.
addIterativeLSI_Traceback.txt

My granges object from my bed file of peaks looks to be formatted the same as the one you created except for the fact that mine is sorted but I doubt that matters (mine wasn't created by sampling from the ArchR project though which I also doubt is a problem):

>bed_grange
GRanges object with 1464 ranges and 0 metadata columns:
       seqnames              ranges strand
          <Rle>           <IRanges>  <Rle>
     1     chr1   11704939-11705441      *
     2     chr1   12049705-12050093      *
    3     chr1   12539885-12540806      *
     4     chr1   15795125-15795665      *
     5     chr1   15849265-15850249      *
   ...      ...                 ...    ...
  1460     chrX 136541059-136541750      *
  1461     chrX 140673019-140673488      *
  1462     chrX 149010069-149011094      *
  1463     chrX 149719930-149720640      *
  1464     chrX 154799391-154800472      *
------- 
  seqinfo: 23 sequences from an unspecified genome; no seqlengths

If there's anything else I could provide that would be helpful, please let me know. Otherwise, if it's still not clear what's happening, I could see if this works on the tutorial data or not.

Thanks!

[rcorces]

Theres a lot going on in this function and the traceback doesnt help pinpoint exactly where the error is occurring. If you could run the code I provided on the tutorial data and make sure that works, that would be a good place to start to make sure it isnt something with your environment that is causing problems.