hms-dbmi / dseqr Goto Github PK

IBD data is too large for github.

Going to just repackage first 1000 or so reads from each sample from IBD data.

Go ahead with:

• GRAS
• SIDER links
• Drugbank links
• Wikipedia links

Also see #4 for further discussion

Probably required exprs_sva to be run through limma::voom.

See issue #2 (detect single/paired end experiments) for some related discussion.

Will put this on hold until we have some example paired data to work with.

Needs to match CMAP02/L1000 annotation as I currently have it (HGNC symbol).

I got about a 0.45 pearson correlation between the logFC values for my analysis run with the current drugseqr pipeline and your Biojupies analysis. This seems kind of low for two different pipelines on the same dataset (maybe not?).

Here are possible differences/unknowns that I can gather based on looking at Sam's report that might contribute:

• quantification of raw data (Alex - salmon, Sam - ?)
• library normalization (Alex - edgeR::calcNormFactors, Sam - logCPM?)
• surrogate variable analysis (Alex - Yes, Sam - No)
• limma analysis (Alex - limma::voom transformation, Sam - ?)

I'll try running a couple variations on my pipeline (e.g. without sva) to see how much difference it makes.

Let me know if you have any ideas samuelfinlayson

┆Attachments: mds-sva.png

Should show CMAP02/L1000 results based on selection

Started with Broad Repurposing Hub Data (Thanks samuelfinlayson!)

Near term:

• SIDER side effects link out
• Drugbank linkout
• GRAS

Longer term:

• info on approved drug with closest RDKIT morgan similarity

samuelfinlayson suggested assessing similarity measures like in my previous blog post but instead of using external GEO datasets, using LINCS data from one cell line to try and retrieve the same perturbation in a different cell line.

This approach has some nice advantages:

• no need to gather/label external data
• creates a larger training dataset for possible neural networks

Related to #9

causes malformed links for #6

Incorrect values include:

• 0
• -666
• MLS003116075

Like #36, plan to fix in private l1000 package and propagate.

e.g. moa column

started happening after 03346dc

Seems to be related to UTF-8 issues

samuelfinlayson expressed interest in training a Siamese neural network that could take as input a LINCS/CMAP02 signature and a query signature and produce a ranked list of compounds. This would contrast with more classical signature-similarity approaches (pearson correlation, cosine similarity, KS).

He mentioned this in the context of his PhD work as well.

The problem is that it takes a bit of time to summarise the ~230,000 signatures by compound name. To resolve this I am going to first subset the results, keeping all entries for compounds that have a correlation value less than the 100th top result with a clinical phase. This way, when filtering by clinical phase, you will still see 100 results. When not filtering by clinical phase, you will also see all compounds that have a correlation value less than (stronger negative correlation) the 100th top results with a clinical phase.

To do:

• show all correlations for a compound
• sort based on all signatures for same pert (see #40)?
• show a tool tip on hover that indicates cell line, dose, etc.

l1000_es and cmap_es are dprime effect size values so they should be queried using equivalent measures.

e.g. 'WZ-4-145' has cid 73707531 for some entries but none for others.

Fill in the missing cids using non-missing entries. This will be done in my private l1000 package so that the fix propagates to other code that uses it.

Sorting measures can be evaluated using the same benchmark/data that I have used in the past.

I have mentioned to samuelfinlayson that performing differential expression on L1000/CMAP02 at the pert level (ignoring dose, time, and cell line) results in worse benchmark results. This suggests that any sort of resorting should be carefully benchmarked as it could also decrease performance.

rnama.com does something like the following:

• perts are first sorted by decreasing percentage of consistent signatures
• ties are first broken by number of total signatures
• ties are then broken by largest absolute correlation value

This sorting on rnama.com has not been benchmarked.

isaackohane1 mentioned this stackexchange as well but I'm not totally clear on the how the stats work when you wants to get a confidence interval across multiple correlations (for the same compound).

Empirical fdr/pvals can also be derived across multiple correlations using fdrtool. This approach could be used to derive pvals and then do some sort of p-value meta analysis (lots of methods in metap package). There are also correlation value meta-analysis methods (e.g. meta::metacor).

I think that unless we benchmark and find that resorting improves connectivity mapping, the cautious approach is to maintain the original sort (i.e. minimum correlation).

Requires porting/adapting code from crossmeta::diff_expr and private repo rnama::diff_expr (initial rnama.com website which would run sva::svaseq on RNA-seq data).

For example, in l1000_es, Cell Migration Inducing Hyaluronidase 2 is annotated as TMEM2 whereas get_tx2gene returns CEMIP2. This is likely related to changes in official HGNC symbols (TMEM2 used to be the official HGNC symbol for what is now CEMIP2).

I'm not sure what the best approach is to resolve this. Possibilities:

1. save a reference entrezid --> HGNC mapping and always use that. The downside of this approach would be that over time the official symbols may differ. I think I will go with this and use the mappings that I used for cmap_es and l1000_es. This was downloaded from NCBI ftp://ftp.ncbi.nih.gov/gene/DATA/gene2accession.gz on 2017-08-12.

2. save everything at the entrezid level and only map to HGNC as needed. The downsides to this are increased complexity and loss of interpretability. Specifically, mapping from entrezid --> HGNC would have to occur at multiple different locations instead of once and all the data would be annotated with entrezid numbers instead of gene symbols. There would also be a loss of integration with other code I have developed that may be useful in the future. If everything was from scratch I might do this but I think it is probably a poor choice now.

This could be useful for datasets where we only have disease samples (e.g. chronic fatigue patient).

isaackohane1 mentioned some previous work where control samples from a different study could be used successfully if gene rank-based differential expression was used.

Some compounds with the same Pubchem CID have a different clinical_phase. I think it makes the most sense to show the most advanced phase only. This will simplify sorting for #37.

isaackohane1 and samuelfinlayson were discussing the contrast between the specificity of cripr and antibodies as compared to small molecules.

The question came up regarding whether targeting a gene at the protein level through a monoclonal anti-body would produce similar changes in gene expression as those caused by cripr at the genetic level.

isaackohane1 mentioned the possibility of applying for funding to create a CMAP-like dataset for monoclonal antibodies.

Rashmi mentioned something like 26 SJIA patients that blood samples have been taken from (expecting ~40 at upcoming family day) including some unaffected sibling control samples.

If running an analysis using multiple patients with controls that are specific to each patient this can be specified in the model matrix.

e.g. atorvastatin, canrenoic-acid

Seeing the value of these ids (Pubchem CIDs used for links to both Pubchem and SIDER) tells you nothing useful other than that they exist.

I think it makes more sense to just have a single External Links column with an image linking to the respective data source. This will reduce table width (single column vs 3), be visually appealing, and be just as informative.

Related to #43

Possible sources:

• FASTQ files of bulk rna-seq reads
• single-cell rna-seq data
• affy microarrays

Will initially assume bulk rna-seq with biological reps as a starting point.

Need a raw data format standard for distinguishing if multiple files belong to the same sample and if the data is single-end or paired-end.

The SRA (GEO) uses the format SRR_1.fastq.gz where:

• SRR is the sample name
• _1 indicates that file contains reads for one end in a paired-end experiment
• absence of either _1 or _2 indicates single-end reads.

Some standard like above would be useful for auto-determination of quantification flags.

Would be nice to see:

• before and after sva plot (side by side)
• ability to hover over samples to inspect outliers
• use ggplot2/plotly to improve appearance/fix legend positioning (e.g. here)
• add some jitter/alpha so that overlapping points are visible.

Related to #8

Caused by difference on mod0 (uses selected and non-selected samples) and mod (uses selected samples only).