Custom genome/gene annotations

[badoi]

Hey y'all,

I want to put this here for anyone who's using rhesus macaque, marmoset, and rat genomes. I'm putting together what I have and anyone in the ArchR community is welcome to use and contribute https://github.com/pfenninglab/custom_ArchR_genomes_and_annotations.

Thanks,
BaDoi

[rcorces]

@badoi - this is really great. thank you for putting this together. Is it ok if I link to this in an upcoming new release of the documentation?

[badoi]

Absolutely! Thank you for the “official” vote of confidence.

[sylestiel]

Hi @badoi

I'm using the genome/gene annotation as per your suggestion to run my Rhesus samples in ArchR. However, I'm encountering an error at the addGeneExpressionMatrix step.
Error in addGeneExpressionMatrix(input = proj, seRNA = seRNA, force = TRUE) :
Detected less than 50% of seqnames for seRNA in chromSizes! Are you sure the seqnames of your seRNA are correct?!
Can you help with this problem?

ArchR-addGeneExpressionMatrix-51bb61760554-Date-2023-12-19_Time-17-12-09.53781.log

[badoi]

Hi @sylestiel,

The way that the data is organized on-disk in .arrow files are chunked by chromosomes. The gene expression matrix is added to the ArchR matrix by chromosome as well b/c this is how the geneScoreMatrix is stored/accessed. To make the seRNA object work w/ the addGeneExpressionMatrix function, the seRNA object needs to have a rowRanges() metadata w/ genomic Ranges characterizing each gene. You can will need to find in the .gff or .gtf annotation file used to make the genome index for CellRanger or STARsolo (what I use), etc. Read this annotation file, and add the appropriate chromosome information for each gene into the rowRanges() for your seRNA object. For me, this is what it looks like. you'll probably need to do something different depending on how you got the seRNA part :).

Also, this is posted somewhere else in the geneExpressionMatrix() usage... it's intended to take raw, unprocessed UMI counts from a cellRanger project. Make sure to not add processed counts in the seRNA or you'll have doubly log-normalized expression data which will funk up downstream RNA-ATAC analyses.

## get the archR project and 
proj <- loadArchRProject(PROJDIR)

## read in gene annotations
annot = import('/home/bnphan/resources/genomes/rheMac10/rheMac10_liftoff_GRCh38.p13_RefSeq.gff3')
names(annot) = gsub('-[1-9]$', '', annot$Name)
names(annot) = gsub('_', '-', names(annot))

## load in raw counts of unlabeled multiomeRNA
sce_fn = file.path('data/tidy_data/rdas', 'raw_multiomeRNA_DLPFC_sce_list_N4.rds')
sce_list = readRDS(sce_fn)
sce = do.call('cbind', sce_list)

mcols(sce)$Symbol = gsub('gene-|-[1-9]$', '', mcols(sce)$Symbol)
sce = sce[!duplicated(mcols(sce)$Symbol)]
sce = sce[mcols(sce)$Symbol %in% names(annot), ]
rowRanges(sce) = annot[match(mcols(sce)$Symbol, names(annot))]

## will use the seRNA chromosome to parallelize stuffing geneExpression data into .arrows
proj <- addGeneExpressionMatrix(input = proj, seRNA = sce, force = TRUE)

## TADA ##