17

Benefits of central repository for Community Having a central repository for packages is very useful. For couple of reasons: It makes very easy to resolve dependencies. Installing all the dependencies manually would be exhausting but also dangerous (point 2). Package compatibility! If I install package with dependencies, I would like to be sure that I ...


15

As suggested, here’s an example showing the relevant lines from a DESCRIPTION file from a CRAN/GitHub hosted project that has Bioconductor dependencies (truncated): Depends: R (>= 3.3.0) biocViews: Imports: methods, snpStats, dplyr The relevant bit is the empty biocViews: declaration, which allows the Bioconductor dependency {snpStats} ...


14

The simplest manner is to not use a wald test, but rather an LRT with a reduced model lacking the factor of interest: dds = DESeq(dds, test="LRT" reduced=~geno+geno:Treatment) The above would give you results for Treatment regardless of level while still accounting for a possible interaction (i.e., a "main effect of treatment, regardless of the type of ...


12

The normalized counts themselves can be accessed with counts(dds, normalized=T). Now as to what the baseMean actually means, that will depend upon whether an "expanded model matrix" is in use or not. Given your previous question, we can see that geno_treat has a bunch of levels, which means that expanded models are not in use. In such cases, the baseMean ...


9

Here is a list of the advantages of having Bioconductor for the bioinformatic community: Outreach: You have a repository for the field, in that language. Some packages related to bioinformatics (in R) are distributed through personal repositories, CRAN, github, bitbucket, sourceforge, but they are less used and harder to find. There are such efforts in ...


9

Save the different scripts with git (seems overkill) Whoa. I did an actual double take when reading this:1 it’s the opposite of overkill. Version controlling your scripts (using Git or something similar) is the absolute minimum, and should become completely automatic. For every new project I begin, one of the very first steps is to issue the git init ...


9

There's a trick to this where one needs to add biocViews: to the package Description. That's the only solution I've ever seen to allowing automatic installation of bioconductor dependencies. If you need a couple of examples, then click through the link I posted and scroll down to pull requests referencing that issue, they will generally include the actual ...


8

According to the manual, all you need to do is: library('GEOquery') gseGSE16146 <- getGEO('GSE16146', GSEMatrix=FALSE) As explanation, getGEO() outputs by default to GSEMatrix=TRUE and returns a list of ExpressionSet objects. You should get what you were looking for with: Table(GSMList(gseGSE16146)[[1]])[1:5,] The manual has also a paragraph about ...


8

You need to specify the number without the version. Instead of "ENSMUST00000178862.1" just "ENSMUST00000178862": You can do this with one more line: g <- gsub("\\..*", "", rownames(txi.kallisto$counts)) (hgnc_symbols <- getBM(attributes = c("hgnc_symbol", "chromosome_name", "ensembl_transcript_id"), filters = "ensembl_transcript_id", values = g, mart ...


8

It seems that the "combining factors" trick described in part 3.3 of DESeq2 current "vignette" (as of may 2017) under the title "Interaction" is a way to access to the desired contrasts. It seems possible do do it directly when building the colData and when calling DESeqDataSetFromMatrix: Let's add a combined "geno" and "treat" factors to the future ...


8

Did you try the fill argument? Something like this: track <- AnnotationTrack(start=c(1,5,7), end=c(2,6,10), strand=c('*','*','*'), stacking="dense", showFeatureId=TRUE, id=c('red','blue', 'red'), fill=c('red','blue', 'red'))


7

According to the FGSEA preprint: We ran reference GSEA with default parameters. The permutation number was set to 1000, which means that for each input gene set 1000 independent samples were generated. The run took 100 seconds and resulted in 79 gene sets with GSEA-adjusted FDR q-value of less than 10−2. All significant gene sets were in a ...


7

There is no golden/standard way to define a gene signature but they are completely different from a gene set enrichment analysis (GSEA). I will start with how to obtain a gene signature: Usually this kind of signatures are defined by comparing a group against another ie doing a differential gene expression analysis. One then selects the most up or down ...


6

It depends what you mean by “normalised”. As Devon said, the normalized = TRUE argument to the count function gives you normalised counts. However, these are “only” library-size normalised (i.e. divided by the sizeFactors(dds)). However, as the vignette explains, downstream processing generally requires more advanced normalisation, to account for the ...


6

I don't believe this is possible using biomaRt, nor using AnnotationHub. I have two suggestions, neither of them very satisfactory. First, you can specify an Ensembl archive for biomaRt, for example: mart72.hs <- useMart("ENSEMBL_MART_ENSEMBL", "hsapiens_gene_ensembl", host = "jun2013.archive.ensembl.org") Of course, that ...


6

Methylation levels have high local correlation, so Fisher's method would be problematic. Having said that, you have no reason to use Fisher's method after a paired t-test. A paired t-test will give you a single p-value per gene, which is what you want. Do be sure to only include CpG with some minimal coverage in both group.


6

It looks like you were using an old annotation. The problematic IDs you posted existed in the GRCh37 annotations, but don't in the most recent GRCh38 annotation. For that reason they were excluded. The IDs that have - as symbols don't have associated symbols, but are present in the database. To use an archived version in biomart: mart = useDataset("...


6

subsetter = function(gr, cname) { return(mcols(gr)[[cname]]) } You can then use things like subsetter(gr, "GC") and subsetter(gr, "score").


6

According to the documentation (?Biostrings::DNAStringSet): width(x): A vector of non-negative integers containing the number of letters for each element in x. Note that width(x) is also defined for a character vector with no NAs and is equivalent to nchar(x, type="bytes"). names(x): NULL or a character vector of ...


6

You can't cbind a bunch of obscure object types. If you want merged count tables you should do this: mdat <- do.call(cbind,lapply(dat,assay)) Where row.names are Ensembl gene IDs and col.names are the SRR accessions. Then run your table writing command. - If you want the coordinates of your genes then do this to make a bed with the genomic locations:...


5

The main purpose of git is to version code, which usually means sequential improvement of the codebase. While it is possible to use branches for multiple variants of the software, permanent branches are traditionally used for gradual integration of new features (i.e. dev/testing/master branches). Supporting multiple independent branches requires some ...


5

I don't think that the issue is the low counts, but rather the number of features without any real variance (the black dots at the bottom). So what the heck is the dispersion plot and why does one need to fit it anyway? In a typical RNAseq experiment, one measures many thousands of genes with only a few replicates per biological group. This then leads to ...


5

You see negative values with your function because you're setting the average of each row to 0 and its standard deviation to 1. In general, I would trust a standard normalization method (rma in this case) more than some random "truncate and then scale the rows" method. Your method isn't even doing any between-array normalization, which is the benefit of rma....


5

There is no one-to-one mapping of gene ids from one database to the other. Ensembl (who maintain Ensembl ENSG IDs), ncbi (who maintain EntrezGene IDs and RefSeq transcript ids) and HUGO (who maintain gene symbols/names) all have different ideas about what a gene is, which part of the genome belongs to which gene etc. In fact, it's worse than that because ...


5

Let's look into this a bit more deeply. For instance: HUGO: SOGA3 Ensembl 1: ENSG00000214338 Ensembl 2: ENSG00000255330 The Ensembl pages (linked above) for both ENSG00000214338 and ENSG00000255330 show that these genes have 3 paralogs and each lists the other as one of the paralogs. Your other two examples are similar, they are also paralogs of each other....


5

As far as I'm aware, Illumina provide CSV annotation files for all their sequencing chips, which can be used when they can't be found in Bioconductor. You can find annotation information for the PorcineSNP60 here, in particular the Manifest file (CSV format). The format is Illumina's weird "we say it's a CSV because there are commas in it" format, so if ...


5

Your request is easy to implement. I wouldn't use any libraries in this case. As you haven't showed your implementation, I will provide one: def gen_cigar(ref, qry): if len(ref) != len(qry): raise Exception('unequal length') cigar = [] for i in range(len(ref)): r, q = ref[i], qry[i]; if r is '-' and q is '-': ...


5

In such cases it is best to check the str of an object gr. Then we could see that meta data is just a dataframe inside S4 object: gr@elementMetadata # DataFrame with 10 rows and 2 columns # score GC # <integer> <numeric> # 1 1 1.0000000 # 2 2 0.8888889 # ... Once we know how to access the dataframe, we can ...


5

You can submit a package to a repo after the publication, but I would say that you really should do it before. Bioconductor accept only packages that are not published on CRAN, however, academic publications are fine. Conversely, Bioinformatics does not enforce any particular platform for sharing the code, but you have to make it available somehow (e.g. ...


5

Ah, looks like I can't even procrastinate on StackExchange anymore without seeing work-related stuff. Oh well. Anyway, the other answers and comments are way off. scran has supported sparse matrices for years, ever since we switched over to the SingleCellExperiment class as our basic data structure. quickCluster does no coercion to dense format unless you ...


Only top voted, non community-wiki answers of a minimum length are eligible