I am working on quantifying RNASeq data at the exon level, which is intended for use in the DEXSeq package. I followed the documentation provided in the DEXSeq tutorial (http://bioconductor.org/packages/release/bioc/vignettes/DEXSeq/inst/doc/DEXSeq.pdf, page 23), as well as the commands on Subread to DEXSeq helper code (https://github.com/vivekbhr/Subread_to_DEXSeq).
First to produce the GTF File, I used the following code:
GTF="/Reference_Data/GRCh37/gtf/Homo_sapiens.GRCh37.87.gtf" python2 dexseq_prepare_annotation2.py -f Homo_sapiens.GRCh37.87_collapsed.gtf $GTF \ Homo_sapiens.GRCh37.87_collapsed.gff
Then I ran subread on the command line as submitted jobs on our shared HPC. We use the SLURM job scheduler.
#Program:featureCounts v1.5.3; Command: featureCounts -p -B -s 2 -T 16 -Q 10 -f -O -F GTF -a Homo_sapiens.GRCh37.87_collapsed.gtf -t exon -g gene_id \ -o sample_withJunctionsOnGenome_dupsFlagged_stranded.exon.cts \ sample_withJunctionsOnGenome_dupsFlagged.bam
I have two questions:
1) Why would I have exons for a gene that are only a few base-pairs long? I have some entries, in the example output below, where the length is calculated as only 1 base-pair. This does not seem representative of the true expression of a real gene, which would not have exons of such short length. Could it be there was some error in the commands run above or in the way the GTF was produced?
2) The number of assigned features for this quantification was particularly low. I have two sets of RNASeq data; 1) created with Poly-A enrichment in the library prep and 2) produced with the ribodepletion (RBS) library protocol. Both sets had paired-end, strand-specific, 75bp reads, and were processed through identical bioinformatic pipeline from raw FASTQ to BAM. I did check that the orientation of the paired reads were the same for both library preparation protocols. It was these BAMs that I input to Subread featurecounts() function.
For my poly-A library, I was getting on average 60% of my total reads assigned to features using the same commands as above, and same GTF reference.
For my RBS library, I am getting on average 30-40% of my total reads assigned to features, again using the same commands as above and same GTF reference. The vast majority of reads, from the summary files produced by Subread and then summarized into one figure by MultiQC, for all samples were classified as "Unassigned_NoFeatures".
I was thinking that this might make sense due to many additional non-protein coding RNAs that the Ribodepletion protocol will sequence. Has anyone seen this type of difference before? Again, should I be optimizing my commands for ribodeplete RNASeq that I am not aware of?
Thank you so much,
Example of my output for one sample is below: