Hi,
I am dealing with extremely polymorphic gene (no alternative splicing involved) in one of Phytophtora non-model species (no reference genome available). I have several cloned variants (~20) of this gene (Target sequences) and I want to check, which of them are really expressed during infection (if any at all).
I have RNAseq data in my disposal (Illumina PE, 76, insert size ~ 220) of a mixed sample (pathogen + plant).
After QC and adaptor trimming I aligned RNAseq reads to the available plant reference genome. Then I extracted the unmapped fraction of reads (assuming that vast proportion of them might belong to the pathogen) and aligned them to my target sequences using BWA.
The problem is that the target sequences are highly similar, so there are lots of multiple mapping reads in my bam file and only 8 uniquely mapping reads (all mapping to the same target sequence). I wanted to rescue some of those multiple multiple mapping reads and to estimate expression more accurately. My first guess here was to use something based on expectation maximization algorithms, so I opted to eXpress (after conversion of bwa produced alignment to ~bowtie format and sorting of bam file).
Generally eXpress results were concordant with estimations based on unique reads only (one major variant expressed), though the numbers of reported fragments uniquely mapping to targets are significantly higher than I get just parsing bwa output. And I also got 'uniq_counts' reported for larger umber of target sequences (6 instead of just 1).
I wonder how does eXpress estimates those 'uniq_counts' and is it correct to use this pipeline to quantify abundance of similar transcripts in my case?
Any suggestions regarding the general approach would be also appreciated.
Here is my eXpress output for target sequences with mapping reads != 0:
bundle_id target_id length eff_length tot_counts uniq_counts est_counts eff_counts ambig_distr_alpha ambig_distr_beta fpkm fpkm_conf_low fpkm_conf_high solvable tpm
7 A 420 318.27736 67 20 60.216406 79.461795 3.940519 0.6646761 1244702 901097.5 1588307 T 396069.1
7 B 420 318.27736 37 36 37 48.825339 2.999997 0.000003 764807.9 746750.9 782865 T 243364.9
7 C 420 318.27736 24 23 23 30.350886 0.000002 1.999999 475421.2 460426.9 490415.4 T 151280.9
7 D 420 318.27736 20 19 19 25.072471 0.000002 1.999999 392739.2 376868.1 408610.3 T 124971.1
8 E 417 315.293273 5 5 5 6.612891 0 0 104330.6 104330.6 104330.6 T 33198.4
7 F 417 315.293273 43 1 3.899894 5.157915 1.774216 23.92226 81375.66 0 176473.1 T 25894.05
7 G 426 324.245535 45 0 3.081816 4.048949 1.233036 16.77149 62530.11 0 145695.7 T 19897.32
7 H 417 315.293273 27 0 0.801884 1.060554 0.4860403 15.87929 16732.2 0 56190.54 T 5324.251
0 answers
No answers yet.
Log in to answer this question.