I have performed an analysis of my RNAseq data using BWA software. The final output of this analysis is number of reads per target gene. I would like to compare these data with my qPCR data, however I do not how to normalize BWA data.
Thanks in advance for your help
Why are you using bwa but not tophat or star for allignment?
how lengthy your reads are?
you would probably loose lot of information by using bwa
Thanks for your reply Manvendra.
I am now finishing Tophat analysis, but I would like to obtain some conclusion from BWA analysis. We used BWA analysis because our reference genonme has many problems (assembly problems, missing data, etc..). We used Illumina Hiseq with paried-ends read. I think that the reads are 100 nt lenght.
which organism's genome you are working with?
you can fetch many informations from bwa allignment. But please remember that it would loose the informations concerning reads mapping onto splice junctions. If genome is not assembled/annotated properly and you have high depth with 100 bp paired end data, then you should probably go for de novo assembly, and identify all putative transcripts and then annotate them.
Thanks for your reply, Mavendra
We are working with Rhodnis prolixus. We also made a de novo assembly and I have validate all gene models I have included in BWA analysis (our target genes). My future plan is try to compare BWA analysis, TOPHAT and qPC results with these set of genes.
I don't think that it would be good idea to compare bwa and tophat and qPC. But , if you are insisting to compare then there is a way out.
If you have assembly and gene models (GTF or GFF).
1. get the list of genes which you have qPC, write the adjusted value in one coloumn.
2. get overlapping features from your alligned/mapped file (you must have got it in bam format e.g. tophat gives "accepted_hits.bam".
3. on your mapped/alligned files (bam format) count the reads in each co-ordinates (htseq-count is wonderful for this), divide by number of million of reads mapped on genome in total.
4. write those values in other coloumns side by side.
5. compute spearman's correlation for each pair (since its ranked correlation so it must work for this)
HTH.
Thanks a lot for the information, Manvendra
We want to use qPCR data because we do not have technical RNAseq replicates, but I know that correlation is not always good. We will try with your approach.
Sorry, I forgot to mention that once you have dataframe ready, please go for quantile normalization, make a boxplot to confirm that.
PS:
believe me, replicates make RNA-seq data analysis little complicated , It adds so much variance that finally you loose much informations , after filtering out on the basis of significance thresholds.
when you compute reads per million for a feature, then its already a normalized and comparable value
I agree with Manvendra that BWA is not the right tool for aligning RNA-seq reads: any that span intron will be not aligned properly. This means you will lose a fair amount of your data unnecessarily and will give spurious read counts as genes with more introns will be more affected than those with fewer.
TopHat is a massive improvement.
I disagree with the comment regarding replicates. Doing an RNAseq experiment without replicates is not going to give much information at all. You wouldn't do a qPCR without replicates? RNAseq is the same. You're not going to be able to say with any confidence whether your RNAseq results are comparable to the qPCR I'm afraid.
With tools like edgeR, Cufflinks, DESeq RNAseq replicates are very easy to deal with. Do more replicates.
As you have a poor quality reference, which will make read alignment somewhat challenging, have you considered doing a de novo assembly of your reads into transcripts? There are several tools available for doing this (e.g. Trinity).
This can be a really powerful way to validate gene models.
@Christian Cole, Thanks.
I agree to all of the above.
But in the case of replicates, if they are lane replicates then its good, but in the case of biological replicated samples, you always gain some and loose some, but if the genemodels exist expressing then anyways you would get in all the replicates .Moreover, during comparison variance level is also high due to replicates.
Exactly! You want to see the biological variance. That's the whole point of doing the replicates! Without them you can't say whether what you're seeing is biologically meaningful and if can't say that, then there's no point in doing the experiment. Unless you're not interested in biology.
If the variance is high, then that is what the biological background is. It cannot be changed it is simply a feature of the system.
Also, when I talk about replicates I mean at least 5. With 2-3 you will get massive variance (which is likely down to sample prep) and possibly poor results.
The reality in many cases is that there isn't the money/material for performing biological replicates to that depth. As you are trying to assess gene expression levels then you really want to get as many as feasible.
In terms of analysis then I agree with what has been stated before. A simple approach would be tophat2 alignment of the reads to the reference assembly (informed with the geneset gtf) followed by cufflinks to produce FPKM values for each locus. This will give you some idea on the variance between samples and you can compare with your qPCR results.
The question of whether to use the genome or a de novo transcript assembly as your reference depends on the representation of your candidate gene list in the genome assembly. Any (further) gene model changes can be accommodated locally in the GTF file used for cufflinks and/or tophat2 analysis.
Mail/call me if you want to discuss this further.
Dan
- Badges
- Science topic
- Similar topics
- Computer Science and Engineering
- Bioinformatics
- Bioinformatic Tools