Currently, I am interested in several (around 100) genes in fish and would like to investigate their expression level using public available RNA-Seq data. My strategy is to build up the reference sequences (interested genes). Index them with bowtie 2 and then align the public available RNA-Seq SRA data (filtered using SRA tool kit) against it. The obtained SAM file was further counted by eXpress for each gene expression level using the FPKM value.
I have several questions about this strategy,
Firstly, when building up the functional gene reference, what kind of sequences should I use if there is no genomic data available? For example, gene A may studied by several scholars and their sequence results can be found in the NCBI Nucleotide database but with difference lengths. Which one should I choose. Besides, RNA splicing proceeded during RNA expression, introns may be spliced out. Therefore, which sequence should I use before or after splicing (this is important because the length of the gene affect the final FPKM value) and how I can identify whether the obtained RNA sequence is spliced or not.
Secondly, is there any problem with the estimated expression level using this strategy? Over or underestimated.
Any other suggestions are strongly welcomed!
Dear Jian,
I understand that you are working with a species without genome and transcriptome sequenced, but in stead you want to use ESTs or other sequenced clones, right?
Now, the available RNA-Seq data that you mention is from this particular species or for other fish?
If it is the first case, I suggest to build a de novo transcriptome using Trinity for instance. With this tool you can also obtain expression level for the recovered transcripts. Then you only have to find which one of these transcripts correspond to your genes of interest.
If the answer to my second question is the second case, I think that you have a very difficult problem. You will try to map sequence from different species, meaning that you have to allow higher level of mismatches. This will introduce a lot of imprecision. Another source of imprecision is that since you won't have the full set of genes to map, but the RNA-Seq is genome-wide, you might end mapping reads to one of your genes that came from a paralog or related gene.
Best,
Juan
Dear Juan,
Thanks a lot for the wonderful suggestion. My study is the second case. 5~6 species were studied, I need to build different fish specific reference sequences. Doing de novo is too much time cost (>100 G data). However, as you mentioned, I may map the read to a paralog or related gene, this can cause over estimation. Can I bring in the parameter of coverage, to see where the aligned reads is, if less than some percentage, then it is determined to be not expressed. And any recommended software to do the coverage calculation?
Best wishes,
Jian
Dear Jian,
To compute the coverage over your genes you can use bedtools (http://bedtools.readthedocs.org/en/latest/content/tools/coverage.html, http://bedtools.readthedocs.org/en/latest/content/tools/genomecov.html), GATK (https://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_tools_walkers_coverage_DepthOfCoverage.php) or simply using samtools.
Best,
Juan
I think Juan's idea of building your own transcript assembly and then mapping to that is a good one. Trinity and Bridger are nice software for that. If you can get access to a large server/cluster you could probably run transcriptome assemblies in a matter of a few days. If not, I would suggest taking the full length genes (with introns) and map with a splicing-aware mapper such as Tophat2 or GSNAP. You could use Bowtie if you're not interested in the splicing. A program like HTseq lets you get raw read counts for your genes of interest, and this can be used to calculate FPKM or be put into differential expression analysis software (eg DESEq). Cufflinks can also be used for quick gene expression analysis.
Good luck
Thanks Cameron, do you know any method to batch download a group of RNA sequences? For example, I am interested in a mineral absorption KEGG pathway, I would like to have all the coding sequences involved in the pathway. Or any other website than can batch download an assembly of functional gene sequences.
Thanks Cameron, do you know any method to batch download a group of RNA sequences? For example, I am interested in a mineral absorption KEGG pathway, I would like to have all the coding sequences involved in the pathway. Or any other website than can batch download an assembly of functional gene sequences.
- Badges
- Science topic
- Similar topics
- Computer Science and Engineering
- Bioinformatics
- Bioinformatic Tools