Error when adding Kmer Bias to Coverage Files in addGroupCoverages

[e75967259]

Hi, everyone:

I am using the custom genome annotation of the T2T.CHM13v2.0.
How I create the annotation files is as follows:

chrs<-c("chr1", "chr2", "chr3", "chr4", "chr5", "chr6", "chr7", "chr8", "chr9", "chr10", "chr11", "chr12", "chr13", "chr14", "chr15", "chr16", "chr17", "chr18", "chr19", "chr20", "chr21", "chr22", "chrX", "chrY","chrMT")

genome <- BSgenome.Hsapiens.NCBI.T2T.CHM13v2.0
genome@seqinfo@seqnames<-chrs

genomeAnnotation <- createGenomeAnnotation(genome = genome)

As you can see I change the seqnames to fit them to my reference file.

The analysis process goes well until I use the addGroupCoverages() to make the pseudo-bulk replicates.

I got an error "chr1 Error in BSgenome[[availableChr[x]]] : no such sequence" when it went to add the kmer bias.

I notice that I may have the same problem as #322 , but the author didn't provide his solution there.

I think that the command didn't use the genome I created but use the seqname from the original "BSgenome.Hsapiens.NCBI.T2T.CHM13v2.0", so it can't match to the chromosome name in my project file.

I want to ask can I let the addGroupCoverages() use the custom annotation file to add kmer bias?

Here is my log file:
ArchR-addGroupCoverages-618b45f340f2-Date-2023-05-14_Time-15-16-12.log

[rcorces]

As far as I know, this:

chrs<-c("chr1", "chr2", "chr3", "chr4", "chr5", "chr6", "chr7", "chr8", "chr9", "chr10", "chr11", "chr12", "chr13", "chr14", "chr15", "chr16", "chr17", "chr18", "chr19", "chr20", "chr21", "chr22", "chrX", "chrY","chrMT")
genome <- BSgenome.Hsapiens.NCBI.T2T.CHM13v2.0
genome@seqinfo@seqnames<-chrs

doesnt do anything to the actual underlying sequences?

[e75967259]

This change is necessary because the geneAnnotation file uses the chromosome label like "chr1" not "1" which is the label used in the BSgenome.Hsapiens.NCBI.T2T.CHM13v2.0.

If I don't do anything I will get an error when I create the arrow file.

Sorry, I think I should show the following steps I took.

I use the "genome" which's seqnames has been changed to create the annotation object by this:

genomeAnnotation .T2T <- createGenomeAnnotation(genome = genome)

and use this custom genomeAnnotation object to create ArchR arrow file (the genomeAnnotation.T2T is customers too):

ArrowFiles <- createArrowFiles(
  geneAnnotation = geneAnnotation.T2T,
  genomeAnnotation = genomeAnnotation.T2T,
  inputFiles = inputFiles,
  sampleNames = sampleNames,
  minTSS = 4, #Dont set this too high because you can always increase later
  minFrags = 100, 
  addTileMat = TRUE,
  addGeneScoreMat = TRUE
)

[rcorces]

As I said before, I dont think that you are actually changing the seqnames of the underlying object, you're just changing some metadata as shown by the below code

> genome
Human genome:
# organism: Homo sapiens (Human)
# genome: hg38
# provider: UCSC
# release date: Feb 2019
# 640 sequences:
#   chr1                    chr2                    chr3                    chr4                    chr5                   
#   chr6                    chr7                    chr8                    chr9                    chr10                  
#   chr11                   chr12                   chr13                   chr14                   chr15                  
#   chr16                   chr17                   chr18                   chr19                   chr20                  
#   chr21                   chr22                   chrX                    chrY                    chrM                   
#   ...                     ...                     ...                     ...                     ...                    
#   chr17_KV766198v1_alt    chr17_KZ559114v1_alt    chr18_KQ458385v1_alt    chr18_KZ559116v1_alt    chr19_KV575246v1_alt   
#   chr19_KV575247v1_alt    chr19_KV575248v1_alt    chr19_KV575249v1_alt    chr19_KV575250v1_alt    chr19_KV575251v1_alt   
#   chr19_KV575252v1_alt    chr19_KV575253v1_alt    chr19_KV575254v1_alt    chr19_KV575255v1_alt    chr19_KV575256v1_alt   
#   chr19_KV575257v1_alt    chr19_KV575258v1_alt    chr19_KV575259v1_alt    chr19_KV575260v1_alt    chr22_KN196485v1_alt   
#   chr22_KN196486v1_alt    chr22_KQ458387v1_alt    chr22_KQ458388v1_alt    chr22_KQ759761v1_alt    chrX_KV766199v1_alt    
# (use 'seqnames()' to see all the sequence names, use the '$' or '[[' operator to access a given sequence)
> temp <- genome@seqinfo@seqnames
> temp[1]
[1] "chr1"
> temp[1] <- "chr1_temp"
> genome@seqinfo@seqnames <- temp
> genome
Human genome:
# organism: Homo sapiens (Human)
# genome: hg38
# provider: UCSC
# release date: Feb 2019
# 640 sequences:
#   chr1                    chr2                    chr3                    chr4                    chr5                   
#   chr6                    chr7                    chr8                    chr9                    chr10                  
#   chr11                   chr12                   chr13                   chr14                   chr15                  
#   chr16                   chr17                   chr18                   chr19                   chr20                  
#   chr21                   chr22                   chrX                    chrY                    chrM                   
#   ...                     ...                     ...                     ...                     ...                    
#   chr17_KV766198v1_alt    chr17_KZ559114v1_alt    chr18_KQ458385v1_alt    chr18_KZ559116v1_alt    chr19_KV575246v1_alt   
#   chr19_KV575247v1_alt    chr19_KV575248v1_alt    chr19_KV575249v1_alt    chr19_KV575250v1_alt    chr19_KV575251v1_alt   
#   chr19_KV575252v1_alt    chr19_KV575253v1_alt    chr19_KV575254v1_alt    chr19_KV575255v1_alt    chr19_KV575256v1_alt   
#   chr19_KV575257v1_alt    chr19_KV575258v1_alt    chr19_KV575259v1_alt    chr19_KV575260v1_alt    chr22_KN196485v1_alt   
#   chr22_KN196486v1_alt    chr22_KQ458387v1_alt    chr22_KQ458388v1_alt    chr22_KQ759761v1_alt    chrX_KV766199v1_alt    
# (use 'seqnames()' to see all the sequence names, use the '$' or '[[' operator to access a given sequence)
> 

