DEE2 project update June 2019

It has been wonderful to be getting great feedback on DEE2, ways it can be improved and directions for future development of open omics data. Here I'll summarise some of these points.

Recent developments

  • iDEP integration. The iDEP service has provided DEE2 counts on it's R Shiny powered page which makes it very easy to download and analyse DEE2 data using the iDEP bioinformatics in the browser platform: http://bioinformatics.sdstate.edu/reads/ 

  • Degust integration. When using DEE2 in the browser you can now choose to send count matrices to Degust for differential analysis. You can send either the STAR gene counts, Kallisto Tx counts or Kallisto Tx counts aggregated to gene counts (recommended). Many thanks to David Powell's team at Monash especially Andrew Perry for working with me on adding this feature.

  • The R package getDEE2 now works for Windows systems. This required a modification to the download.file() options to account for default behaviour on these OSs. Thanks to Steven Ge for the bug report.
  • Each DEE2 download now contains a README to help end users understand the contents of the folder. Thanks to Peter Karp for the suggestion.
  • Backend dee_pipeline.R issue: I initially thought there was a bug with the backend R script dee_pipeline.R but it was just a memory overuse problem with a call to mclapply, so by reducing the parallel threads, the issue went away :) In addition, whenever load() is used it overwrites functions and variables like CORES, so in future I'll be loading individual objects from the RData file instead.
  • Data updates: Having pinned down the backend script issue, there are about 200,000 newly added mouse and human runs since last week. If there has been missing data in your project, it would be a good idea to try again.

What were working on now

  • Filling in missing datasets: You might have seen that there may be some datasets missing for particular runs in a project. We have scripts running to to search for these and and process them. 
  • SRAdbV2 is sadly not working anymore. This means that although SRA is growing rapidly, the database of corresponding metadata is static. We are looking to replace our dependence on SRAdbV2 but it will take a bit of time to integrate new tools. It could be that https://github.com/seandavi/OmicIDXR/ is the replacement.
  • More integration: Given the success with iDEP and Degust we are looking to integrate with other tools for downstream analysis, especially pathway level tools. 

Future directions

  • Suggestions for more species are welcome, especially if those suggesting can contribute resources (eg: computational, financial, etc)
  • Incorporating NGS data from sources outside SRA, such as GSA (China)
  • Transcript annotations: We've found that Ensembl cDNA used by Kallisto doesn't include ncRNA, which means most lncRNAs won't be properly represented. In future we may opt for the Gencode annotation as suggested by Alex Predeus. This can't be done immediately as it would require reanalysis of the entire dataset but is on the radar for future versions. 
  • Single cell data remains a tricky situation. Some submitted datasets appear to be processed very well by DEE2 pipeline, especially when the SRR runs correspond to individual cells which is typical of fluidigm preparation. However pooled runs from Drop-Seq or 10X seq seems to provide an average of all reads that make it through the pipeline. I think scRNA-seq will require a greater degree of manual curation due to the many varieties of library preparations which strongly impact the suitability of any analysis pipeline. Doing this properly will be a longer term direction for future development.
    Location: Laverton VIC 3028, Australia

    Popular posts from this blog

    Mass download from google drive using R

    Google drive is great for sharing documents and other small files but it's definitely not suited to moving  many large files around. For example I just received 170 fastq files that are about 200 MB in size. If you use the browser to download the whole folder, the web app will zip the contents for you which will take a LOOONG time. Alternatively you can download each and every one of those files one by one, which is annoying and prone to human error. You can insist to your collaborators to transfer in a different way, but there are not that many user-friendly and economic approaches. Your biologist collaborators probably won't be able to use rsync  to get the data to you safely. And fast convenient tools for moving around large files like Hightail cost a lot of money. A good solution to this problem is to use the R package googledrive  which enables command line automation of tasks that might take a long time manually. The package vignette has a good overview of the main comma
    Read more

    Data analysis step 8: Pathway analysis with GSEA

    Image
    In our RNA-seq series so far we've performed differential analysis and generated some pretty graphs, showing thousands of differentially expressed genes after azacitidine treatment. In order to understand the biology underlying the differential gene expression profile, we need to perform pathway analysis. We use Gene Set Enrichment Analysis ( GSEA ) because it can detect pathway changes more sensitively and robustly than some methods. A 2013 paper compared a bunch of gene set analyses software with microarrays and is worth a look. Generate a rank file The rank file is a list of detected genes and a rank metric score. At the top of the list are genes with the "strongest" up-regulation, at the bottom of the list are the genes with the "strongest" down-regulation and the genes not changing are in the middle. The metric score I like to use is the sign of the fold change multiplied by the inverse of the p-value, although there may be better methods out there
    Read more

    Extract properly paired reads from a BAM file using SamTools

    For a few different types of sequence analysis I need to extract read pairs that are properly mapped and satisfy some map quality filter, but I always forget the command line. So here it is: To read pairs from a bam file (-b) that map with mapQ≥30 including the bam file header (-h). The -f 0x2 option corresponds to the bitwise flags that specify that reads need to be properly paired. Proper pairing means reads are in Read1 forward, Read2 reverse orientation or Read1 reverse, Read2 forward orientation. $ samtools view -q 30 -f 0x2 -b -h in.bam > out.bam To extract single end reads from a bam file (-b) that map with mapQ≥30 including the bam file header (-h). $ samtools view -q 30 -b -h in.bam > out.bam
    Read more