Create genome pkg for Lotus japonicus from NCBI

[francisco-franco-86]

Hello,

I'm trying to create a genome package for Lotus japonicus from NCBI and I'm getting the following error:

forgeBSgenomeDataPkgFromNCBI(assembly_accession="GCA_012489685.2",
pkg_maintainer="Francisco Franco <ffrancopiccioli@gmail.com>",
organism="Lotus japonicus",
circ_seqs=c("MT"))

The file GCA_012489685.2_LjGifu_v1.2_genomic.fna.gz is already in the current directory, so it will be used.

Error in .sort_and_rename_fasta_sequences(dna, assembly_accession):

Incompatible assembly report and FASTA file for assembly GCA_012489685.2: the number of sequences in the FASTA file is not equal to the number of sequences in 'getChromInfoFromNCBI("GCA_012489685.2")' that have a non-NA GenBank accession.

I think there's an inconsistency between the ChromInfo obtained from NCBI (getChromInfoFromNCBI("GCA_012489685.2")) and the downloaded FASTA file, but I don't know how to fix it.

[https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_012489685.1/]

I hope you can help me create the package.

[hpages]

There's an inconsistency between the assembly report for assembly GCA_012489685.2 and the corresponding FASTA file available at NCBI: the latter contains the chloroplast sequence even though that sequence is not included in the assembly report. This inconsistency confuses forgeBSgenomeDataPkgFromNCBI().

If you're ok to not include the chloroplast sequence in your BSgenome data package, you can remove it from the FASTA file with (assuming you already have the file locally):

library(Biostrings)
dna <- readDNAStringSet("GCA_012489685.2_LjGifu_v1.2_genomic.fna.gz")
names(dna)[67]
# [1] "AP022636.1 Lotus japonicus B-129 chloroplast DNA, complete genome"
writeXStringSet(dna[-67], "GCA_012489685.2_LjGifu_v1.2_genomic.fna.gz", compress=TRUE)

Then forge the BSgenome data package with:

forgeBSgenomeDataPkgFromNCBI("GCA_012489685.2",
    pkg_maintainer="Francisco Franco <ffrancopiccioli@gmail.com>",
    organism="Lotus japonicus", circ_seqs="MT")
# File GCA_012489685.2_LjGifu_v1.2_genomic.fna.gz is already in current directory so will be used.
# Creating package in ./BSgenome.Ljaponicus.NCBI.LjGifuv1.2 

Hope this helps,
H.

[francisco-franco-86]

Hi Hervé Pagès,

Thank you for your quick response.

With this information I can now work, but I would like to know how to include the chloroplast sequence in the BSgenome data package.

Best regards,
Francisco

[hpages]

With this information I can now work, but I would like to know how to include the chloroplast sequence in the BSgenome data package.

You can use forgeBSgenomeDataPkgFromTwobitFile() for that. But first you'll need to download the original FASTA file containing all the genomic sequences, load the sequences, rename/reorder them, and save them in 2bit format.

Here's how to proceed:

library(BSgenomeForge)
assembly_accession <- "GCA_012489685.2"

## Download the FASTA file.
downloadGenomicSequencesFromNCBI(assembly_accession)

## Load the DNA seqences.
dna <- readDNAStringSet("GCA_012489685.2_LjGifu_v1.2_genomic.fna.gz")
names(dna) <- sub(" .*$", "", names(dna))  # clean the names

## Rename/reorder the DNA seqences.
chrominfo <- getChromInfoFromNCBI(assembly_accession)
seqnames <- chrominfo[ , "SequenceName"]
accns <- chrominfo[ , "GenBankAccn"]
seqnames <- append(seqnames, "chloroplast", after=7)
accns <- append(accns, "AP022636.1", after=7)
dna <- setNames(dna[accns], seqnames)

## Save the DNA sequences in 2bit format.
rtracklayer::export.2bit(dna, "GCA_012489685.2_LjGifu_v1.2_genomic.2bit")

## Call `forgeBSgenomeDataPkgFromTwobitFile()` on the 2bit file.
forgeBSgenomeDataPkgFromTwobitFile("GCA_012489685.2_LjGifu_v1.2_genomic.2bit",
    organism="Lotus japonicus", provider="NCBI", genome="LjGifu_v1.2",
    pkg_maintainer="Francisco Franco <ffrancopiccioli@gmail.com>",
    circ_seqs=c("MT", "chloroplast"))

H.

[hpages]

@francisco-franco-86 Did this work for you? Can we close this issue?