bmmalone
commented
7 years ago Hi Tonu,
Given the counts of ORF types, there could be several things going on.
I'm not sure what sequencing protocol you are using, but it could be that the relevant adapter sequences are not present in the file given to flexbar. I believe the one included with the distribution includes standard TruSeq and ArtSeq adapters, but if you are using something else, that could be a problem.
Something else (which you may have already checked) is that the number of uniquely-mapping aligned reads in the <riboseq_data>/without-rrna-mapping/<sample-name>[.<note>]-unique.bam file are roughly the same.
Likewise, while it isn't amazingly helpful for riboseq, you can also try FastQC and just see if anything comes up. I would run it on the fastq file after removing ribosomal rna: <riboseq_data>/without-rrna/<sample-name>[.<note>].fastq.gz. Of course, the read length distribution throws up flags, but most of the other results there should look relatively normal.
It could also be that the semantics of start/end are different on the reverse strand. To check that, you can load the annotated bed file (the <genome_base_path>/<genome_name>.annotated.bed.gz file) into something like IGV and make sure the "thick" parts match up with starts and stops in the genome sequences. In particular, the start codon should be included in the "thick" part and the stop codon should not be included. I attached example images from IGV to show what I mean.
When this "off by one" kind of problem has turned up before, we usually ended up with many "within" and "suspect" ORFs predicted as translated, but not so many "canonical_truncated" ones. The "truncated" ORFs are in-frame, but there is better signal starting from a down-stream start than the canonical one. Here, you could open the uniquely-mapping read bam file in IGV and just take a look at the signal. Unfortunately, we do not currently create a bam file which removes the non-periodic reads and adjusts for the P-site offsets (we just directly do that when creating the ORF profiles), so there is not an easy way to view that. However, there still could be something here.
Depending on how that looks, you could also look into the "unfiltered" predictions bed file. In ipython, you can look at the count of the different ORF types here with something like:
import misc.bio_utils.bed_utils as bed_utils
unfiltered_predictions = bed_utils.read_bed("/path/to/my/unfiltered-predictions.bed.gz")
orf_type_groups = unfiltered_predictions.groupby(['orf_type', 'strand'])
orf_type_groups.size()(unfiltered_predictions is just a pandas data frame here). The (so far undocumented) create-orf-types-pie-chart script can also be used to visualize these counts.
These are the predictions before the filtering described in the "Final prediction set" paragraph (the last one in the methods part) in the paper. If there are many more canonical ORFs on the negative strand here than in the final set, then it means the periodic signal is stronger in the truncated (or "within") ORFs than the canonical one, so those are selected instead.
There are other things to look at, as well, but I think this covers most of the easy and not-so-difficult things.
We have run data from a number of different species and using several different protocols, and we haven't encountered any protocols that have a noticeable strand bias; however, if that were the case, it could still be handled within the current framework. However, it would require some code to handle alignments from the different strands separately. In essence, though, alignments from different strands could be treated as replicates (with independent periodic read lengths, P-site offsets, etc.).
Please let me know how this goes, or if there is any trouble checking anything. Hopefully something here gives some clues, though.
Have a good day, Brandon
reverse strand end:
reverse strand start:
tmargus
Hi Brandon,
I found something unexpected when plotted out metagene profiles using included IPy-notebook file "create-orf-type-metagene-profiles.ipynb" by using my own data from yeast.
Most of ORFs are tagged as "canonical", canonical_truncated", or "within". In forward strand majority ORFs are "canonical" (2,563) and "truncated/within" (121/114) forms smaller group. Surprisingly in reverse strand there is a much less "canonical" (239) while "truncated" and "within" forms the major class of ORFs (1,851 and 836 respectively). orf-metaplot
For reverse strand there is no ribosomes pilup at ends when in forward strand its clearly visible. Could it be that coordinates become a bit shifted in reverse strand?
What kind test could be best to run to figure out what went wrong if any?
Cheers, Tonu