Question: Generating FPKM matrix accross all samples after stringtie
0
gravatar for ever_wudi
10 weeks ago by
ever_wudi10
ever_wudi10 wrote:

Hi,

I finished the "stringtie --merge" step and used "stringtie -e -B -G" to reestimate fpkm for each sample and generated the ballgown folder. Now I am trying to get FPKM for all genes and across all samples into a matrix. I know there is a prepDE.py script, but it generates gene/transcript counts instead of FPKM values, and I really need an FPKM matrix. Is there a good way to output such matrix in ballgown, stringtie, or a software that can be directly connected to the outputs of stringtie?

Thanks! Di

rna-seq • 467 views
ADD COMMENTlink modified 10 weeks ago by manuel.belmadani60 • written 10 weeks ago by ever_wudi10

Why do you need FPKM? It's an outdated form of normalisation for RNA-seq, and effectively renders your samples incomparable in terms of differential expression analysis.

The idea of using prepDE is to generate 'raw' counts for the purpose of using other [more robust] differential expression analysis tools, like EdgeR and DESeq2.

ADD REPLYlink written 8 weeks ago by Kevin Blighe28k

Thanks for reminding. I didn't intend to feed the count matrix to EdgeR or DESeq2, so I just needed the matrix in the FPKM form. Is there something wrong with the FPKM normalization method?

ADD REPLYlink written 6 weeks ago by ever_wudi10

Yes, it renders the samples non-comparable in terms of differential expression analysis (you can compare them, obviously, but the false-positives / -negatives will be higher than other methods). Not a life or death situation, for now, but just be aware that there are proven better methods.

ADD REPLYlink written 6 weeks ago by Kevin Blighe28k

Thanks for the tips. What are the better metrics used now? TPM or some others?

ADD REPLYlink written 6 weeks ago by ever_wudi10
1

TPM is better but still not ideal. You should use those methods employed by either DESeq2 or EdgeR. Here is the manuscript that I had wanted to share with you (it is a very well cited manuscript):

In their key points:

The Total Count and RPKM normalization methods, both of which are still widely in use, are ineffective and should be definitively abandoned in the context of differential analysis.

Note - FPKM is essentially the same as RPKM

ADD REPLYlink modified 6 weeks ago • written 6 weeks ago by Kevin Blighe28k

I see, thank you very much!

ADD REPLYlink written 6 weeks ago by ever_wudi10
3
gravatar for ATpoint
10 weeks ago by
ATpoint7.5k
Germany
ATpoint7.5k wrote:

If you did use stringtie properly, there should now be a folder "ballgown" that contains subdirectories for each sample, and each of these subdirectories should have any kind of common pattern in its name, like "Sample1_ballgown", "Sample2_ballgown". If so, load everything into R and get FPKM, which is as easy as:

require(ballgown)

## Create ballgown object:
my.data = ballgown(dataDir=/Path/To/ballgown_Folder/, samplePattern='*_ballgown', meas='all')

## get FPKM:
gexpr(my.data)

There are a lot more one-liner commands to extract common things from the data, such as exon coverage and stuff, which can easily be looked up in the BioC tutorial. A tip for the future: Really try to find these kind of answers yourself, including proper googling and a lot of try and error if you have the time. You learn so much by trying-out things and get good practice in finding solutions. A provided solution may help you know, but does not really challenge you to improve yourself.

ADD COMMENTlink modified 10 weeks ago • written 10 weeks ago by ATpoint7.5k

Thanks ATpoint, the commands worked. I am very new to R and RNA-seq, did try to search for solutions before but didn't find a satisfactory way. Thanks for the help.

ADD REPLYlink written 9 weeks ago by ever_wudi10

You're welcome. As you'll see when you get more experience with R and Bioconductor, most commonly there are simple commands like this for standard tasks. You only have to find them ;-D Please consider to mark the answer as accepted to help others in the future.

ADD REPLYlink modified 8 weeks ago • written 8 weeks ago by ATpoint7.5k

I just saw this message, and just realized there is a "tick" sign for marking an accepted answer. I already marked it accepted.

ADD REPLYlink written 6 weeks ago by ever_wudi10

Hi ATpoint, I have another question. When I use gexpr() to output the matrix, it generates a matrix with MSTRG.x gene ids as identifiers. However, I really need gene_name (or gene symbol) as my gene identifiers for the next step. I tried really hard to fiddle it around and search online but still can't do it. Would you have any suggestions for it? Thanks.

ADD REPLYlink written 6 weeks ago by ever_wudi10
2
gravatar for manuel.belmadani
10 weeks ago by
Canada
manuel.belmadani60 wrote:

Edit: The other answer is probably a more proper way to get the gene-level matrix, using gexpr https://www.rdocumentation.org/packages/ballgown/versions/2.4.2/topics/gexpr and should probably be marked as accepted. I'll leave mine here in case the info is useful, but I wrote my answer thinking of transcript per sample matrices. texpr seems to do that as well.

I think you can get your FPKMs with Tablemaker per sample with your ballgown output, but you might have to assemble your matrix yourself (should be very few lines of Bash and/or R).

See https://bioconductor.org/packages/devel/bioc/vignettes/ballgown/inst/doc/ballgown.html

Under Running Tablemaker

tablemaker -p 4 -q -W -G merged.gtf -o sample01_output read_alignments.bam

The output is 5 files, written to the specified output directory:

t_data.ctab: transcript-level expression measurements. One row per transcript. Columns are:

  • ...
  • FPKM: Cufflinks-estimated FPKM for the transcript (available for each sample)

Here's an example t_data.ctab: https://github.com/alyssafrazee/ballgown/blob/master/inst/extdata/sample02/t_data.ctab

Note that each sample should have it's own t_data.ctab file, so all you need to do it grab the FPKM column for each sample and merge them together.

You can get the relevant info from the ctab with cut

$ cut -f6,12 sample1.ctab
t_name  FPKM
TCONS_00000010  800.706
TCONS_00000017  715.775
TCONS_00000020  579.569
TCONS_00000598  2205.3
TCONS_00000024  304.873
TCONS_00000613  101.848
TCONS_00000029  87.7438
TCONS_00000032  0.0249079
TCONS_00000637  323.543

And join with join

$ join -j1 <(cut -f6,12 sample1.ctab | head | sort -k1) <(cut -f6,12 sample2.ctab | head | sort -k1) 
TCONS_00000010 800.706 800.706
TCONS_00000017 715.775 715.775
TCONS_00000020 579.569 579.569
TCONS_00000024 304.873 304.873
TCONS_00000029 87.7438 87.7438
TCONS_00000032 0.0249079 0.0249079
TCONS_00000598 2205.3 2205.3
TCONS_00000613 101.848 101.848
TCONS_00000637 323.543 323.543

For more than 2 files, you can join the first 2, then join the joined output with the next file etc. Or this thread had more elaborate solutions: https://stackoverflow.com/questions/10726471/join-multiple-files

Or R could do something similar if you load them as dataframes and used merge(all.x=T, all.y=T) iteratively.

Whatever you end up doing, make sure that you test your results. Check how many rows you end up with. I'm not familiar with the stringtie ctab files, so I'm not sure if they would all have the same number of rows or different; that's something to check too. Spot check some transcripts, ideally a random sample for top, bottom, middle of your files; make sure they match with the appropriate sample values.

Hope it helps!

ADD COMMENTlink modified 10 weeks ago • written 10 weeks ago by manuel.belmadani60
1

Thanks for the detailed steps, manuel.belmadani! I really appreciate your help.

ADD REPLYlink written 9 weeks ago by ever_wudi10
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 775 users visited in the last hour