CINXE.COM
Unable to restore chromosome sequence lengths or other metadata in BSgenome.Mfascicularis.NCBI.6.0
<html> <head> <title> Unable to restore chromosome sequence lengths or other metadata in BSgenome.Mfascicularis.NCBI.6.0 </title> <meta charset="UTF-8"> <meta name="viewport" content="width=device-width, initial-scale=1"> <link rel="icon" href="/static/favicon.ico" type="image/x-icon"/> <script src="https://ajax.googleapis.com/ajax/libs/jquery/3.2.1/jquery.min.js"></script> <script async src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-chtml.js" id="MathJax-script"></script> <script src="https://ajax.googleapis.com/ajax/libs/jquery/3.2.1/jquery.min.js"></script> <script src="https://code.jquery.com/ui/1.12.1/jquery-ui.min.js"></script> <link disabled rel="stylesheet" href="//cdn.jsdelivr.net/gh/highlightjs/cdn-release@9.18.0/build/styles/default.min.css"> <script async src="//cdn.jsdelivr.net/gh/highlightjs/cdn-release@9.18.0/build/highlight.min.js"></script> <!-- Global site tag (gtag.js) - Google Analytics --> <script async src="https://www.googletagmanager.com/gtag/js?id=UA-55275703-1"></script> <script> window.dataLayer = window.dataLayer || []; function gtag() { dataLayer.push(arguments); } gtag('js', new Date()); gtag('config', "UA-55275703-1"); </script> <link rel="stylesheet" href="/static/CACHE/css/output.b3e21268cb4d.css" type="text/css"> <script src="/static/CACHE/js/output.8d9edd4e014b.js"></script> <link rel="stylesheet" href="/static/autocomplete/tribute.css" /> <link rel="stylesheet" href="/static/autocomplete/at.who.css" /> <link href="/static/pagedown.css" type="text/css" media="all" rel="stylesheet"> <link href="/static/prism.css" rel="stylesheet"> <link href="/static/pagedown/demo/browser/demo.css" type="text/css" media="all" rel="stylesheet"> <script src="/static/autocomplete/tribute.js"></script> <script type="text/javascript" src="/static/pagedown/Markdown.Converter.js"></script> <script type="text/javascript" src="/static/pagedown-extra/pagedown/Markdown.Converter.js"></script> <script type="text/javascript" src="/static/pagedown/Markdown.Sanitizer.js"></script> <script type="text/javascript" src="/static/pagedown/Markdown.Editor.js"></script> <script type="text/javascript" src="/static/pagedown-extra/Markdown.Extra.js"></script> <script type="text/javascript" src="/static/pagedown_init.js"></script> <script src="/static/prism.js"></script> <script src="/static/markdown-it.js"></script> <script src="/static/autocomplete/at.who.js"></script> </head> <body itemscope itemtype="https://schema.org/QAPage"> <div class="ui inverted container main"> <span class="menus"> <div class="ui top attached text menu" id="menu-topics"> <div class="ui inverted dimmer"></div> <div class="right menu" id="login-opts"> <a class="item " href="/accounts/login/"> <i class="sign-in icon"></i>Log In</a> <a class="item " href="/accounts/signup/"> <i class="user icon"></i> Sign Up</a> <a class="item " href="/info/about/"> <i class="info icon"></i>about</a> <a class="item " href="/info/faq/"> <i class="question icon"></i> faq </a> </div> </div> <div class="ui labeled icon attached pointing menu" id="menu-header"> <div class="header item" id="logo"> <a href="https://support.bioconductor.org/"><img class="ui image" src="/static/transparent-logo.png"></a> </div> <a class="item format add-question" href="/new/post/"> <span class="">Ask a question</span> </a> <a class="item format " href="/"> <span class="">Latest</span> </a> <a class=" item format" href="/t/news/"> <span class="">News</span> </a> <a class=" item format " href="/t/jobs/"> <span class="">Jobs</span> </a> <a class=" item format " href="/t/tutorials/"> <span class="">Tutorials</span> </a> <a class=" item format " href="/t/"> <span class="">Tags</span> </a> <a class=" item format " href="/user/list/"> <span class="">Users</span> </a> </div> </span> <span class="phone-menus"> <div class="header item" id="logo"> <a href="/"><img class="ui image" width="220px" src="/static/transparent-logo.png"></a> </div> <div class="ui labeled icon attached fluid pointing menu" id="menu-header" style="background-color: white"> <div class="ui left simple dropdown item"> <i class="bars icon"><i class="dropdown icon"></i></i> <div class="menu"> <a class="item " href="/new/post/"> <i class=" edit icon"></i>New Post </a> <a class="item " href="/"> <i class=" list icon"></i> Latest </a> <a class="item " href="/t/news/"> <i class=" newspaper icon"></i> News </a> <a class="item " href="/t/jobs/"> <i class=" briefcase icon"></i> Jobs </a> <a class="item " href="/t/tutorials/"> <i class=" info circle icon"></i> Tutorials </a> <a class="item " href="/t/"> <i class="tag icon"></i>Tags </a> <a class="item " href="/user/list/"> <i class=" world icon"></i> Users </a> </div> </div> <div class="ui right simple dropdown item"> </div> <div class="ui right simple dropdown item"> <i class="user icon"><i class="dropdown icon"></i></i> <div class="menu" id="login-opts"> <a class="item " href="/accounts/login/"><i class="sign-in icon"></i> Log In</a> <a class="item " href="/accounts/signup/"><i class="arrow circle up icon"></i> Sign Up</a> <a class="item " href="/info/about/"><i class="info circle icon"></i>About </a> </div> </div> </div> </span> <div class="ui bottom attached segment block"> <div class="ui stackable grid"> <div class="fit twelve wide content column"> <span itemprop="mainEntity" itemscope itemtype="https://schema.org/Question"> <div class="ui vertical segment"> <a name="9161293"></a> <div class="post view open" data-value="9161293"> <div class="title" itemprop="name">Unable to restore chromosome sequence lengths or other metadata in BSgenome.Mfascicularis.NCBI.6.0</div> <span itemprop="answerCount" style="display: none">0</span> <div class="body" > <div class="voting"> <button class="ui icon button" data-value="upvote" data-state="0"><i class="thumbs up icon"></i> </button> <div class="score" itemprop="upvoteCount">0</div> <button class="ui icon button" data-value="bookmark" data-state="0"><i class="bookmark icon"></i> </button> </div> <div class="content"> <div class="droppable inplace"> <div class="ui inverted dimmer"> <div class="ui text loader"> <div class="muted">Entering edit mode</div> </div> </div> <div class="wrap magnify"> <div class="user_box tablet right"> <div> <a class="" href="/u/79579/"> <img class="ui centered circular image" src="https://secure.gravatar.com/avatar/182054bebbc7802ad4798c4c06ee8e42?s=80&d=mp"> </a> </div> <div class="muted"> <div><a href="/u/79579/">Genevieve</a> • 0 </div> @fc87388f <div>Last seen 7 days ago</div> <div>United States</div> </div> </div> <span itemprop="text"> <p>Hello,</p> <p>I was using BSgenome.Mfascicularis.NCBI.6.0 genome for the <a href="https://stuartlab.org/signac/articles/pbmc_multiomic" rel="nofollow">Signac</a> multiomics tutorial. It was working fine until the I was trying to change the name of the chromosomes to match the GTF I have (macFas6-hg38-gencode.v47.basic.sortedgtf).</p> <p>This code worked before</p> <pre><code>mf_genome <- BSgenome.Mfascicularis.NCBI.6.0 # see 1st chromosomes head(seqlengths(mf_genome)) # see how many base pairs are in chromosome 1 mf_genome[["1"]] # chromosome 1 had 223606306 base pairs. # if you look at the chromsosomes table on bottom of page, chromsome 1 is 'CM021939.1' and it has 223606306 base pairs (no longer working???) # also lists the gc content % for each chromsome </code></pre> <p>But then when I tried to change the chromosome names...</p> <pre><code># # Store mapping in R as vector # chr_map <- c( # "1" = "CM021939.1", "2" = "CM021940.1", "3" = "CM021941.1", # "4" = "CM021942.1", "5" = "CM021943.1", "6" = "CM021944.1", # "7" = "CM021945.1", "8" = "CM021946.1", "9" = "CM021947.1", # "10" = "CM021948.1", "11" = "CM021949.1", "12" = "CM021950.1", # "13" = "CM021951.1", "14" = "CM021952.1", "15" = "CM021953.1", # "16" = "CM021954.1", "17" = "CM021955.1", "18" = "CM021956.1", # "19" = "CM021957.1", "20" = "CM021958.1", "X" = "CM021959.1" # ) # # # # Get the original chromosome names # original_chr_names <- seqnames(BSgenome.Mfascicularis.NCBI.6.0) # original_chr_names # # [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14" "15" "16" "17" "18" "19" "20" "X" # # # # # Check which names need to be changed (need to remove NA's from here --- rest of the scaffolds that aren't chromosome 1-20 or X) # new_chr_names <- chr_map[original_chr_names] # new_chr_names # # # 1 2 3 4 5 6 7 8 9 10 11 12 13 # # "CM021939.1" "CM021940.1" "CM021941.1" "CM021942.1" "CM021943.1" "CM021944.1" "CM021945.1" "CM021946.1" "CM021947.1" "CM021948.1" "CM021949.1" "CM021950.1" "CM021951.1" # # 14 15 16 17 18 19 20 X # # "CM021952.1" "CM021953.1" "CM021954.1" "CM021955.1" "CM021956.1" "CM021957.1" "CM021958.1" "CM021959.1" # # # # # # Keep only non-NA values ( chrom 1-20 & X) # # new_chr_names <- new_chr_names[!is.na(new_chr_names)] # # # # # Remove NA values to ensure lengths match # # valid_indices <- !is.na(new_chr_names) # # original_chr_names <- original_chr_names[valid_indices] # # new_chr_names <- new_chr_names[valid_indices] # # # # # Confirm lengths match # # length(original_chr_names) == length(new_chr_names) # # # # Assign the new chromosome names # names(BSgenome.Mfascicularis.NCBI.6.0@seqinfo) <- new_chr_names # # # Error in `seqnames<-`(x, value) : # # length of supplied 'seqnames' vector must equal the number of sequences # # # # # # Verify the update # seqnames(BSgenome.Mfascicularis.NCBI.6.0) 1 2 3 4 NA NA NA NA 5 6 7 8 NA NA NA NA 9 10 11 12 NA NA NA NA 13 14 15 16 NA NA NA NA 17 18 19 20 NA NA NA NA X Super-Scaffold_100058 Super-Scaffold_100060 Super-Scaffold_100064 NA NA NA NA Super-Scaffold_100065 Super-Scaffold_100066 Super-Scaffold_100068 Super-Scaffold_100069 NA NA NA NA Super-Scaffold_100072 Super-Scaffold_100073 Super-Scaffold_100078 Super-Scaffold_100080 NA NA NA NA </code></pre> <p>It removed all the other metadata. I've tried uninstalling and reinstalling BSgenome.Mfascicularis.NCBI.6.0 & BSgenome but it doesn't recover the seqlengths or the other metadata that was originally there.</p> <pre><code>remove.packages("BSgenome.Mfascicularis.NCBI.6.0") remove.packages("BSgenome") BiocManager::install("BSgenome") BiocManager::install("BSgenome.Mfascicularis.NCBI.6.0") </code></pre> <pre><code> print(names(BSgenome.Mfascicularis.NCBI.6.0@seqinfo)) [1] "CM021939.1" "CM021940.1" "CM021941.1" "CM021942.1" "CM021943.1" "CM021944.1" "CM021945.1" "CM021946.1" "CM021947.1" "CM021948.1" "CM021949.1" "CM021950.1" "CM021951.1" [14] "CM021952.1" "CM021953.1" "CM021954.1" "CM021955.1" "CM021956.1" "CM021957.1" "CM021958.1" "CM021959.1" </code></pre> <p>I tried this too <a href="https://support.bioconductor.org/p/83588/" rel="nofollow">linked here</a>. The metadata that was originally in BSgenome.Mfascicularis.NCBI.6.0 is no longer there and I can't restore it? I succesfully changed the chromosome names I needed and created a subsetted BS genome object called <code>new_genome</code> but again I no longer have seqlength info, etc</p> <pre><code># hack from https://support.bioconductor.org/p/83588/ to keep only some chromosomes within a BSgenome object. keepBSgenomeSequences <- function(genome, seqnames) { stopifnot(all(seqnames %in% seqnames(genome))) genome@user_seqnames <- setNames(seqnames, seqnames) genome@seqinfo <- genome@seqinfo[seqnames] genome } # load macfas6 NCBI BSgenome object library(BSgenome.Mfascicularis.NCBI.6.0) # assign to R object called genome genome <- BSgenome.Mfascicularis.NCBI.6.0 # check if seqlengths info is in original genome head(seqlengths(genome)) # 1 2 3 4 5 6 # NA NA NA NA NA NA # see how many base pairs are in chromosome 1 genome[["1"]] # Error in if (length(ans) != seqlengths(x)[[user_seqname]]) { : # missing value where TRUE/FALSE needed # # Store mapping in R as vector chr_map <- c( "1" = "CM021939.1", "2" = "CM021940.1", "3" = "CM021941.1", "4" = "CM021942.1", "5" = "CM021943.1", "6" = "CM021944.1", "7" = "CM021945.1", "8" = "CM021946.1", "9" = "CM021947.1", "10" = "CM021948.1", "11" = "CM021949.1", "12" = "CM021950.1", "13" = "CM021951.1", "14" = "CM021952.1", "15" = "CM021953.1", "16" = "CM021954.1", "17" = "CM021955.1", "18" = "CM021956.1", "19" = "CM021957.1", "20" = "CM021958.1", "X" = "CM021959.1" ) # create new R object called gew_genome with only chromosomes 1-20 & X # this is the new BSgenome object new_genome <- keepBSgenomeSequences(genome, names(chr_map)) seqnames(new_genome) <- chr_map # print out genome # | BSgenome object for Crab-eating macaque # | - organism: Macaca fascicularis # | - provider: NCBI # | - genome: Macaca_fascicularis_6.0 # | - release date: 2020/03/10 # | - 936 sequence(s): # | 1 2 3 4 # | 5 6 7 8 # | 9 10 11 12 # | 13 14 15 16 # | 17 18 19 20 # | ... ... ... ... # | tig00001423_obj tig00001424_obj tig00001425_obj tig00001428_obj # | tig00001430_obj tig00001431_obj tig00001432_obj tig00001433_obj # | tig00001434_obj tig00001435_obj tig00001436_obj tig00001437_obj # | tig00001438_obj tig00001440_obj tig00001441_obj tig00001442_obj # | tig00001443_obj tig00001444_obj tig00001445_obj tig00001446_obj # | # | Tips: call 'seqnames()' on the object to get all the sequence names, call 'seqinfo()' to get the full sequence info, use the '$' or '[[' operator to access a given sequence, see # | '?BSgenome' for more information. new_genome # | BSgenome object for Crab-eating macaque # | - organism: Macaca fascicularis # | - provider: NCBI # | - genome: Macaca_fascicularis_6.0 # | - release date: 2020/03/10 # | - 21 sequence(s): # | CM021939.1 CM021940.1 CM021941.1 CM021942.1 CM021943.1 CM021944.1 CM021945.1 CM021946.1 CM021947.1 CM021948.1 CM021949.1 CM021950.1 CM021951.1 CM021952.1 CM021953.1 CM021954.1 # | CM021955.1 CM021956.1 CM021957.1 CM021958.1 CM021959.1 # | # | Tips: call 'seqnames()' on the object to get all the sequence names, call 'seqinfo()' to get the full sequence info, use the '$' or '[[' operator to access a given sequence, see # | '?BSgenome' for more information. chr_map # 1 2 3 4 5 6 7 8 9 10 11 12 13 14 # "CM021939.1" "CM021940.1" "CM021941.1" "CM021942.1" "CM021943.1" "CM021944.1" "CM021945.1" "CM021946.1" "CM021947.1" "CM021948.1" "CM021949.1" "CM021950.1" "CM021951.1" "CM021952.1" # 15 16 17 18 19 20 X # "CM021953.1" "CM021954.1" "CM021955.1" "CM021956.1" "CM021957.1" "CM021958.1" "CM021959.1" seqinfo(new_genome) # Seqinfo object with 21 sequences from an unspecified genome; no seqlengths: # seqnames seqlengths isCircular genome # CM021939.1 <NA> <NA> <NA> # CM021940.1 <NA> <NA> <NA> # CM021941.1 <NA> <NA> <NA> # CM021942.1 <NA> <NA> <NA> # CM021943.1 <NA> <NA> <NA> # ... ... ... ... # CM021955.1 <NA> <NA> <NA> # CM021956.1 <NA> <NA> <NA> # CM021957.1 <NA> <NA> <NA> # CM021958.1 <NA> <NA> <NA> # CM021959.1 <NA> <NA> <NA> # see 1st chromosomes head(seqlengths(new_genome)) # CM021939.1 CM021940.1 CM021941.1 CM021942.1 CM021943.1 CM021944.1 # NA NA NA NA NA NA # see how many base pairs are in chromosome 1 new_genome[["CM021939.1"]] # Error in if (length(ans) != seqlengths(x)[[user_seqname]]) { : # missing value where TRUE/FALSE needed </code></pre> <p>I tried manually adding the sequence length info in the seqlengths column of the <code>new_genome</code> BS genome object using chromosomes info I retrieved from NCBI for Macaca_fascicularis_6.0 <a href="https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_011100615.1/" rel="nofollow">here</a>, but I still wan't able to add tha to the BS genome object.</p> <pre><code># Read the TSV file (assuming tab-separated and columns: seqname, length, gc_content or gc_count) seq_info <- read.table("~/macaque_multiomics/macfas6_ncbi_sequence_report_chromosomes.tsv", header = TRUE, sep = "\t", stringsAsFactors = FALSE) # Check column names to ensure correct assignment head(seq_info) # Print column names print(colnames(seq_info)) [1] "Assembly.Accession" "Assembly.unit.accession" "Chromosome.name" "GC.Count" "GC.Percent" "GenBank.seq.accession" [7] "Molecule.type" "Ordering" "RefSeq.seq.accession" "Role" "Seq.length" "UCSC.style.name" [13] "Unlocalized.Count" "Sequence.name" # Extract chromosome names from BSgenome object bs_chroms <- seqnames(new_genome) # Subset seq_info to only include relevant rows seq_info_subset <- seq_info[seq_info$GenBank.seq.accession %in% bs_chroms, ] # Check structure str(seq_info_subset) head(seq_info_subset[, c("GenBank.seq.accession", "Seq.length")]) #tried assigning sequence lenghts to new_genome # Convert column types seq_info_subset$GenBank.seq.accession <- as.character(seq_info_subset$GenBank.seq.accession) seq_info_subset$Seq.length <- as.numeric(seq_info_subset$Seq.length) # Create a named vector for sequence lengths seq_lengths <- setNames(seq_info_subset$Seq.length, seq_info_subset$GenBank.seq.accession) # Assign sequence lengths to BSgenome object seqlengths(new_genome) <- seq_lengths # Error in .check_new2old_and_new_seqinfo(new2old, value, seqinfo(x), context) : # seqlengths() and isCircular() of the supplied 'seqinfo' must be identical to seqlengths() and isCircular() of the current 'seqinfo' when replacing the 'seqinfo' of # a BSgenome object </code></pre> <p>I kept getting the same error when trying many different ways to add the sequence lengths to the BS genome object. Is it true that once a BSgenome object is created, you cannot directly modify the seqinfo slot, including the seqlengths? I don't know why the original BS genome object BSgenome.Mfascicularis.NCBI.6.0 doesn't have any metadata besides chromosome names even after unloading & reloading object. I've also tried uninstalling and reinstalling BSgenome.Mfascicularis.NCBI.6.0 and the BSgenome R packages but no luck.</p> </span> </div> <span class="inplace-tags"> <a class="ptag" href="/tag/BSgenome.Mfascicularis.NCBI.6.0/"> BSgenome.Mfascicularis.NCBI.6.0 </a> • 47 views </span> <time itemprop="dateCreated" datetime="2025-02-14"></time> <div class="actions top-level"> <a class="add-comment ui tiny label">ADD COMMENT</a> • <a itemprop="url" href="/p/9161293/#9161293">link</a> <span class="status muted user-info"> <span> 11 days ago <a href="/u/79579/"> Genevieve </a> • 0 </span> </span> </div> <span class="diff-cont"></span> </div> </div> </div> </div> </div> </span> <div class="ui info message"> <p><a class="ui small label" href="/accounts/login/"> <i class="sign in icon"></i>Login</a> before adding your answer.</p> </div> <script> $(document).ready(function () { $('#subscribe').dropdown(); drag_and_drop(); $('.ui.dropdown').dropdown(); //var users = "**MISSING**".split(','); autocomplete_users(); //init_pagedown(); $(this).on('click', '#inplace .save', function () { event.preventDefault(); var post = $(this).closest('.post'); edit_dropdown_post(post); }); $(this).on('click', '#inplace .create', function () { event.preventDefault(); create_comment(); }); $(this).on('click', '.edit-button', function (event) { event.preventDefault(); inplace_form($(this)); }); $(this).on('click', ".add-comment", function () { inplace_form($(this), true); }); }); </script> </div> <div class="four wide column sidefeed"> <div class="ui large wrap-text" id="sidebar"> <div class="ui large wrap-text" id="sidebar"> <div class="title">Similar Posts</div> <div class="ui basic segment similar" id="similar"> <div id="similar-feed" post_uid="9161293"></div> <div class="ui inverted dimmer"> <div class="content"> <div class="ui text loader"> <div class="muted">Loading Similar Posts</div> </div> </div> </div> </div> </div> </div> </div> </div> <div id="traffic">Traffic: 601 users visited in the last hour</div> </div> <div class="ui three center aligned column stackable tiny grid"> <div class=" left aligned column" style="padding-right: 13%"> <b>Content</b> <a href="/#search-anchor">Search</a><br> <a href="/user/list/">Users</a><br> <a href="/t/">Tags</a><br> <a href="/b/list/">Badges</a> </div> <div class="left aligned column" style="padding-right: 12%"> <b>Help</b> <a href="/info/about/">About</a><br> <a href="/info/faq/">FAQ</a><br> </div> <div class=" left aligned column"> <b>Access</b> <a href="/info/rss/">RSS</a><br> <a href="/info/api/">API</a><br> <a href="#">Stats</a> </div> </div> <div class="ui divider"></div> <div class="ui vertical center aligned basic segment"> <p>Use of this site constitutes acceptance of our <a href="/info/policy/">User Agreement and Privacy Policy</a>.</p> <div class="smaller muted"> Powered by the <a href="https://github.com/ialbert/biostar-central" class="ui image"> <img src="/static/images/badge-forum.svg"></a> version 2.3.6 </div> </div> </div> </body> </html>