Correlating back to cancer type or other annotation fields would be important. Visualization as a bar or pie-chart to show % of samples with some germline mutations.

Should we attempt this for the paper or tackle that in OpenPBTA2? 😄

We have generated CNV output from ControlFreeC and CNVKit, but are seeking individuals to determine consensus focal calls and/or identify additional algorithms we can run to instill high confidence in focal CNV calls from the WGS dataset.

Hi @cbethell! I think this is a great start and one thing we have been struggling with is that some of the RNA-Seq was a poly-A prep while the majority was stranded. I am not sure anyone has yet to figure out how to normalize across these two batches (eg I think UCSC Treehouse folks are also very interested in this). In some of our discussions, it has been suggested to subset the dataset for only those genes which will be captured in a polyA prep, and see if that helps, but you will lose a lot of data in that way.

Originally posted by @jharenza in #83 (comment)

The idea is to demonstrate that samples in the PBTA cluster by cancer type / molecular subtype using a dimensionality reduction technique such as PCA/T-SNE/or UMAP. Please see figure 2 (https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1006599) as an example of using this to look at GTEx and differentiate different tissue types.

Characterize PBTA with COSMIC Mutational Signatures and correlate back to cancer type, genetic lesions, and pathway activity. Visualization of this would be important as well.

Scientific goals

