This is a test version of Biostars. For the public version, visit https://www.biostars.org.
eXpress uniq_counts estimation

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
express rna-seq non-model-organism bwa

0 answers

No answers yet.

Log in to answer this question.