Comments (3)
Hi Maria,
Thanks for the follow up!
Just to give you a little more background as I'm sure you're aware from reading the paper (https://genome.cshlp.org/content/early/2018/06/13/gr.228080.117), the expression-based model in HoneyBADGER is able to identify large-scale copy-number alterations based on consistent deviations in expression from a normal reference. For example, if your cancer sample has a duplication of a chromosome arm, on average, you should see increased expression for the genes within that chromosome arm compared to the normal tissue. Likewise, for a deletion of a chromosome arm, you should see decreased expression compared to the normal tissue. And HoneyBADGER uses a hidden Markov model to identify these potential regions as well as a Bayesian hierarchical approach to quantify the probability of these alterations within individual cells.
Now what this means practically is that we need single-cell RNA-seq cancer data (in your case pancreatic adenocarcinoma) and a normal RNA-seq analogue (in your case normal pancreas). And we need to normalize these two dataset in as much as the same way as possible such that we can interpret deviations in the cancer expression data as signal from underlying CNV changes rather than simply artifacts from normalization differences. In the manuscript, we used a Counts Per Million (CPM) type normalization. You may also use other normalizations as long as they are consistent.
To have HoneyBADGER acknowledge your own files, you will just need to pass them in:
hb <- new('HoneyBADGER', name='pancAdeno')
hb$setGexpMats(YOUR_SCRNASEQ_CPM_NORMALIZED_DATA,
YOUR_NORMALSEQ_CPM_NORMALIZED_DATA,
mart.obj, filter=TRUE, scale=TRUE, verbose=TRUE)
Here, YOUR_SCRNASEQ_CPM_NORMALIZED_DATA
would be your CPM normalized single-cell RNA-seq data with genes as rows, cells as columns. And likewise YOUR_NORMALSEQ_CPM_NORMALIZED_DATA
would be your CPM normalized bulk RNA-seq GTEX data with genes as rows, samples as columns.
Note that unlike the tutorial, here, we are using your normalized data matrices AND we are setting the filter
and scale
parameters to TRUE
. The built-in datasets have already been pre-normalized, filtered, and scaled.
There are additional filter parameters that you will likely need to toggle. For example, it is well known that single-cell RNA-seq and bulk RNA-seq data differ systematically such that simply summing up expression from single-cell RNA-seq data does not yield the same expression quantifications as bulk RNA-seq data, especially for lowly expressed genes due to technical limitations. Therefore, if we see a gene that's lowly expressed in your normal reference's bulk RNA-seq and not at all expressed in cancer single-cell RNA-seq, this is likely due to technical limitations rather than because there is a deletion in the cancer single-cell RNA-seq. So there are a few parameters minMeanBoth
minMeanTest
and minMeanRef
that simply attempt to filter out such poorly detected genes. See ?setGexpMats
for more information.
In the manuscript, for one 10X multiple myeloma dataset we actually had single-cell RNA-seq data for both the cancer and the normal reference, sequenced on the same platform using the same chemistry, etc. In such a case, this kind of filtering was no longer necessary. So do keep in mind what is the purpose of each step in the tutorial and see if it is still appropriate for your analysis.
Hope that helps. Let me know if you have any additional questions,
Jean
from honeybadger.
Hi Jean!
Thank you SO much for the help here! Really appreciate it!
One upstanding issue- I have normalized 2 data sets (hpne_means is normal pancreas tissue, bnc2 is my sample data), and went to run the hb$setGexpMats command to generate the gene expression profile plot. This is the error I received:
hb$setGexpMats(bnc2, hpne_means, mart.obj, filter=TRUE, scale=TRUE, verbose=TRUE)
Initializing expression matrices ...
23 genes passed filtering ...
Scaling coverage ...
Normalizing gene expression for 23 genes and 400 cells ...
Done setting initial expression matrices!
hb$plotGexpProfile()
Error in seq_len(nrow(d)) :
argument must be coercible to non-negative integer
In addition: Warning message:
In seq_len(nrow(d)) : first element used of 'length.out' argument
I have not gotten this error before, even as I was troubleshooting HoneyBADGER on various data sets over the past couple weeks, and I'm not sure what I can do to solve.
Any help??
Thanks,
Maria
from honeybadger.
Hi Maria,
Great! Sounds like you're making progress.
So as I mentioned previously, there are a few threshold parameters minMeanBoth
minMeanTest
and minMeanRef
that attempt filter out such poorly detected genes. You can imagine how these threshold parameters may need to be tweaked depending on your dataset. Perhaps your sequencing depth isn't as deep as the dataset I used to set these default parameters.
It looks like there are only 23 genes that remained after filtering! This is a very very small number of genes so most likely the default filter thresholds are way too stringent for your data.
You should look at the distribution of gene expression in your data. For example:
# histogram of average gene expression
hist(rowMeans(bnc2))
# you may want to plot on a log scale
hist(log10(rowMeans(bnc2)+1)
As well as how well these expression measurements correlate with your bulk. For example:
# correlation of average gene expression
plot(hpne_means, rowMeans(bnc2))
# you may want to plot on a log scale
plot(log10(hpne_means+1), log10(rowMeans(bnc2)+1))
# add lines for your thresholds
abline(v = minMeanTest, col='red')
abline(h = minMeanRef, col='red')
abline(v = minMeanBoth, col='blue')
abline(h = minMeanBoth, col='blue')
My sense is there will be some systematic differences, particularly for lowly expressed genes due to drop-out artifacts in your single-cell data. You can use these plots to choose a more sensible threshold for your data.
Best,
Jean
from honeybadger.
Related Issues (20)
- lt$setGexpDev HOT 1
- Warning in setGeneFactors and error in retestIdentifiedCnvs
- Allele-mode for selecting normal cells HOT 2
- 10X + Honeybadger question HOT 2
- Error in summarizeResults HOT 3
- Error in calcGexpCnvBoundaries when running with numeric chromosome names HOT 1
- Error: subscript contains invalid names HOT 9
- read bam files HOT 4
- Filtering of identified CNVs HOT 2
- speed of running setAlleleMats step
- What is the reason for only including snps in HoneyBADGER? HOT 2
- gene filtering different in HoneyBADGER object HOT 3
- Error: Failed to install 'HoneyBADGER' from GitHub HOT 2
- Showing error when trying Getting_Started.Rmd
- no method for coercing this S4 class to a vector HOT 3
- showing NULL in the step of calcGexpCnvBoundaries - getting started tutorial HOT 12
- Applying bayesian hierarchical model in integrating analyses tutorial
- error in hb$summarizeResults
- Error in summarizeResults for the allele-based approach
- Error in calcCombCnvProb HOT 1
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from honeybadger.