As previously noted (e.g. in #83 (comment) ) the RNA-seq data was compiled from experiments that used two different techniques, poly-A prep vs. stranded. This is confounded with the source of the data: PNOC vs. CBTTC. If we are going to be able to compare these two sources, we will want to investigate if there is a way to identify comparable parts of the data set.

Proposed methods

Filtering of expression data to identify gene sets suitable for cross-method comparison. It may be that no such comparison can be made at this point, in which case all analysis between the two data sets should be expected to be performed in separately, in parallel.

Required input data

RNA-seq data (kallisto and RSEM outputs)

Proposed timeline

What is the timeline for the analysis?

Relevant literature

If there is relevant scientific literature, put links to those items here.

Issue

I made a PR checklist inspired by the OpenPBTA project guidelines. This has been helpful for me personally, and I therefore believe that it can be helpful to other contributors. It would serve as an additional and optional resource.

The checklist is in markdown format and can be found below.

OpenPBTA-analysis

PR checklist

• Fork the repo.
  You only need to do this step once, continue to the third step of this list if you have already done this.
• Make the original repo, AlexsLemonade/OpenPBTA-analysis, an upstream repository.
• Checkout a new branch to make changes on.
  WAIT Is your master branch up to date with the upstream repo?
  Try git status to check.
  If it's not, go here before creating your new branch.
  Remember to git push when you are done.
• Did you make numerous changes, too many for one pull request?
  Then let's create stacked pull requests.
• Now that you have committed changes, be sure to write a detailed summary of the changes by following this template.

Now check to ensure you have done the following before filing the PR (remember to push any additional changes):

• Run a linter (using styler for example)
• Set the seed (if applicable)
• Structure of analysis directory meets the Project Structure guidelines
• Comments and/or documentation up to date
• Double checked paths
• Spell check any Rmd file or md file (using mdspell for example)
• Restart R and run all notebooks fresh and save (be sure to run in the Docker container)
• Connect the pertinent issues on ZenHub
• Updated Dockerfile and built Docker image successfully
• Added analysis to continuous integration

If you have completed all of the applicable steps above, you can request a suitable reviewer and file the PR/draft PR.

Note: Upon approval of all Pull Requests for an analysis, create a final Pull Request for a shell script that runs the entirety of scripts in your analysis directory.

Question

Is it necessary to add this checklist to the repo?

Determining genetic lesions (Mutation, CNV, Fusion) and/or pathways which co-occur or are mutually exclusive across the PBTA. This could help associate lesions with pathways or define potential synthetic lethality relationship.

This file will also need to be added since #45 has landed.

In service of making the project more approachable for new contributors, I'm creating this issue to discuss the use of our issue labels.

From the Planned Analysis section of the README:

We have tagged a subset of these with the label "good first issue". These "good first issues" are ones that have a limited scope, that can be done with the software already present on the Docker container, and that we expect to have few dependencies with other issues.

We do not currently have a 'fleshed-out' Docker image. As a result, I've untagged all issues that were labeled as good first issue for the moment. When we have the container building and ready for use, we should focus on making limited-in-scope issues and tagging those.

I've also recently added the blocked label, which is meant to indicate that an issue is blocked by something that is external to or upstream of this project. We should apply this label to existing issues where appropriate.

Another consideration: should we use issues to indicate priority?

Context is #76 (comment), to quote that here

• As individual notebooks or scripts are coming in for an analysis module, they should be added to CI as multiple runs. This validates that they can be executed at the time of review.
• Once an analysis has been reviewed and validated, the last PR for an analysis module should be a shell script used to run everything in the module. That PR also replaces the multiple runs in the config file with a single run that runs the shell script.

I propose that we try this model out with @cbethell's addition of #5 (relevant PRs: #54, #55). If successful, we can update the instructions in the README.

File(s)

What data file(s) does this issue pertain to?
pbta-sv-lumpy.tsv.gz

Release

What release are you using?
V2/V3

Link to OpenPBTA-manuscript

Put a link to the relevant section of the OpenPBTA manuscript here.

Question/issue

Put your question or report your issue here.
LUMPY file only contains 171 participants and should include ~700 - we are remerging this file and will include with V4 data release. cc: @yuankunzhu @kgaonkar6 @jaclyn-taroni

Propose dividing primary_site into:

midline = Thalamus|Pons|Spinal Cord|Brain Stem|Midbrain|Hippocampus
hemispheric = Lobe|Cerebellum
other = all other brain regions

2019-08-13-brain_regions.xlsx

@awaanders will help assess accuracy of this prior to updating the clinical file on S3.

Determine immune/stromal composition using software such as xCell and correlate back to genetic and phenotypic characteristics.

File(s)

What data file(s) does this issue pertain to?

Release

What release are you using?

Link to OpenPBTA-manuscript

Put a link to the relevant section of the OpenPBTA manuscript here.

Question/issue

Put your question or report your issue here.
We had some confusions about how to use symlinks with R, and the R.utils package installed with the docker should have been able to read symlinks, but was unable to (with #92), so file.path had to be used instead. May be good to add some of the detail around use of symlinks for specific types of files with R to the README.

Scientific goals

In #73 it's noted that the reported sex of the participants in the study didn't align with information from germline sequencing in 11 cases. It could be interesting to understand how accurately this can be called from gene expression data. For some studies, gene expression data is all that is available. If it can also be accurately called directly from gene expression data, these studies will be better positioned to evaluate their metadata.

Proposed methods

• Construct an elastic-net logistic regression classifier using gene expression values as features and reported sex as the labels. I suggest elastic net because the signal is expected to be relatively linear, most genes are expected to play little role, but a few sets of genes are likely to be predictive and highly correlated, and it makes sense to spread weights across them to produce a robust predictor.
• Evaluate using cross-validation with reported labels across the full set.
• Evaluate using cross-validation with germline-based labels across the full set
• Evaluation prediction accuracy using germline sequencing within each histology.

Required input data

For classifier construction:

• The reported sex in the PBTA histologies file.
• Gene expression estimates from kallisto, RSEM, or both.

For evaluation:

• The germline-based sequencing calls.
• Histologies, so that performance can be broken out by histology as well.

Proposed timeline

I am proposing this analysis but don't have time to do it, so I would leave this as an estimate for someone who decides to take this on.

Relevant literature

There are a number of reports of this being readily discoverable, even with unsupervised methods. These are two from our group where we noticed this, but there are also others from other groups:

• Unsupervised VAE -> feature aligned with sex: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5728678/
• Multiple unsupervised methods -> feature aligned with sex: https://www.biorxiv.org/content/10.1101/573782v1

The idea would be to visualize trends of deletions / amplifications in each of the cancer types (where relevant). Pointing out highly amplified/deleted regions with oncogenes/tumor suppressors may be interesting.

The idea here is calculate pathway scores on each of the samples and summarize via a visualization to demonstrate associations of different pathways with phenotypic or genetic elements. This visualization be in the form of a clustered heatmap and should include some sort of statistical analysis when identifying "strong" associations.

Working on #2 - @yuankunzhu - do you have a checksum for these files that we're downloading that can be downloaded at the same time? I'd like to make sure that when folks download they are getting the expected contents.

While working on #33, I read the Arriba file into R like so (removing the full path info for simplicity):

arriba <- readr::read_tsv("pbta-fusion-arriba.tsv.gz")

And the output in the console was:

Parsed with column specification:
cols(
  `gene1 gene2 strand1.gene.fusion. strand2.gene.fusion. breakpoint1 breakpoint2 site1 site2 type direction1 direction2 split_reads1 split_reads2 discordant_mates coverage1 coverage2 confidence closest_genomic_breakpoint1 closest_genomic_breakpoint2 filters fusion_transcript reading_frame peptide_sequence read_identifiers tumor_id` = col_character()
)

Running

arriba <- readr::read_delim("pbta-fusion-arriba.tsv.gz", delim = " ")

yields

Parsed with column specification:
cols(
  .default = col_character(),
  split_reads1 = col_double(),
  split_reads2 = col_double(),
  discordant_mates = col_double(),
  coverage1 = col_double(),
  coverage2 = col_double()
)
See spec(...) for full column specifications.
1 parsing failure.
  row col   expected     actual                                   file
37269  -- 25 columns 24 columns '../../data/pbta-fusion-arriba.tsv.gz'

and

> dim(arriba)
[1] 54279    25

Also checked by gunzipping the file and using less -U. I get a similar pattern with pbta-fusion-starfusion.tsv.gz.

I would expect files with a tsv.gz extension to be tab-delimited. Is it possible to change these files to be TSV files in the interest of smoothing the way for folks using a TSV parser on them?

EDIT: If this is a limitation of the upstream tools, I would suggest changing the extension to txt.gz and noting the format in https://github.com/AlexsLemonade/OpenPBTA-analysis#data-formats + possibly adding something to docs/format.

Someone noted in person that changing the files with something like sed /s/ /\t/g might be preferrable for downstream analysts over making the documentation change.

Evaluate concordance and determine overall set of high-confidence somatic variants for PBTA.

This should probably be partitioned by cancer type and/or molecular subtype. This could be a stacked bar chart (one bar per fusion) for instance with colors representing different cancer types.

With the additions in #45, all the data files will need to be represented in the md5sum.txt file. This issue pertains to rsem.fpkm.rds, first mentioned here: https://github.com/AlexsLemonade/OpenPBTA-analysis/pull/43/files#diff-04c6e90faac2675aa89e2176d2eec7d8R125

This will need to be added to md5sum.txt.

Summary of druggable/actionable genes via mutation, over-expression, fusion, synthetic lethality and actionable features i.e. TMB, mutational signatures, PD-L1 expression, etc..across the PBTA. Analysis of what's currently used and what may be some potential opportunities.

This figure typically focuses on a set of high frequency and cancer associated lesions (mutations / CNV) and will typically have corresponding phenotypic data so as to visualize associations. See Figure 1 (https://www.biorxiv.org/content/10.1101/656587v1.full) for an example.

We have planned to include an analysis of tumor mutation burden as well as a comparison to the tumor mutation burden of adult tumors in TCGA. This issue is for fleshing out the analysis. The manuscript issue for the methods description is: AlexsLemonade/OpenPBTA-manuscript#7

It is often helpful to have a part of the first figure for a dataset landscape paper that shows how the samples are distributed across cancer type to characterize the overall dataset. We would like a figure that summarizes the content of the dataset "at a glance."

File(s)

What data file(s) does this issue pertain to?

Release

V2/V3

What release are you using?

Link to OpenPBTA-manuscript

Put a link to the relevant section of the OpenPBTA manuscript here.

Question/issue

Put your question or report your issue here.
BS_6DCSD5Y6 is duplicated for its RNA-Seq entry in the clinical file RNA-Seq and this is due to two different ages at diagnosis. Will look into which is correct - we may have a discrepancy between KF DRC and Pedcbio portal. Will also fix this with V4. cc: @yuankunzhu

File(s)

All files

Release

release-v3-20190829

Question/issue

I was attempting to increase the number of matched participants in the subset file to 100 over on #91. I map the biospecimen IDs in each file to the participant IDs:

I find 91 participants in all files

Here are the number of mapped participant ids for each file:

> lapply(participant_id_list, length)
$cnvkit
[1] 939

$controlfreec
[1] 941

$arriba
[1] 1029

$starfusion
[1] 720

$kallisto
[1] 1029

$rsem
[1] 944

$mutect2
[1] 973

$strelka2
[1] 973

$manta
[1] 940

$lumpy
[1] 147

I feel like the number of participants in all was higher in release v2 but I don't remember off the top of my head. Alternatively, this strategy for finding matched participants should be changed.

File(s)

README.md
release-notes.md

Release

release-v4-20190909

Question/issue

This file does not seem to be present in the S3 bucket at the location specified in download-data.sh. Note that the script refers to release-notes.md, but the release notes suggest it may be named README.md. Neither seems to be present, however.

File(s)

What data file(s) does this issue pertain to?

Release

What release are you using?

Link to OpenPBTA-manuscript

Put a link to the relevant section of the OpenPBTA manuscript here.

Question/issue

Put your question or report your issue here.
We are currently generating lancet and vardict calls and will release in V5

We should flush out this process - giving people access to CAVATICA to analyze data, ability to use virtual notebooks to do analyses, then submission of their code/figures/notebooks to github.

It would be really handy if there was a script that could be run over the input data files to produce files containing data for only a small, specified number of matched participants. I'm imagining that we would use this to produce data files with information for 10 or so participants. We would then use this code via CircleCI to test the analytical code that people contribute without having to run over the full set of input data.

Without methylation data, the goal of this analysis is to best molecularly subtype all brain tumors based on defining features.

• Medulloblastomas (MBs) into SHH, WNT, Group 3, and Group 4 - @PichaiRaman has created a classifier here and classified the samples with RNA-Seq.
• ATRTs into SHH, MYC, TYR. See: 1-s2.0-S1535610816300356-main.pdf
• Ependymomas into ST-RELA, ST-YAP1, PF-A, and PF-B. Will require some genome-wide copy number examination. See: 0D556E30-D832-4D6F-9F2F-49554F68D56F.pdf
• High-grade gliomas, according to brain region (if known, eg: midline vs hemispheric) and defining lesion (eg H3F3A K27/28M, G34/5R/V, HIST1H3B K27/28M, IDH-mutant, ACVR1-mutant, PDGFRA, and mesenchymal subtypes). See: WHO 2016 classifications
• Non-MB embryonal tumors (old PNET classified tumors) into ETMRs and CNS embryonal NOS
• Oligodendrogliomas as IDH-mutant and 1p/19q codeleted
• Chordomas - subtype into poorly differentiated (INI1/SMARCB1 loss)

CNV

1. We have ControlFreeC and CNVkit. CNVkit output is VCF – do users want raw output for these or should we convert ControlFreeC output to VCF, then annotate all with AnnotSV (which is what we are using for SV annotation) and then provide annotated tab files?

SV
2. Do we want Novobreak as a 3rd algorithm? Have Manta and and will have Lumpy.

Annotated/Merge Files
3. We ran Strelka2 and Mutect2 for somatic calls. Do we want to use VEP for SNVs (MAFs are currently annotated with this) or SNPeff? We have both, but VEP was provided.
4. Do we want to provide a limited merged tab file in maf format that includes both Mutect2 and Strelka2 variant calls with a column for “Algorithm”? If so, we need to define the fields - Is this better for users?
5. Should we do the same for CNVs and SVs – provide one merge file of all data? Note: these will be large – (22GB +) and won’t be able to be processed locally
6. Should we set up a CAVATICA project for people to work on interactive analyses (R studio/Jupyter notebooks) or do we think they will do that on their own? We can probably provide some cloud credits for this if they want to do on own – maybe we need instructions?
7. Current plan is for CBTTC and PNOC to be separate merge files, but so we want these to remain separate or be merged? If separate, we need instructions for downstream users to merge all files they download, so I prefer merging all from the start.

SNV
8. Do we want only tumor/normal paired samples or should we also include tumor-only samples? If yes to tumor only, what should we use as PON - normal from cross-brain tumors, normal within disease type (if high enough N), non-cancer pediatric normals?

@gonzolgarcia @samesense @Yiran-Guo @gwaygenomics @jpfeil @kgaonkar6 @LauraEgolf @apexamodi @nathankendsersky @afarrel @jaclyn-taroni @PichaiRaman please weigh in. cc: @yuankunzhu

Currently planned for 24-Sept-2019

Planned addition + changes from @jharenza :

• Vardict - VEP annotated MAF (related: #103)
• Lancet – VEP annotated MAF (related: #103)
• Updated arriba fusion file (related: #92 (comment))
• RSEM gene level count matrix
• RSEM transcript level count matrix (related: #14 (comment))
• Updated lumpy file - to fix T/N columns (related: #27 (comment))

somehow the .gitignore file keeps ignoring changes in data/ path, i saw you @cgreene only skips the DS_Store file, not sure what happened. I ended up using git add -f for force adding my commits.

Comparison of genomic/transcriptomic features and derivates such as TMB, pathway activity, etc... to survival. Determine markers that could add prognostic values to different PBTA cancer types. KM Plots for significant and interesting features.

We plan to produce a MultiPLIER LV matrix. We then plan to use that analysis to identify LVs associated with survival as well as those associated with samples that contain alterations in selected genes. This issue is for more fully fleshing out the analysis.

File(s)

What data file(s) does this issue pertain to?
pbta-gene-expression-rsem.fpkm.RDS

Release

What release are you using?
V2/V3

Link to OpenPBTA-manuscript

Put a link to the relevant section of the OpenPBTA manuscript here.

Question/issue

Put your question or report your issue here.
I noticed that there are 1123 BS_IDs in this file, but our data has 1028 unique RNA BS_IDs. Additionally, BS_ID BS_XM1AHBDJ is missing from the RSEM file, but present in the kallisto file. We will have to update this in V4.

Determine differential isoforms across PBTA or within specific PBTA cancers and correlate back to genomic or phenotypic characteristics.

Scientific goals

What are the scientific goals of the analysis?
Accurately predict sex of PBTA samples for downstream analyses

Proposed methods

What methods do you plan to use to accomplish the scientific goals?
Ratio of Y:X+Y chromosomes
XIST expression in RNA

Required input data

What input data will you use for this analysis?
normal DNA BAMS
RNA-Seq FPKM

Proposed timeline

What is the timeline for the analysis?
One week

Relevant literature

If there is relevant scientific literature, put links to those items here.

Here, we will filter potential artifacts, filter fusions observed in normal tissue, retain high-confidence calls, and annotate with several databases to create a final list of putative driver fusions.

Scientific goals

Identify potential variants of unknown significance that may play a role in tumorigenesis.

Proposed methods

Look at frequently recurrent VUS, determine if there is a functional effect on a gene (or it is in a regulatory region), look at VAF and focus on samples without other relevant driver mutations.

Required input data

VCF or MAF files

Proposed timeline

2 weeks

Relevant literature

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6359859/

We should put together a script to download the data into a defined folder. We should require folks to not modify that folder. We should include dummy germline data in the download if we are not able to distribute it without restrictions so that at least the CI works.

File(s)

pbta-histologies.tsv and to a lesser extent, pbta-snv-strelka2.vep.maf.gz

Release

release-v2-20190809

Link to OpenPBTA-manuscript

The methods section about data generation is where I would have expected to see information about the "Panel" samples.

Question/issue

What does it mean in pbta-histologies.tsv when an experimental strategy is "Panel"?

For the entire set of metadata, there appears to be only two samples that are called "Panel". Here is their information:

More specifically, there is only one tumor sample in the Strelka2 data that is noted as using this strategy: BS_7KR13R3P. The rest of the samples in the Strelka2 data are noted as "WGS" and "WES" which is what I would have expected for all the samples.

What does "Panel" mean regarding these samples? Is it incorrectly labeled? If the experimental strategy for these two samples is actually different, should it be discarded from the data in general?

Compare primary samples to relapse samples to determine genomic and transcriptomic differences and identify trends. See (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5172533/).

Problem

13 samples that have data in the MAF files (pbta-snv-mutect2.vep.maf.gz and pbta-snv-strelka2.vep.maf.gz) are missing from the metadata. I am using data release-v2-20190809 and downloaded it via the bash script this morning.

Samples that are missing metadata but have data reported in the MAF files.

"BS_0BA5TZND" 
"BS_0HYD1VHH" 
"BS_0PQGSCJA" 
"BS_3GSGMV4T" 
"BS_3RJ3CDE6" 
"BS_AH51BZ8D" 
"BS_BF95AEXC" 
"BS_D0T6V861" 
"BS_KVPJVJR7" 
"BS_TC8R5HY4" 
"BS_TKF4H7CZ" 
"BS_V1DPR5TD"
"BS_X2G3JMM1"

Here are the steps I took to obtain this list of samples were missing from the metadata. Let me know if there is something I missing or misunderstanding.

# Read in the metadata
metadata <- readr::read_tsv(file.path("..", "..", "data", "pbta-histologies.tsv"))

# Read in Mutect2 data
mutect2 <- maftools::read.maf(file.path("..", "..", "data", "pbta-snv-mutect2.vep.maf.gz"))

# Read in Strelka2 data
strelka <- maftools::read.maf(file.path("..", "..", "data", "pbta-snv-strelka2.vep.maf.gz"))

# Find out what samples don't have data in the metadata
missing_mutect2_samples <- mutect2@data %>% 
  dplyr::filter(!(Tumor_Sample_Barcode %in% metadata$Kids_First_Biospecimen_ID)) %>%
  dplyr::distinct(Tumor_Sample_Barcode, .keep_all = TRUE) %>%
  dplyr::pull(Tumor_Sample_Barcode)

missing_mutect2_samples <- strelka2@data %>% 
  dplyr::filter(!(Tumor_Sample_Barcode %in% metadata$Kids_First_Biospecimen_ID)) %>%
  dplyr::distinct(Tumor_Sample_Barcode, .keep_all = TRUE) %>%
  dplyr::pull(Tumor_Sample_Barcode)

# These are the same 13 samples in both datasets
sort(intersect(missing_mutect2_samples, missing_strelka2_samples))

Next steps:

Were these samples supposed to be removed from the MAF files? (And if so, why?)
OR
Where is the metadata for these 13 samples and how do we get it?

The goal of this analysis is to identify recurring, putative oncogenic SV lesions/genes with recurrent lesions within and/or across brain tumor subtypes. Coupled to this, investigation of chromothripsis and correlations with survival should be done. See this paper for chromothripsis analysis using WGS.