I was working on a material that I put together as an example of studying the Ebola data and strangely code that ran fine in the Fall now did not produce any results. Quite the head scratcher ...
I ended up troubleshooting it for a while only to realize that most of the new data deposited into the main Ebola project does not actually map to Ebola. At all. In fact 541 files map at rates below 1%!
Below is the code used to evaluate the mapping percentages for all data (891 files) in this project, see the results at the end). Hopefully it might save some time to someone else.
0.25 bam/SRR1735115.bam 100 + 0 properly paired (0.25%:nan%)
I have repeated the mapping for a read-set you have found only 100 reads for. I mapped SRR1735115 to KR817241.1 with bwa mem. I found 7690 reads mapping. This is a 40-fold coverage in average (the Ebola genome is only 18000 nt in size). I inspected the mapping with Tablet. The reads show a very high error rate. I have never inspected RNA sequencing experiments before. Maybe that is normal? Nevertheless, I could clearly identify a few SNP with respect to the reference genome (which is from the 2014 outbreak).
In my opinion the challenge is sample preparation (as always). Ebola is a RNA virus. You have to prepare RNA from human blood serum, then reverse-transcribe it into DNA. Isolating RNA from blood is tricky and error-prone.
You stated that all samples were positive in qPCR (quantitative PCR). PCR is very sensitive. If PCR is positive, this does not mean that there is enough RNA in the serum for sequencing. With qPCR you can determine the viral load. The viral load varies by several orders of magnitude depending on the stage of the disease, but we currently know little about this for Ebola. It would be very interesting to compare the viral loads estimated from qPCR with the results from read mapping.
I've selected the first 20K (paired) reads from the samples so that I could map all samples (891). I wonder was the mapping percentage 0.25% a good approximation?. But you are right in that this may be sufficient. I have also noted low sequencing accuracy for the samples with low coverage.
The ratio of mapped reads is 0.3 % in my mapping, which is very close to your result obtained with only the first 20K.
Yep, it is a very small genome, you don't need much to get enough.
Comparing qPCR levels with coverage would be interesting but you'd have to be careful about how you calculated coverage, Ebolaviruses makes positive (antigenome) and negative sense genomic RNAs as well as positive sense monocistronic mRNAs. Although since this is serum I'm guessing you won't find much in the way of replicative intermediates.
If RT-qPCR and RNA-seq are run on the same specimen, it will not really matter, if it is plus or minus strand or both. I agree that you have to be very careful when comparing results from differing sample types (full blood, plasma, sweat) or if differing RNA extraction methods have been used.
It does matter, Ebola is a single-strand RNA virus. With a ssRNA genome, there is no positive and negative strand. The PCR assay may only be specific to viral genome, while RNA-Seq will detect genome, anti-genome and viral mRNAs. The number of copies of genome are not the same as the number of copies of anti-genome intermediate or viral mRNAs.
You'd want to see how viral genomic coverage and PCR levels play out, you don't want total viral coverage. qPCR is used to estimate the amount of virus present, if you were to add counts from anti-genome and viral mRNAs you could over-estimate this count.
The second point had to do with the life cycle of the virus. I would imagine that in serum you wouldn't find much in the way of products of replication (anti-genome, viral mRNAs) but instead would likely find genomic material since serum would presumably contain only mature virus. This is different than if you were dealing with tissues where you'd be sequencing intracellular RNA, i.e. where the virus is replicating.
I was curious to find out where the more than 2.5 million reads not mapping to Ebola come from. First I mapped them to some human genes (dhfr, actb, 18S rRNA), but does not got any hits. The serum seams to be free of human RNA. Then I mapped to some bacterial rRNA operons from several phyla (proteobacteria, firmicutes, actinomycete). I got diffuse mapping with thousands of reads for all of them. Thus the major fraction of RNA in this sample is bacterial rRNA from a broad spectrum of species. Presumably the serum was stored a room temperature for an elevated period of time.
I also noted that all of the read sets with high amount of reads mapping to Ebola belong to the 99 fully assembled sequences submitted to Genbank along with the Science paper. The other read sets seam to be just the trash of that great study.
That's pretty cool. Thanks for sharing, I wish more studies had posts like this were we can talks about it.
Also since we are talking about it, here is the original post that I was trying to reproduce: Mission Impossible: you have 1 minute to analyze the Ebola Genome As it turns out this time one needs to pick the right data.
Log in to answer this question.
Don't quote me on this, but I think there was a non-trivial number of cases where secondary infections or coinfections were seen, that might reduce the relative amount of EBOV reads.
Still, seems strange. Is there a way to organize these by date to see if there's a trend there?
Yes there are about five dates for data deposits, for example 195 datasets from Aug 2014, many of these show high mapping rates, then there a many hundreds of runs deposited in April of 2015 and those are mostly lacking all viral reads.
We can also see this from the SRR ids, the numbers associated with the run increase as time passes. The runs from August 2014 all start with 15 and the runs from April 2015 all start with 19.
The sample identifiers provided in the SRA entries seem to be associated with a viral genome under the associated BioSample. Although they list specific isolates for all of the samples there clearly are plenty of samples where there's no corresponding isolate genome.
Assuming they're still using the same criteria as in the original paper, all of these samples come from cases confirmed with qPCR. I wonder if these were places where samples may have been degraded or the depletion of host material didn't wory properly. Coinfection or secondary infections could confound this.
Excluding SRRs starting with 16 (only a few samples), the averages are all low, but the "15-Series" average is ~35% while the 17 and 19 series are both under 10%. The "15-Series" has the lowest percent (~36%, 82 samples) of samples under 5%, while 17 and 19 are ~90% and ~80% under 5%, respectively.
So it seems like something did change over time. Still, I wonder if it is possible to extract reasonable amounts of viral genomic data in poor mapping quality cases, but seeing how the number of genomes deposited is so much less it might not be possible to do reliably.
Interesting analysis. I think the discordance between the number of genomes and samples is due to the fact that most published runs do not contain useful information - but that in turn begs the question why is all that seemingly useless data uploaded and stored in a public repository under the same type entries as valid data?
Well I think I know the answer to that - many life scientists prefer to store everything "just in case" something big is hiding there ... but then this mentality was OK when datasets were tiny and the cost of the infrastructure and dealing with its complexity was localized to that person. It should not be OK to upload massive amounts of nonsense into a repository just because that nonsense happened to be collected.
Today we are being choked by data - alas I am afraid it may be mostly the useless and pointless data - as demonstrated by this project.
I was thinking the same thing, but it might just be a bit of semantics and a bit of awkwardness with how bioprojects work:
The project title: "Zaire ebolavirus Genome sequencing"
Under that: "Zaire ebolavirus sample sequencing from the 2014 outbreak in Sierra Leone, West Africa."
Of the two, the second best fits the data. The data is more "host-depleted blood/plasma sequencing samples from patients with qPCR confirmed ebolavirus infections in the 2014 Sierra Leone outbreak". I guess you could even call these metagenomic samples.
I think problem is that bioprojects are supposed to group how you got from A to B in your sequencing. If your goal is genomes, all of these samples might not be worth having.
I would surely not call most of this data useless. You could use this to establish the prevalence of secondary/coinfections, something that I would wager plays a major role in EVD. You could also do more mundane things like compare techniques. E.g. study how sample collection/transport/etc methods in an outbreak setting impact things, allowing you to better prepare for the future (although this would require more metadata).
I think there is a tremendous cost in lost opportunities and invisibly lost effort. Perhaps there is some major insight that this data has - or perhaps there isn't any. Clearly there should be some level of information embedded there with the data to help people make decisions beyond pure faith. It may be that putting data up the way it is done today may be more detrimental than helpful. It could be a bit like early screening for diseases - overall is more harmful than beneficial - it leads to wild goose chases.
Of course it does not help that the way SRA works is somewhat surreal - they have an interface where, if you wish so, you could choose to view every single read in each run (or you can page through the data, all the millions of pages ...) as if this was a use case that is useful in any way. But if you have a 32 bit computer you cannot download the data from command line (no really!! fastq-dump requires 64 bit computer) - as if downloading data would ever require over 2GB of RAM ...
Sure, it was a fishing expedition (I'd argue many sequencing projects are), but no one had ever done that before and given the challenging conditions/requirements to get this done, they managed to do very well.
In the supplemental for the paper from the first batch, they provide some additional metadata (location, survival, date). I'm not sure what the exact details are, they may not have been allowed to release more detailed data except under specific conditions. They also released sequencing stats, such as the number of EBOV reads, percentage of human reads remaining after depletion, and so on.
I don't really see the big issue, this is a very complex data set from a very complex situation.
Well just to clarify it - my beef is not with the project and whether their samples contain or not the ebola genome. I have nothing but respect for everyone involved and the work that was done. It is perfectly fine to have samples with no viral loads etc. that was and is a heroic task to put it all together.
My beef is with SRA, the rudimentary interfaces and information stored therein and limitations of it all, the way the current data sharing requirements are implemented and the hopelessness of something better coming about and how that sets back science.
Also a second observation - semi empirical. Some of the data for the runs that do not have good mapping percentages seems really bad. Initially I thought it must be that - so I trimmed and corrected the reads - but the mapping percentage did not change still zero -
the content of a few reads blasted against nr indicated a jumble of hits to everything: human, environmental, plants etc. notably none of the hits are full alignments rather than partial hits within the query.