[e75967259]

Sorry, I miss a line of code.

genome <- BSgenome.Hsapiens.NCBI.T2T.CHM13v2.0
genome@seqinfo@seqnames<-chrs
genome@user_seqnames<-chrs

This is my output:

genome <- BSgenome.Hsapiens.NCBI.T2T.CHM13v2.0
genome
Human genome:
# organism: Homo sapiens (Human)
# genome: T2T-CHM13v2.0
# provider: NCBI
# release date: 2022/01/24
# 25 sequences:
#   1   2   3   4   5   6   7   8   9   10  11  12  13  14  15  16  17  18  19  20  21  22  X   Y   MT     
# (use 'seqnames()' to see all the sequence names, use the '$' or '[[' operator to access a given sequence)


chrs<-c("chr1", "chr2", "chr3", "chr4", "chr5", "chr6", "chr7", "chr8", "chr9", "chr10", "chr11", "chr12", "chr13", "chr14", "chr15", "chr16", "chr17", "chr18", "chr19", "chr20", "chr21", "chr22", "chrX", "chrY","chrMT")
chrs.fac<-factor(chrs, levels = chrs)

genome@seqinfo@seqnames<-chrs
genome@user_seqnames<-chrs

genome
Human genome:
# organism: Homo sapiens (Human)
# genome: T2T-CHM13v2.0
# provider: NCBI
# release date: 2022/01/24
# 25 sequences:
#   chr1  chr2  chr3  chr4  chr5  chr6  chr7  chr8  chr9  chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17
#   chr18 chr19 chr20 chr21 chr22 chrX  chrY  chrMT                                                      
# (use 'seqnames()' to see all the sequence names, use the '$' or '[[' operator to access a given sequence)

[rcorces]

I dont have enough information to help. Your error is saying that your chromsomes arent present in your arrow files.
what does .availableSeqnames(getArrowFiles(ArchRProj)) show?

[e75967259]

This is my output:

.availableSeqnames(getArrowFiles(proj))
[1] "chr1"  "chr10" "chr11" "chr12" "chr13" "chr14" "chr15" "chr16" "chr17" "chr18" "chr19" "chr2" 
[13] "chr20" "chr21" "chr22" "chr3"  "chr4"  "chr5"  "chr6"  "chr7"  "chr8"  "chr9"  "chrX"  "chrY"

[rcorces]

I think that the command didn't use the genome I created but use the seqname from the original "BSgenome.Hsapiens.NCBI.T2T.CHM13v2.0", so it can't match to the chromosome name in my project file.

This is definitely not the case. The BSgenome object is obtained from the project:

ArchR/R/GroupCoverages.R

Lines 248 to 255 in c61b064

	o <- .addKmerBiasToCoverage(
	coverageMetadata = coverageMetadata,
	genome = getGenome(ArchRProj),
	kmerLength = kmerLength,
	threads = threads,
	verbose = FALSE,
	logFile = logFile
	)

You'll have to step through the code to figure this out on your end. There isnt much more I can do. This is where the availableChr is coming from
https://github.com/GreenleafLab/ArchR/blob/c61b0645d1482f80dcc24e25fbd915128c1b2500/R/GroupCoverages.R#LL608C4-L608C64

[e75967259]

Thanks for your help!!

I'll try what I can to.