-
Notifications
You must be signed in to change notification settings - Fork 8
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Forge BSgenome data package for UCSC genome canFam6 #46
Conversation
Thanks @Priceless-P Can I ask you to rename the Also there's no need to put GenomicFeatures in Does the BSgenome.Cfamiliaris.UCSC.canFam6 package that you forged contain examples? Start R, load the package, and open the only man page in the package with Note that in more recent versions of the seed files, we've renamed the Here is an interesting experiment. Do
This is no coincidence. When you forged the package with By the way, specifying Another even more interesting experiment is the following. Right now the sequences in the
However, the sequences can be renamed to use the NCBI names. Try:
This works because of all the work you've done with registering You should be able to switch back to the UCSC names with:
The ability to switch easily from one naming convention to the other is a very popular feature in Bioconductor. There has always been a strong demand for it. Finally note that all the above is not restricted to BSgenome objects. It also applies to other Bioconductor objects, like TxDb objects:
|
Oh looks like you renamed |
Lol. I wish I can. Actually, I just happened to be on the repo when I saw your update to the text file in the extdata/Outreachy folder telling us to use the package name with -seed attached to it as the name of the seed file, so I did. Didn't even know you were writing a comment for me.
It's beautiful how these functions are linked behind the scene and how you configured them to work together. When I took on this issue, I thought it would take so much to complete. But I was marveled to discover that just one function
I didn't realize
What might the use cases be? Why prefer one to another when the information provided by both UCSC and NCBI are publicly available?
I came across I learn so much every single time you drop a comment for me. It's a huge privilege to have this much access to you. I will make the most out of it. Thank you for all you do @hpages |
Hi @Priceless-P
It's great that you felt that way 😊 However, note that this is an ideal situation, because the genome is already registered in GenomeInfoDb. Also the fact that we're dealing with an UCSC genome helps (those genomes are the easiest to deal with). Forging a BSgenome data package can get a a little bit more complicated when the genome is not registered, or if we're dealing with a non-UCSC genome (e.g. an NCBI assembly).
Of course it's allowed. It's even encouraged! We have nothing to hide 😉 This is open source. That's why open source is so great! You'll find the codebase for
The last line also shows in what package the function is defined. For
This tells us that S4 is a class system in R that allows OO programming (Object Oriented programming). In this class system, everything is about S4 classes, S4 generic functions, and S4 methods. We use those things a lot in Bioconductor. FWIW the S4Vectors package contains some slides that give a quick overview of the S4 class system. A link to those slides is available on the landing page of the S4Vectors package here. So when we display the body of We can display the list of methods defined for that generic function with:
Only 3 methods. One defined for BSgenomeDataPkgSeed objects, one for character objects (a.k.a. character vectors), and one for list objects. To see a particular method, we can use
This shows us the body of the method, and in what package the method is defined. There's a lot going on here. S4 is a rich and powerful class system. But it has some important differences with more "traditional" OO programming in other languages like Python, Java, or C++, that can make it a little bit harder to learn. It also has some unique (and very useful) features like S4 coercion methods for turning objects of a given class into objects of another class. There's a lot to discuss about it, but I will stop here for now. This comment is starting to become quite long so I'll take a break and post it. But I'm not done commenting and answering your questions. I'll do this in the next comment. |
Okay. Thanks for the clarification. I definitely have a lot of studying to do. How about my PR? Did I do something wrong? |
Sorry for the wait @Priceless-P, you did nothing wrong, I just got busy with other duties. Have you tried to remove the If you can confirm that having About this field:
I wonder if this is really the file that you downloaded. The extension in
They don't say this is a FASTA file, but we know this because of the If we look at the directory listing:
we see that Let's take a closer look at the file. Note that we can uncompress it with
The file is big so the command will take some time to uncompress it. When it's done,
Not surprisingly, the file is even bigger now. About 3x bigger! We didn't really need to uncompress the file. The only reason we did it is that now we can take a look at what's inside the file:
This is a typical FASTA file. It's organized in records. Each record starts with a header line that contains the name of a sequence, immediately followed by the sequence itself. The first sequence in our file is the DNA sequence of the first Dog chromosome, named Let's extract all the header lines from the file:
We should see one record per chromosome in the Let's count them:
That's exactly what we expected. So all seems to be good with this file: it contains all the sequences of the
In other words, the forging process only supports the two following options for the sequence data: a single twoBit file (like I wonder how you handled this @Priceless-P. Did you manage to split the big |
Okay. I understand. Yes, I removed
Yes, that is the file I downloaded. I unzipped it the same way you did above and it. I spilt and added them to a folder with this command I found on the Bioconductor support page
It was used for another genome but I did the necessary substitution. I also studied the code and I now understand how it works. |
Awesome! It's amazing that you figured how to do this. Here are some notes about this code:
So here's what I'll ask you to try:
As you can see, I'm using this task as an opportunity to introduce some key concepts and discuss a few other important things with you. There will be more. Hope you don't mind. (And I will answer your question about |
Okay, @hpages . I will redo it like you said.
Surprisingly
Also,
I think they are indeed identical. |
This code cannot produce 147 files. Maybe you've tried several other things while you were trying to figure out how to split the big FASTA file, and you probably managed to produce 147 files at some point. But that When I run this code on my laptop, I end up with the following files in
Then when I try to run
|
You are absolutely right. Indeed I tried several other things while trying to figure it out, I'm sorry for the mix-up.
All checked. And yes, I'm sorry for the stress🥺 |
No stress
Excellent! Make sure to also try all the other things discussed earlier like Earlier you asked what might the use cases be for switching between UCSC and NCBI names. This is a great question, thanks for asking! The typical use case is to harmonize the sequence names between two datasets. There's often the need to do this when a bioinformatician is in possession of two datasets that are associated with the same genome assembly, but one dataset uses the UCSC names and the other one the NCBI names. For example, suppose that a bioinformatician has forged and installed BSgenome.Cfamiliaris.UCSC.canFam6, but a coworker sent them a GFF file that uses the NCBI sequence names. GFF files are files that contain annotations for a given genome assembly, like gene positions and other genomic features. See here and here to learn more about these files. More precisely, let's say that a coworker decided that they want to work with the
The resulting object is a TxDb object:
Let's take a closer look at this object:
We know that the GFF file is associated with the
The names in the
But it doesn't matter because we can switch them to the official NCBI or UCSC sequence names with either Now, as part of their analysis, our bioinformatician needs to obtain the gene sequences for this assembly. The standard way to do this is to first obtain the gene positions (sometimes called "gene ranges"):
Then, load the BSgenome data package that contains the DNA sequences of this assembly and use
The reason this doesn't work is because the sequence names in
Note that this harmonization can also be done the other way around, that is, by doing Bottom line: The need to harmonize the sequence names between two datasets is a very common situation. And registering and linking a UCSC/NCBI pair like you did for
The Back to the process of forging BSgenome.Cfamiliaris.UCSC.canFam6: Have you considered using
The 2bit file format is a modern alternative to the FASTA format for storing DNA sequences. 2bit files tend to be smaller than their FASTA counterpart, and they allow much faster random access. As mentioned in the "How to forge a BSgenome data package" vignette, the 2bit file provided by UCSC can be used to forge a BSgenome data package. It actually makes the process easier because, unlike with the big FASTA file, the 2bit file doesn't need to be split into a collection of smaller files (one small file per sequence). Using a 2bit file instead of a FASTA file, whenever possible, is actually the preferred way to go. |
Yes, I tried them. They all worked as expected. I used
This makes a lot of sense. Thank you for explaining it to me.
You mentioned it in this vignette, and actually it was the first file I downloaded. But somehow i was still expecting to make a file per chromosome from the 2bit. Apparently i didn't really understand how it worked. I also got Thank you for opening my eyes to this. |
🚀
Yes, what gets included in the BSgenome data package is a 2bit file, independently of whether the input sequence data was in a collection of FASTA files or in a 2bit file. This file is produced by Note that when the input sequence data is in a 2bit file, the Would you mind updating your seed file to use Thanks! |
Yes, gladly! |
Thank you. About the The Some background about how the The Click on the Click on the As you can see, the
In fact, those values (e.g. Let's go one level up and click on the Click on the Every man page should have an
So when we call
Note that the Back to the So that is what the
Which becomes this line after injection:
Finally note that:
|
Thank you very much for the detailed explanation. It's very helpful. |
I added the canFam6 seed file.