text:
Handling sample identity in aggregated 10x libraries?
cellranger aggr can combine multiple libraries (samples), and appends each barcode with an integer (e.g. AGACCATTGAGACTTA-1). The sample identity is not recorded in the combined matrix.mtx file.
https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/aggregate#gem_groups
How can we keep and assign sample information to each cell after loading the data into R? (e.g. Seurat::Read10X())

An example solution for Seurat:
Retrieve sample IDs from the .csv used with cellranger:
samples <- read.csv(file.path("/path/to/csv", "nameof.csv"), stringsAsFactors=F)

Load the 10x dataset and initialize the Seurat object:
cells.data <- Read10X("path/to/filtered_gene_bc_matrices")

cells <- new("seurat", raw.data=cells.data)

Get barcodes and suffix:
cellcodes <- as.data.frame(cells@raw.data@Dimnames[[2]])
colnames(cellcodes) <- "barcodes"
rownames(cellcodes) <- cellcodes$barcodes

cellcodes$libcodes <- as.factor(gsub(pattern=".+-", replacement="", cellcodes$barcodes))
cellcodes$samples <- as.vector(samples$library_id[cellcodes$libcodes])

Create dataframe for meta.data argument and set up object:
sampleidentity <- cellcodes["samples"]

cells <- Setup(cells,
               meta.data=sampleidentity,
               min.cells=3,
               min.genes=200,
               do.logNormalize=T, total.expr=1e4, project="projectname")

Use the group.by argument of various functions.

The new Seurat 2.0 has the MergeSeurat() and AddSamples() functions. See this question and seurat vignette

You have to keep track of the file order you used for aggr.  The suffix number represents which represents the aggregated sample. You could store this information in a text file and load it into R independently of the Read10X function and combine it with the cell names to get their sample names. From the page you link:

Gem Groups
This is an integer that is appended to each barcode in the
  gene-barcode matrix. For example, AGACCATTGAGACTTA-1 and
  AGACCATTGAGACTTA-2 are distinct cell barcodes from different
  libraries, despite having the same nucleotide sequence.
The numbering of the GEM groups will reflect the order that the
  libraries were provided in the Aggregation CSV.


text:
What are the ways to process a list of differentially expressed genes?
We are studying six different human macrophage/dendritic cell types isolated from healthy skin. They all differ from each other in a few cell surface markers.
We are interested in the characteristics of each cell type ("marker" genes or biological processes), the differences between them (especially cell surface proteins), and their relations (e.g. "which are the closest relatives?") as can be inferred from the transcriptome, from an immunological point of view. The wider context is HIV-infection, thus HIV infection-related differences (and similarities) are of particular interest.
One naive approach is to contrast two of the most similar cell types (as inferred from previous knowledge, the sorting strategy and the PCA plot), and produce a list of differentially expressed (DE) genes.
I have now a list of ~500 DE genes between two cell types, together with fold changes, expression means etc., produced with DESeq2 on bulk RNA-seq data.
What are the best practices to process this list to answer some of the above?
I've found PPI analysis on NetworkAnalyst, and  GO and pathway analysis on InnateDB useful.
What other options are available?
I'm particularly interested in approaches using R.

The question you should be asking isn't "what can we do with this?", but rather "what biological question do we want to answer?" You sound like you're aimlessly fishing for anything, which is likely not worth your time.

What do you mean by "process"? Process the list in order to do what? What question are you asking? Please [edit] your post and clarify what you are trying to do. As it stands, a valid answer would be that you can use wc to count the number of lines in the file. That is a type of "processing" but, obviously, not very helpful for you. The more specific you make your question, the higher the chances of your getting a useful answer.

Hi Peter, thanks for your question and welcome to Bioinformatics Stack Exchange. Bioinformatics is a very large subject area, frequently with many different solutions to the same problem. It would be helpful if you could give more story / context surrounding your problem to guide answers. For example, what organism are you working with? Are these conditions treatment vs control, or separate treatments (without a control)? What kind of results are you interested in getting out of this analysis?

Thank you all for the feedback. I understand that each dataset has to be studied according to the questions of interest, just wanted to get a few ideas. Also tried to make the question general so that it will be useful for more people, but I now added the specifics.

I think the question is still too broad. What knowledge do you expect to learn? Could you narrow more to one of your interests (Process, pathways, or relations between cells lines)?  Which is your "immunological point of view" ?

It's best to stick to questions that are as specific as possible, and let other people decide on their own generalisations about the answers. There are so many different things that can be done in bioinformatics even for a single task that it's no useful to speculate about what other similar things could be done.

You originally had asked a very broad question, so I'll try to demonstrate why that is such a hard question to answer. I've done two fairly large differential analysis studies (and a few smaller ones) covering very different areas of research, and the approaches that other researchers used subsequent to my differential expression calculations were unsurprisingly also very different.
In the first study, I did a de-novo transcriptome assembly of Schmidtea mediterranea (Smed) from Illumina RNASeq reads, then carried out a DE analysis (DESeq) on a β-catenin knockout mutant population compared to wildtype. The other researchers had some very impressive tools available at their disposal, allowing them to subsequently visualise the positional expression of differentially-expressed genes in the knockout and wildtype populations after amputation. They found a few candidate genes that had graded expression along the anterior-posterior axis in the wildtype population, and this expression disappeared in the knockout mutants. One of these genes was a teashirt (tsh) family member, and a good candidate for further discovery based on the researchers' existing knowledge of the gene pathways associated with cellular regeneration. To cut a very long story short, they experimentally validated the teashirt link through further RNAi studies, generating a double-headed Smed as a result. This was followed up in Zebrafish (Danio rerio), because there happened to be another building about 100m away that specialised in Zebrafish development.
The second study was/is a bit more similar to your current situation. We looked at sorted dendritic cell expression in mice after the application of two Th2-inducing treatment conditions (injected Nippostrongylus brasiliensis; applied DBP+FITC) and one control condition (injected PBS). I did a differential expression analysis with DESeq2, and developed a Shiny App for result exploration after I got tired of getting asked to do a whole bunch of mechanical follow-up tasks (mostly generating graphs, lists, and heatmaps for different gene subsets). We ended up with thousands of differentially expressed genes, so it wasn't feasible to look individually at every gene for research evidence and pathways -- one scientist spent a couple of weeks looking at one interesting-looking candidate gene. The differentially-expressed genes were fed into the oh-so-expensive IPA (which gives results, but not necessarily good ones) and DAVID, among a few others. These helped a little bit, and led to what was the central idea of our paper (i.e. that there were at least two different Th2 responses), but we still had the problem of too many "significant" pathways, and lots of differentially-expressed genes that were members of a whole bunch of other pathways. We tried WGCNA; we did a bit of GSEA; and we compared Limma and DESeq2, but still couldn't find anything concrete or obvious to demonstrate what was triggering the Th2 response. There were experimental studies as well: blocking and mutant studies to investigate the two types of Th2 response (and experimentally confirm the differential expression). Our experiments and exploration are still ongoing, but we felt it was necessary to get one paper out at the intermediate stage, so just ended up taking the top 25 differentially-expressed genes in each experiment (by fold change), and finding some way of grouping some of them together (via expert knowledge and literature searches). It kind-of-sort-of worked, but we're left with a bit of frustration in that we (and the tools we used) were basically stumped by having too many biologically-relevant results.

We should move to a chat,... but study 2 looks like a study I did: how did you make your enrichment analysis (outside IPA), specially of GO? How did you try WGCNA? Did you do a GSEA of the genes up-regulated in the comparison of a treatment vs control in the comparison of treatment1 vs treatment2?

WCGNA: We had about 80 different biological replicates comprising about 25 different treatment conditions, so there was enough expression variation to be able to pull out a few shared gene sets

GO enrichment analysis wasn't done; I don't know of any GOEA tools that use the hierarchical structure of the ontology map, and without that it's very prone to selection bias. For GSEA, we were mostly looking at treatment vs control (for various treatments and cell populations), and I think using the gene sets in the database that seemed to be immune-related. I'm not completely sure on that; it was other scientists that did the GSEA.

For the WGCNA, I first did a VST transform of the expression data in DESeq, normalised by the length of the longest transcript for that gene, then determined the 5,000 genes that were most variable (by MAD) across all our samples. This matrix of "lowFiltered" genes was fed into WGCNA by another scientist, and we spent a fair amount of time trying to understand the process and get the results looking right.

topGO package of Bioconductor uses the DAG structure of GO to find the enriched terms. WGCNA requires homogeneous samples, so if you mixed the 25 treatment conditions in those 80 samples to build a network, it would not detect the network for each treatment. You would need to do a consensus analysis (with multiset.* functions) to retrieve the network shared between each treatment, but WGCNA requires above 15 samples per condition to effectively work.

Thanks, I might give topGO a try once I've sorted out the proportionality calculations

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

sort: votes


text:
How to automate NCBI genome download
I need to download all the completely assembled cyanobacterial genome's GenBank file(.gbff) from NCBI(RefSeq or INSDC ftp data).
For this I think, the steps are:

Need to find the completely assembled genomes.
find the GenBank file URL based on the taxonomic name.
download the GenBank file(.gbff file)

Is their any way to do this by using any python module or any other idea??

Biopython has tools that could help you for steps 2 and 3, but I'm not sure how to get the information of your step 1.

can you please tell the tool name for step 2.

@Arijit Searching for "Download biopython Genebank" lead you here

Outline of solution:

get this file: ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/bacteria/assembly_summary.txt
filter for lines having "Complete Genome" in column 12
filter for lines having a taxid (column 6) that is a descendant of taxon id 1117 (phylum Cyanobacteria), you can use nodes.dmp from ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz
The file download URL is the concatenation of column 20 and the last field of column 20, after separating it at '/', followed by appending '_genomic.gbff.gz'.  For example: ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/007/365/GCF_000007365.1_ASM736v1/GCF_000007365.1_ASM736v1_genomic.gbff.gz

Edit:
Since I will need something like that soon, I made a Perl script for downloading genomes.
For Cyanobacteria, you would do ./download-refseq-genomes.pl 1117

As an alternative, you can use the NCBI Datasets CLI, in a two step process to generate a dehydrated dataset for later hydration:
# Download metadata for all Cyanobacteria
# And pull out accession for 'Complete' genomes
./datasets assembly-descriptors taxon Cyanobacteria --limit ALL \
     | jq -r '.assemblies[].assembly | "\(.assembly_accession) \(.assembly_level)"' \
     | grep Complete \
     | cut -f 1 -d ' ' > accs.txt

# Download a dehydrated dataset for the provided set of accessions
./datasets download assembly -i accs.txt -f 1117_complete.zip \
     --dehydrated \
     --exclude-protein \
     --exclude-rna \
     --exclude-seq \
     --exclude-gff3 \
     --include-gbff

You may be interested in some of the other files that are excluded.
Once you have the dehydrated dataset, you can rehydrate it when you are ready to analyze the data with:
unzip 1117_complete.zip -d 1117_complete
./datasets rehydrate -f 1117_complete

Use the NCBI Datasets command line tool to download genomes by taxon (NCBI Taxonomy ID, scientific or common name at any tax rank).
datasets download genome taxon cyanobacterial \
  --assembly-source refseq --assembly-level complete_genome \
  --exclude-genomic-cds --exclude-gff3 --exclude-protein \
  --exclude-rna --exclude-seq --include-gbff \
  --filename genomes.zip

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Does GISTIC (v 2.0) estimate amplified/deleted probabilities on a single sample basis?
Does GISTIC 2.0 estimate the background model:
G = -log(Probability | Background)
by permuting within the sample or across all samples in the set?
The paper describes the probabilistic scoring method based on permutations, but I could not understand if this permutation is performed within the sample only. The documentation page seems to suggest across samples, but this would mean that different sized sets might lead to different outcomes on the sample sample.
Basically, does it matter if the set is constituted by, for example, 10 samples of the same tissue of origin or ~1000 from multiple tissues?
Thanks,

I expanded a bit my question.

There doesn't seem to be much differences between GISTIC 1.0 and 2.0 as it says:

As with GISTIC 1.0, we obtain P-values for each marker by comparing the score at each locus to a background score distribution generated by random permutation of the marker locations in each sample

But on the supplementary material of the GISTIC 1.0 there is a more detailed explanation of the method. See the "Stage 2" section:

Second, we compare these G-scores to the distribution of scores expected if only random aberrations were observed. This distribution can be determined by rescoring the genome after permuting marker locations within each sample; we instead derive a semiexact estimate.

Furthermore in another section ("Stage 2: Aggregation of Data from Different Tumors to Differentiate Between Driver and Passenger Aberrations") it says:

To determine which of the aberrations identified in Stage 1 are likely to represent driver events, we aggregate the data from all tumors used in the analysis to generate summary scores for amplifications, deletions, and LOH. The statistical significance of each score is determined by comparison to the distribution of scores obtained by all permutations of the data (using a semiexact approximation), with correction for multiple hypothesis testing.

However, the relevant section (of the supplementary materials) seems to be "Null Hypothesis Generation: An Analytic Derivation of the Null Distribution", where it describes the semiexact approximation used.
In a supplementary file it describes it as "Generate all permutations of SNP labels within each sample to simulate datasets with random aberrations"
Conclusion:
It seems that the background probabilities are calculated within sample.
I can't say if it matters how the sets are constituted but I would say that the broader the sets are, the better estimation of the structural variations it will perform but the same estimation of the background will be done.
Ultimately you can check the code of the program, or test with 10 samples and replace one of them to see if the results are changing accordingly.

Thanks, yes it seems to be within sample. I'll do a test of the same sample with different data-sets and report here later.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
Find gene at position from gff or gbk file
I have a VCF file with SNPs from a bacterial genome and want to find if the SNPs are located inside genes, is there some CLI-tool where you can pass a VCF file and a gff or gbk file and it returns the name of the genes?

Please show us an example of your inputs. What format are your SNPs in? Are they in a proper VCF file? Are they in some other format one per line? Many per line? Do you want a command line solution? an online graphical tool? Please [edit] your question and clarify.

Via BEDOPS:
$ gff2bed < annotations.gff > annotations.bed
$ vcf2bed < snps.vcf > snps.bed
$ bedmap --echo --echo-map-id-uniq snps.bed annotations.bed > answer.bed

This can be reduced to a one-liner if you're using bash:
$ bedmap --echo --echo-map-id-uniq <(vcf2bed < snps.vcf) <(gff2bed < annotations.gff) > answer.bed


text:
Standardizing Variant IDs best practise
I'm currently pondering on the best approach to handle and standardize variant IDs within our department to mitigate the challenges associated with rsids and the potential ambiguity of just using chromosomal positions.
The GWAS glue (https://mrcieu.github.io/gwasglue2/articles/Strategy.html) packagae, state oa method where a variant ID is crafted by integrating the chromosomal position (chr), effect allele (ea), and non-effect allele (nea). It assumes the format: chr:pos_ea_nea. For alleles that exceed 10 characters, like certain indels, they apply the murmur32 hashing algorithm from the R digest package.
Given our objectives:
We likely prefer not to have excessively lengthy IDs.
Yet, using hashed IDs for certain long variants might be perplexing for some who aren't informed about the approach.
Hashing every ID could yield consistency, but perhaps might not be as intuitive for regular use.
I'm keen on collecting insights:
Thoughts on the described method from Bristol?
What might be the maximum length for an indel we could encounter?
Is there a broader effort underway to address this concern?
Would adopting the Bristol approach set a potential standard, or should we explore or develop alternative methods? For example should the genome build be built into the id? or species?
I appreciate the community's feedback as we strive for clarity and efficiency in our variant naming conventions.
For clarity:
I want to store data on variants. For example GWAS summary statistics.
I want to be able to match variants from differnt studies so I want to have an identifier for the variants.
I could use rs numbers
I could use chr:pos
I could use chr:pos:ref:alt
I could use chr:pos:effect:noneffect
I could use chr:pos:population:effect:noneffect
or studyid:chr:pos:populatoin:effect:noneffect
or etc
There are pros and cons for each of these but I would like to here what people on here think these are.
I am also intrested in is there work out there to do this already. e.g. to update the idea of rsids to consider multiple alleles and effects in populations?

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

I have tried to clarify

I would not hash. Variant IDs should be static across builds, similarly to how rsIDs are handled between b37 and b38. I do not think hashing is worth it because there's not that much downside for a long ID

You need this to work on multiple species, right? And across different reference genome builds? Are you only interested 8n species with a sequenced genome or do you also need to handle variants defined with respect to a specific transcript or protein only, in the absence of a reference genome?

Also, are you open to paid tools? The company I work for has an API that will, among other things, provide unique ids for arbitrary (human) variants. If you are open to non-free tools check out https://landing.varsome.com/varsome-api.

Thanks ill check that out.

Thoughts on the described method from Bristol?

I don't immediately see the benefit to this over rsIDs.

What might be the maximum length for an indel we could encounter?

Potentially very long. However, why does it matter? This will be rare.

Is there a broader effort underway to address this concern?

dbSNP comes to mind.

Would adopting the Bristol approach set a potential standard, or
should we explore or develop alternative methods? For example should
the genome build be built into the id? or species?

It depends on the type of analysis you're planning on doing. Are you doing cross-species analyses? And if so, what types? Do you expect to be comparing genomes with different genome builds?

For example GWAS summary statistics. I want to be able to match
variants from different studies so I want to have an identifier for
the variants.

Both chr:pos as well as rsID work for this, but just make sure that the genome build matches if you go with chr:pos.

I could use rs numbers I could use chr:pos I could use chr:pos:ref:alt
I could use chr:pos:effect:noneffect I could use
chr:pos:population:effect:noneffect or
studyid:chr:pos:populatoin:effect:noneffect or etc

If you're sure your input data has rs numbers and always will, that might be easiest. However, you may find that sometimes it doesn't have rs numbers (e.g. some Affymetrix data), and it's surpassingly annoying to add them. I'd choose chr:pos, personally, and then consider the effect/non-effect/ref/alt/ on a per-analysis basis, because I would want to still identify the same position if there was a different ref/alt for example.
Edit:
Based on the comment (OP wants to match effect direction), I would suggest chr:pos:effect_allele, but of course this means there is an effect allele (i.e. the position has been studied).

The vast majority of variants (remember that variants are effectively  infinite since you can have any length of indel at any position) won't have rsIDs. Novel variants won't, by definition. So if the OP needs to work with arbitrary variants, rsIDs won't help. If not, rsIDs are indeed the best approach.

@terdon Yes. I think chr:pos would be ideal/most simple, and OP would just have to be careful regarding builds

Thanks for the response and your thoughts.
I guess the thing is in the context of merging GWAS results for coloc and mr etc, chr:pos does not seem sufficeint. I can make sure everything in on the same build, but how do I know if an effect of a snp on say bmi in a japanese study can be matched to a  study which finds an effect at the same position in a european study. I need to take into account ref/alt to effect non effect , direction..

@user27815 I see. I would try chr:pos:effect_allele in that case. Would that work?

gnomAD uses chr-pos-ref-alt, not as an official ID but as a unique way to access web pages detailing each variant. That might be worth looking into. Example: https://gnomad.broadinstitute.org/variant/1-55051215-G-GA?dataset=gnomad_r3

chr-pos-ref-alt is a good way to match variants, but may become problematic (without further tinkering) if the same variant has opposite effect sizes in different studies

yes, which is why I think some times chr:pos:effect:noneffect is often used, then that clears that up. But you still have the issue of population if the effect is differnt in differnt pops..

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Correct for gene length or read counts in GO enrichment analysis
It is a well reported fact that GO analysis of RNAseq results is affected by a number of biases, including length bias and expression level bias.
The bioconductor package goseq allows you to correct for these biases.
By default it corrects for length bias, but you can also get it to do read count bias. Using read counts to do the correction is attractive because in theory it should account for both sources of bias ($read counts\approx expression \times length$).
I'm doing an enrichment analysis were I have tried both options (length and read counts) and get very different answers. If I run a binomial regression on expression and length vs probability of being differential, I can see that both are independently important.

> model <- glm(sig ~  expression + log(length), data=retained_genes,  family=binomial(link="logit"))
> print(anova(model, test="Chisq"))

Analysis of Deviance Table

Model: binomial, link: logit

Response: sig

Terms added sequentially (first to last)

                       Df Deviance Resid. Df Resid. Dev  Pr(>Chi)
NULL                                    6676     4507.1
expression              1  114.998      6675     4392.1 < 2.2e-16 ***
log(length)             1  102.553      6674     4289.5 < 2.2e-16 ***
expression:log(length)  1   34.094      6673     4255.4 5.252e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

So I'm know unsure what to do, should I use the analysis corrected for length or read count. Or perhaps take only terms significant in both? Or only in one?

Did you test for the interaction of them (expression*log(length))?

I had not. Now I have. See updated question.

My understanding of the subject is that the bias of the gene length (and other bias) should be taken care when analyzing the expression and before enrichment analysis. The enrichment analysis should be done once the corrections has been performed. Because as the abstract of the GOseq paper states:

GO analysis is widely used ... but standard methods give biased results on RNA-seq data due to over-detection of differential expression for long and highly expressed transcripts.

So first, take care of the differential expression bias by length and then use the GO to reduce complexity. How you take care of the bias in the RNA-seq data is another question. But the cqn package of Bioconductor can correct the expression by gene length and GC content. However, this correction might hurt the differential tool used (see this discussion in Bioconductor), so it might be better in some cases to use GOSeq.
Now, to the question itself:

So I'm know unsure what to do, should I use the analysis corrected for length or read count. Or perhaps take only terms significant in both? Or only in one?

Use whatever correction method that yields better differentially expressed genes (DEG. If you find that the correction for length improves the accuracy of the predictions of DEG better than correcting by length and GC, then use that one.
Another option to obtain accurate GO terms, then you could use other testing procedures which don't rely uniquely in the Fisher test, such as those that take into account the structure of the GO graph. TopGO use this approach (note that it is a bit difficult to work with this package),this will reduce the role of the gene length bias (and probably other bias) in the resulting significant GO.

These bias' work through power. You cannot magically generate power out of thin air, so there are two ways you could remove the bias pre or at the DE stage: 1) reducing variance between samples 2) down-weighting more highly powered genes. As Mike Love notes in that link, modeling in the row-wise read counts is not compatible with DESeq2s methods (it were possible to adjust for extra power of highly expressed genes at the DE stage, presumably DESeq would already do it). Also who does one measure the accuracy of DEG prediction with no ground truth?

Taking into account the length of the genes won't reduce the variance between samples, so your only option is to down-weight highly powered genes. Maybe there are other methods to take it into account (or you could develop them). To find the true DEG for a given comparison: You can explore the FDR values of the genes in the comparison. The FDR values are assumed to follow an uniform distribution, because each gene has the same probability to be DEG. By exploring how uniform is that distribution you can find if the model captures the biases or not. You could also use the SeqGSA method

But I'm not interested in finding a better set of DEGs, rather, in better assessing whether their division in categories is what you would expect by chance. GOSeq allows you to do this while correcting for bias, but only one bias at once.

The problem in this approach is that the DEG you test, are affected by this length and GC bias. If you want to to assess their division in categories use topGO, because it is the only method which does an ORA analysis taking into account the structure of the GO graph.  GOSeq is biased by not taking into account the relationships between GO terms.

I am not sure if it makes sense to use read counts as bias instead of gene length (and I certainly wouldn't expect the same results).
Do you use total read counts of all your samples (library size)?
The correction for gene length is a pure technical one, the longer a gene the more reads will align (and higher read count genes are easier significant, since they are way above the noise threshold). If you use read counts, you also have a biological factor (expression) in there, which (I think) is the stuff you test with statistics (e.g., with edgeR) and thus not the bias you want to correct for.

Imagine you have two genes: Gene A has 10,000 counts in Condition A and 5,000 in condition B. Gene B has 1,000 counts in Condition A and 500 in condition B. The lfcs are the same, but Gene A is more likely to be call differential than gene B. Now consider two GO categories. Category X is made up of 10,000/5,000 genes and Category Y is made up of 1,000/500 genes. Cat X is more likely to be enriched than Cat Y, despite the real LFC in the genes being the same.

I don't know if the second example is true (that category X is more likely enriched). Maybe you should ask real statisticians on bioconductor?


text:
How to correlate two zero inflated bedgraph-like signals?
This question pertains to iCLIP, but it could just as easily be ChIP-seq or ATAC-seq or mutation frequencies.
I have iCLIP read counts across the transcriptome and I wish to know if the signals are correlated - that is, where one of them is high, is the other likely to be high.
Often when dealing with such data (e.g. iCLIP data) we know that the data is generally sparse - that is at most positions both signals are zero and this is correct, and also zero-inflated - that is many bases that "should" have a signal are missing that data. So just calculating the Spearman's correlation is likely to give an artificially low value.
What might be a way to asses the association? I should add that the aim is the assess association of binding patterns within genes, rather than (or as well as) between genes.
Things I have thought of:

Apply some sort of smoothing to the data (eg a rolling mean). Remove any bases with 0 in both samples. Compute the spearmans.
Calculate the average pairwise distance between every read in sample one and every read in sample two. Compare this to data where the reads have been randomised within genes.

In the first case removing all bases with 0 in both samples seems wrong. But if 99.99% of all bases have zero in both samples, then this seems like its necessary for Spearman.
In the second case, the result seems like it would be non-intuitive to interpret. And also calculating this would be massively computationally intensive.

Honestly I'd just use multiBigwigSummary and then plotCorrelation from deepTools for this, but I'm a bit biased. There, the idea would be to consider each gene as a unit (you could instead use bins, but I don't think that would as nicely do what you want), namely by giving the tools a BED or GTF file input. It would then calculate the average signal in each gene/transcript and you could do your spearman's correlation. Features with 0 in all samples could optionally be removed (plotCorrelation --skipZeros).
While you certainly could go the whole 9 yards and use per-base comparisons, that seems a bit overkill and I suspect that it won't really yield appreciably more information (especially when one considers the additional time overhead).

See my clarification. - I want to look at the spatial patterns within genes, not between them.

@IanSudbery Ah, good to know, my answer isn't applicable then.

It depends whether you want to treat the peak intensities as binary (comparing presence/absence of peaks in the sets) or continuous (comparing the relative magnitudes of the peaks).
Binary
For starting out, a simple binary comparison may be appropriate. You can use a peak caller of your choosing to identify peaks in each sample according to your desired criteria. Then you can use a similarity metric such as the Jaccard index to quantify the level of agreement among the peaks in the two samples.
One potential obstacle is that defining the boundaries of your peaks won't be totally straightforward. For example, a peak in one sample might have 2 overlapping peaks in the other sample, one on each end. A rough solution for this is to divide the genome into bins (maybe around 100-1000 bp, depending on your desired resolution). You can treat a peak as present in a bin if more than half of the peak lies in the bin. That way, bins in one sample can be directly compared to the corresponding bins in the other sample. Obviously, this isn't the only way to do this; other appropriate methods exist too.
Continuous
If you want to treat the peak intensities as continuous, you could apply a similar binning method, taking the "score" of a bin to be the average peak intensity at positions within that bin. You could then throw away all bins with no peaks or only low-intensity peaks throughout the genome. Then you could compute the Spearman's correlation for the remaining bins. I'm guessing it will be harder to find a strong correlation for continuous intensities, because of the amount of experimental variability that is inherently present.
If, after following these steps, the Spearman's correlation is still "artificially low" as you suggested, then this is likely a problem with the underlying data, not the overall analysis; maybe your two datasets actually don't agree that well.

I was worried about artificially low values when most observations were 0. Excluding cases where observations are 0 in both signals just didn't "feel" right.

@IanSudbery I see. So you're also interested in the extent to which the "non-peaks" agree with one another. One possibility then is not to throw away the bins where both values are zero, and to report the Spearman correlation. However, since you said the peaks are quite sparse, your zero values will make it hard to see the extent of correlation between the true positives. I agree that it's not "right" to exclude the zero values in Spearman correlation scores, but if you instead use a metric for set intersections (like the Jaccard index) excluding zeroes would be okay.

Looking at the intersection solves the sparsity problem, but not the zero-inflated one: signals from similar looking tracks tend to be close to each other, but not right on top of each other. Honestly i'm starting to think your second suggest is the best there is.

That's true, but if you use the binning approach to define intersections instead of looking for strict overlap, this could help. (i.e. it counts as an intersection if the same bin contains peaks in both tracks, even if the peaks themselves don't directly overlap)

Rather than working on the base level, you could probably work on say gene level counts. Kendall's tau, an ordinal association metric, can then be used as an appropriate correlation measure.
If $X$ and $Y$ are your iCLIP replicates, $i$ represents gene index and $(x_i, y_i)$ represents the number of RBP binding sites in $X$ and $Y$ respectively for the $i^{th}$ gene, Kendall's tau is defined as :
$$
\frac{\text{#(concordant pairs)} - \text{#(discordant pairs)}}{n(n-1)/2}
$$
Where any two pairs $(x_i, y_i)$ and $(x_j, y_j)$ are concordant if:

$x_i > x_j$ AND $y_i > y_j$

OR

$x_i < x_j$ AND $y_i < y_j$

Correspondingly they are discordant if:

$x_i < x_j$ AND $y_i > y_j$

OR

$x_i > x_j$ AND $y_i < y_j$

Not sure why tex rendering is off.

See my clarification: I'm interested in the patterns within genes, rather than between them. Is tau better at dealing with data where many observations are zero?

If your interest is to study the pattern within genes, does your second sample act as a replicate? I am not sure I completely understand what "within genes" would mean here.

No, I'm trying to determine if two proteins show similar within gene binding patterns.

Mathjax support hadn't been activated. It now has been so your answer renders correctly.

It is one of my favorite stories.
Drop a glance to StereoGene software, it for genomic track correlation, it described in a preprint.
You also can run MACS or another peak caller and estimate the correlation of two interval sets using the GenomtriCorr package.


text:
Getting two class of promoters : ATAC seq data
In this paper they classify the promoters into two categories, highlighted in italics below.

Two Types of Promoters
Mammalian promoters can initiate transcription across broad or narrow
genomic regions, and these promoter shapes, broad or narrow (sharp),
correlate with distinct transcriptional regulatory mechanisms.33 We
evaluated cardiac promoters predicted from CAGE clusters for these 2
major types of promoters by calculating the interquartile range of
promoter CAGE clusters by determining the base pair distance between
10% and 90% of a promoter’s total signal. We observed the expected 2
distinct populations, defined as sharp (interquartile range <10 bp)
and broad promoters (interquartile range ≥10 bp; Figure 2A). Broad
promoters were those associated many different cellular functions,
including housekeeping functions. In contrast, genes with sharp
(narrow) promoters were those encoding proteins critical for
tissue-specific functions seen by the presence myofilament, muscle
tissue, and muscle contractions genes (Figure 2B). Thus,
tissue-specific genes important for LV specification and function were
more likely to have sharp promoters.

Now My question is can I do this only based on ATAC seq data. If yes how do i get that information of narrow and broad promoters.

Its highly unlikely that you will be able to categorize promoters on this basis using only ATAC data. This is a consequence of people using the terms "promoters" and "transcription start sites" interchangeably, when they are actaully different things. CAGE data measures transcription start sites which is the specific location where transcription starts. CAGE does not identify promoters, which are the complete region around the start of a gene that contains, at a minimum, all the sequences necessary to assemble the RNA polymerase II pre-initiation complex. Some definitions of promoter are even broader, including all regulatory sequence proximal the the gene start.
ATAC identifies regions of open chromatin, and I would expect the open chromatin regions around the starts of genes to be more like the broadest definition of promoter, and therefore to be significantly larger than the transcription starts sites they contain, whether they be broad or narrow TSSs.
You could of course look at a distribution of the peak sizes over TSSs and see if you found a bimodal distribution. If you did you could correlate the genes with broad/narrow TSS peaks with their known functions so see if it fitted the idea that narrows TSSs

thats a good start for me to start


text:
Strange p-value histogram for differential gene expression analysis
I'm trying to perform pretty standard differential expression analysis using RNA-seq data.
I've used Kallisto to perform RNA quantification and am using Sleuth to perform the differential expression analysis.
I have 2 conditions that I'd like to compare, which is stored in the metadata:
> metadata
    sample cell replicate batch                   path
1:  Base-1 Base         1   1   /path/to/Batch1/Base-1
2:  Base-2 Base         2   2   /path/to/Batch2/Base-2
3:  Base-3 Base         3   2   /path/to/Batch2/Base-3
4:  Test-1 Test         1   2   /path/to/Batch2/Test-1
5:  Test-2 Test         2   2   /path/to/Batch2/Test-2
6:  Test-3 Test         3   3   /path/to/Batch3/Test-3

I'm attempting to test for differentially expressed genes between the Base and Test cells while controlling for batch (the sequencing batch). The R code to do this is adapted from Sleuth's manual:
sleuth_obj = sleuth_prep(
    metadata,
    extra_bootstrap_summary = TRUE,
    num_cores = 1,
    target_mapping = transcripts_genes, #ENSEMBL information mapping transcripts to genes
    aggregation_column = "ens_gene"
)

# smooth raw Kallisto abundances
sleuth_obj = sleuth_fit(sleuth_obj, ~cell + batch, "full")

# fitting error measurements on reduced model
sleuth_obj = sleuth_fit(sleuth_obj, ~batch, "reduced")

# compare models to calculate differentially expressed transcripts
sleuth_obj = sleuth_lrt(sleuth_obj, "reduced", "full")

# extract results
sleuth_table = as.data.table(sleuth_results(
    sleuth_obj,
    'reduced:full',
    'lrt',
    show_all = FALSE
))
sleuth_table_tx = as.data.table(sleuth_results(
    sleuth_obj,
    'reduced:full',
    'lrt',
    show_all = FALSE,
    pval_aggregate = FALSE
))

To check the quality of the null hypothesis, I've checked the p-value histograms of the gene-level and transcript-level analyses, and I get the following:

Neither of these plots show well-behaved p-values, which makes me question the validity of the reduced model. How can I go about adjusting my analysis to get more reliable results?
Edit following from @swbarnes2's answer:
I've used fdrtool to create an empirical null model of Sleuth's likelihood ratio test.
library(fdrtool)
local_fdr_txs = fdrtool(so$tests$lrt$`reduced:full`$test_stat, statistic="normal")

Which gives the following plots:

Plotting the distribution of these p-values produces:

This still has a hump, but looks more anti-conservative than before.
My remaining question is about whether I can use fdrtool to model the null hypothesis as normal, since I'm not certain about the distribution for Sleuth's likelihood ratio test

what do mean by null hypothesis as normal ? Anyway, a good question +1.

I found that using the Wald test instead of the Likelihood Ratio Test produced very different results with respect to the p-values.
sleuth_obj = sleuth_prep(
    metadata,
    extra_bootstrap_summary = TRUE,
    read_bootstrap_tpm = TRUE,
    num_cores = 1,
    target_mapping = transcripts_genes,
    aggregation_column = "ens_gene"
)

# same model of covariates
sleuth_obj = sleuth_fit(sleuth_obj, ~cell + batch, "full")

# using Wald test instead of LRT
sleuth_obj = sleuth_wt(sleuth_obj, "cellLT1", "full")

# extract results
sleuth_table_genes <- as.data.table(sleuth_results(sleuth_obj, 'cellLT1', 'wt', show_all = FALSE, pval_aggregate = TRUE))
sleuth_table_txs <- as.data.table(sleuth_results(sleuth_obj, 'cellLT1', 'wt', show_all = FALSE, pval_aggregate = FALSE))

The p-values are different, but surprisingly the rankings of the genes aren't too different from each other.

Try the protocol here for generating a different (smaller) dispersion estimate
https://www.huber.embl.de/users/klaus/Teaching/DESeq2Predoc2014.html#inspection-and-correction-of-pvalues

I found this previously, but where DESeq2 has its test statistic be normal, Sleuth detects differentially expressed transcripts using a likelihood ratio test. fdrtool can create an empirical null model on normal data, but is the test statistic from Sleuth's test normal?


text:
How to deal with correlated SNPs in GWAS?
I have population SNP data in VCF format, and I found that some SNPs have a great similarity across samples(> 99%). For example:
CHROM_POS   s0  s1  s2  s3  s4  s5  s6  s7  s8  s9
chr1_1  A   G   G   A   G   G   A   A   G   A
chr1_2  A   G   G   A   G   G   A   A   G   T
chr1_3  C   C   C   A   C   A   C   C   C   C

The similarity between chr1_1 and chr1_2 is 0.9, because they have a different SNP in sample s9.
Is there a good way to remove these similar SNPs before I put the data into another pipeline?
Added: These SNPs will be used for GWAS analysis.

Do you mean that some SNPs have a MAF of <1%?

@Emily_Ensembl MAF is about allele frequency. But I want to remove similar SNP sequence across sampls .

The question is about if removing common SNPs a good thing.

@SmallChess, I think similar snps may will provide duplicate information and cost more resource to compute. Maybe similar snp in long distance should be kept?

It wasn't clear from your original question what you meant. Now you've edited it I can see that you're looking at haplotypic blocks with SNPs in LD. No, don't remove them. There may be causal SNPs within these blocks – you won't able to tell which SNP is causal but you can identify that the block is important.

@Emily_Ensembl, I am not familiar with haplotypic block. SNPs in this block is positional adjacent？ If they are adjacent, I would consider the distance between them, beside their similarity.

Variant that are adjacent are more likely to have alleles inherited together, due to less crossing over between them in meiosis. Pairs of variants are referred to as being in Linkage Disequilibrium (or LD) and blocks of them are called haplotypic blocks.

What you are attempting to do is known as LD-pruning.
As @Emily_Ensembl said, it is not customary to do this for standard association tests: it is possible that one of the SNPs you remove is causal, or a better proxy for the causal locus, and would give (slightly) better association signal than the other. Even for SNPs in perfect LD, pruning is unwise because it would complicate interpretation: if SNP chr1_2 had some obvious function, but you removed it and are left only with its proxy chr1_1 (which could be far away), you will have a hard time connecting the dots. Typical GWAS tests are handled pretty quickly on modern computers, and are linear in SNPs, so there is no computational reason to prune either.
On the other hand, many fancier procedures would require LD-pruned data, because they either
a) are very computationally expensive, or
b) assume independent SNPs.
I suggest reading PLINK's LD pruning documentation, as it pretty much covers the standard ways to do this.


text:
snakemake restricting resources for specific rule
I am trying to generate several snakemake rules. However, I would like to include some restrictions due to storage limitations
I first would like to copy the the bamfiles, 1 for each sample, to a staging area. However, while I have over 50 samples, I can only process three samples at a time.
Then I want to compare the bam file to chunked candidate files, 1 for each chromosome. This however, I would like to run on 23 cores per sample so that once each bam file is copied 69 cores are working.
finally, if all chunks for a given sample have been processed I want to clean up the staging area, copy the bamfile for the next sample and process the chunks for that sample.
My current implementation, see below, just continues copying bamfiles and crashes when I run out of strorage;
rule copy_bam:
    input:
        bamfile=os.path.join(config['bam_folder'], "{sample}.gatk.sorted.MD.bam"),
        bamindex=os.path.join(config['bam_folder'], "{sample}.gatk.sorted.MD.bam.bai")
    output:
        staging_bamfile="{staging_area}/{sample}.bam",
        staging_bamindex="{staging_area}/{sample}.bam.bai"
    wildcard_constraints:
        sample="[A-Za-z0-9_]+"
    resources:
        bam_copy=1
    shell:
        """
        cp {input.bamfile} {output.staging_bamfile}
        cp {input.bamindex} {output.staging_bamindex}
        """

rule process_chunks:
    input:
        chunkfile="{staging_area}/chunks/{chunk}_with_header.bcf.gz",
        bamfile=rules.copy_bam.output.staging_bamfile
    output:
        dynamic("{staging_area}/calls/{sample}/{chunk}.bcf")

    threads: 23
    shell:
        """
        process chunks
        """

rule cleanup:
    input:
        bamfile="{staging_area}/{sample}.bam",
        bamindex="{staging_area}/{sample}.bam.bai",
        processed=lambda wildcards: expand("{staging_area}/calls/{sample}/{chunk}.bcf", chunk=get_chunks()[wildcards.sample], staging_area=config['staging_area'], sample=wildcards.sample)
    output:
        touch("{staging_area}/{sample}/cleanup_done.txt")
    shell:
        """
        rm {input.bamfile}
        rm {input.bamindex}
        rm {input.processed}
        """

run with: snakemake --resources bam_copy=3 --cores 70

One thing to note about any approach here that uses Snakemake's resources is that they're defined per running job, versus with the "in use until something happens in a later job" idea here. Jame's Hawley's suggestion of temp() should help, just note that I don't think the resource tracking even matters in that case (the resource is only "in use" while copy_bam is actively running.)

One Snakemake feature you can use is temporary outputs.
If you make your outputs from the copy_bam rule temporary, they'll be cleared up after no other rule requires that file as an input.
Your Snakefile would look something like this:
rule copy_bam:
    input:
        bam = os.path.join(config['bam_folder'], "{sample}.gatk.sorted.MD.bam"),
        idx = os.path.join(config['bam_folder'], "{sample}.gatk.sorted.MD.bam.bai"),
    output:
        bam = temp("{staging_area}/{sample}.bam"),
        idx = temp("{staging_area}/{sample}.bam.bai"),
    wildcard_constraints:
        sample="[A-Za-z0-9_]+"
    resources:
        bam_copy=1
    shell:
        """
        cp {input.bam} {output.bam}
        cp {input.idx} {output.idx}
        """
# omitted for brevity

You can then specify how many jobs run concurrently with the --cores/--jobs options to keep the space usage small.


text:
Download multiple fastq files using fastq-dump
I want to download the following fastq files at the same time in Salmon:
 - SRR10611214
 - SRR10611215
 - SRR10611215
 - SRR10611216
 - SRR10611217

Is there a way to do this using a bash for loop or fastq-dump? or prefetch

https://www.biostars.org/p/347050/ you can refer to this answer. Here you are trying each read separately. it would be better if you can download and split the read.

A sample code is given in the salmon documentation as follows. Source
#!/bin/bash
mkdir data
cd data
for i in `seq 25 40`;
do
  mkdir DRR0161${i};
  cd DRR0161${i};
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/DRR016/DRR0161${i}/DRR0161${i}_1.fastq.gz;
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/DRR016/DRR0161${i}/DRR0161${i}_2.fastq.gz;
  cd ..;
done
cd ..

This could be modified as follows.
#!/bin/bash
mkdir data
cd data
for i in `seq 14 17`;
do
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR106/0${i}/SRR106112${i}/SRR106112${i}_1.fastq.gz;
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR106/0${i}/SRR106112${i}/SRR106112${i}_2.fastq.gz;
done
cd ..

You can save the code as a shell script and run it from the linux terminal. For example bash download_fastq.sh

Make a list.txt file containing a single column of SRA numbers to download.
then:
for i in $(cat list.txt); do echo $i; date; fasterq-dump -S $i; done
It works well to use NCBI's web interface to find SRA samples of interest, download and open findings in Excel, then copy single column containing SRA numbers and paste into list.txt using document editor such as vim.
After downloading including the "R" can be nice:
for i in *_1.fastq; do mv $i ${i%_1.fastq}_R1.fastq; done
for i in *_2.fastq; do mv $i ${i%_2.fastq}_R1.fastq; done
and zip:
pigz *fastq
If needed a conda option for downloading fasterq-dump:
conda install -c bioconda sra-tools

You can use parallel.
parallel -j 3 fastq-dump {} ::: SRR10611214 SRR10611215 SRR10611215 SRR10611216 SRR10611217

The option -j says how many jobs should maximal run parallel. So in this case maximal 3 identifier would be handled at the same time.
How many jobs you can run parallel depends on your machine.
You can also take a look at parallel-fastq-dump.


text:
Designing a lab NGS file database schema
I am the resident Bioinfo Geek in a hospital academic lab that routinely employs NGS as well as CyTOF and other large volume data producing technologies. I am sick of our current "protocol" for metadata collection and association with the final products (miriad excel sheets and a couple poorly designed RedCap DBs).
I want to implement a central structured, controlled datastore that will take care of this. I know that the interface to the technicians how will be inputing the data is crucial to its adoption, but this is not the focus of THIS particular question: Does there exist a schema or schema guidelines for this type of database?
I would rather use a model that has been developed by people who know how to do this well. I know of BioSQL but it seems more geared towards full protein/nucleotide records like those found in uniprot or genbank. That is not what we have here. What I want is something similar to the system touched on in this preprint: http://biorxiv.org/content/early/2017/05/10/136358
Alternatively, can anyone provide links to where I might find relevant guidelines or supply personal advice?

Are you looking to store processed or unprocessed data?  What would an example file format be that you would be trying to capture?

This is mostly for primary data organization: we get 800 BAMs of WES and I want the file location of each BAM associated with metadata like: PROJECT, READ_LENGTH, SAMPLE_NAME, FAMILY_ID, DATA_TYPE, DIAGNOSIS, etc.

Hey Gus, we are also doing the same research and created this question in Biostars, let us know if you found something! https://www.biostars.org/p/350514/

For metadata, I would use a SQL schema something like the following:
CREATE TABLE Project (
    ac TEXT, -- project/Study accession
    PRIMARY KEY (ac)
);
CREATE TABLE Sample ( -- biological sample/biopsy
    ac TEXT,
    PRIMARY KEY (ac)
);
CREATE TABLE AnalysisSample (
    prj_ac TEXT, -- project acccession (Project.ac)
    symbol TEXT, -- a short name unique in the project
    sample_ac TEXT, -- sample accession (Sample.ac)
    PRIMARY KEY (prj_ac, symbol)
);
CREATE TABLE Collection ( -- a BAM file
    ac TEXT, -- collection/alignment file accession
    prj_ac TEXT, -- project accession (Project.ac)
    PRIMARY KEY (ac)
);
CREATE TABLE ReadGroup (
    cl_ac TEXT, -- collection accession (Collection.ac)
    rg_id TEXT, -- @RG-ID
    sample_sym TEXT, -- @RG-SM; matching AnalysisSample.symbol
    PRIMARY KEY (cl_ac, rg_id)
);
CREATE TABLE VariantSet ( -- a VCF file
    ac TEXT, -- VCF file accession
    prj_ac TEXT, -- project accession (Project.ac)
    PRIMARY KEY (ac)
);
CREATE TABLE VariantSample (
    vs_ac TEXT, -- VCF file accession (VariantSet.ac)
    sample_sym TEXT, -- sample symbol in the VCF file; matching AnalysisSample.symbol
    PRIMARY KEY (vs_ac, sample_sym)
);

In the schema, you have Project and biological Sample tables, which are independent of each other at the high level. An AnalysisSample describes a sample used in BAM or VCF and connects Project and biological Sample. Importantly, each AnalysisSample has a symbol unique in a project (see the primary index). This is the symbol on a BAM read group line or on a VCF sample line. A Collection is in effect a BAM/CRAM file. In theory, a BAM file may contain more than one samples (though rare in practice), which is addressed by a separate ReadGroup table. Finally, a VariantSet is a VCF file. VariantSample tells you which samples are included in each VCF file.
This is the skeleton of a full schema. You can add extra fields to appropriate tables (e.g. file path and hg19/hg38/etc to Collection, read length to ReadGroup and family ID to Sample). You also need indices for efficient table joining and perhaps more tables for complex structures (e.g. pedigree).
For the projects I have participated, this schema should work most of time. It is inspired by GA4GH's JSON schema, but my version is in SQL, is simpler and also has a slightly different structure which I think is better.

The Global Alliance for Genomics and Health has been working on the issue of representing sequencing data and metadata for storage and sharing for quite some time, though with mixed results.  They do offer a model and API for storing NGS data in their GitHub repository, but it can be a bit of a pain to get a high-level view.  I am not sure if any better representation of this exists elsewhere.
I can say from personal experience (having built over a dozen genomic databases), there is no ideal data model and storage best practices.  Genomic data comes in many shapes and sizes, and your needs are going to vary from every other organization, so what works for one bioinformatics group won't necessarily work for you.  The best thing to do is design and implement a model that will cover all of the data types in your workflow and downstream analyses you might do with the data and metadata.

I agree that there is no ideal data model that is going to be stable for very long in a quick-moving field like genome informatics. Perhaps a schema-less (NoSQL or some other document-based system, such as MongoDB) database approach would work better? This gives you ultimate flexibility to attach whatever information is relevant to database entries you're adding to your database now, without the need to rebuild the database later if you want to attach more/different information to subsequent database entries.


text:
Is there a way to retrieve ENSEMBL IDs from a search query?
I would like to be able to batch download FASTA files from ENSEMBL. I normally would use the API to download them from the ENSEMBL gene IDs. However, the IDs aren't ones I know, only ones I would get from searching on the ENSEMBL website and collecting the IDs manually.
Is there a way I can use the API to query search terms and collect all results that are returned?
So for example, if I use the ENSEMBL website to search "HLA" I could collect the list of IDs (ENSG00000204252 etc. etc.) manually. HLA works okay using the API to query gene IDs because there's only a few terms (A, B, C plus non-classical for class I) but I was wondering if there was a way to directly access the search query one the homepage programmatically, as I may need to do it for some messier examples.

Can you give an example of the sort of search you're running? There's lots of things in the API that might do what you want but don't think there's a generic "search" call.

Yes, sorry I didn't think to do that. I want to search a term associated with a gene family, say "HLA" and retrieve ENSEMBL IDs (or the sequences associated with, but I already have some code to do that) of the results that come up from the search.

What kind of terms are you looking for: Diseases terms, protein family names (KRT7,, KRT19, KRT11), ...? Why do you need to do so (maybe we can work around this problem) ?

Please [edit] your question and give us a few example IDs and the output you would like to get from them.

Yes, this is possible but you really need to give us specific examples. You say you query for HLA, OK, but what do you want to getr from it? The FASTA sequence of what? The gene? Primary transcript? All transcripts? Peptide? All peptides? What species? All species? This will almost certainly be possible either directly from the API or, at worst, by processing the HTML but we can't help unless you tell us what you need.

Sorry I've been trying to keep it relatively agnostic. I just need the ENSEMBL gene IDs for the human alleles. The problem is not what happens downstream of that, I just wanted a way of getting IDs of genes and their alleles in one go rather than searching each indivdual gene and getting the gene ID.

Yes, but it's hard to understand what you need without an example. You say you want to download FASTA, but also that you want to retrieve IDs. If you give us an example, we can show you how to do it in a way that is applicable to anything but we need an example to understand exactly what you need.

I would download the complete geneset FASTA file from ensembl and then parse out the genes you want with a script.
Alternatively, you could do similarly with API. Get all the gene names, select the ones you want, and then pull down each sequence.


text:
Why is this makeblastdb command not working?
I am trying to make local databases for the ncbi-blast+ package (version 2.60). I am doing so for 4 T-cell receptor genes. 3 of the 4 (TRAV, TRAJ, TRBV) have worked fine, but I am having problems with TRBJ. I say this because they should be the same format and if one works the rest works.
They take the format below:
>K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48 nt|3| | | |15 AA|15+0=15| | |
NTEAFFGQGTRLTVV

I used the following command (like I say, it was fine for the rest, oddly enough):
makeblastdb -in TRBJ_AA.fa -parse_seqids -dbtype prot

and I get the following error message:
BLAST options error: TRBJ_AA.fa does not match input format type, default input type is FASTA

One issue I checked was the length of the sequences being too short, but I believe 11 is the minimum length required (for inputs, don't know if that is the same for databases)- but would someone be able to confirm, that strikes me as the most likely candidate?
I'm happy to use some other metric, but I was hoping on keeping BLAST for consistency's sake for the time being.
Thank you.
edit: I know it's to late, but for reference I'm on Ubuntu 16.04.2 LTS x86_64 for all steps.

What operating system are you running this on and what operating system was the input file created on? Also, please [edit] your question and show us the actual file you are using.

This looks like a bug in makeblastdb. Removing the final | | | from your sequence's description makes it work:
>K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48 nt|3| | | |15 AA|15+0=15

So does removing everything after the 1st space:
>K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48

There is nothing in the FASTA definition that would explain this and, even more weirdly, the presence of the | | | alone isn't enough to trigger the issue. I played around with this a little and it looks like it doesn't like the combination of + or -along with empty | |, but there's no clear pattern I can see:
These fail:
>| | |15-| | |
>| | |15+| | |
>| | |15-1| | |
>| | |15+1| | |

These work:
>| |15+| |
>| | |151| | |

So yeah, sounds like a bug. I suggest you report it to the NCBI devs. In the meantime, a simple fix would be to just keep the first space-separated field of your fasta header. Pass your fasta file through this:
awk '{print $1}' fasta.fa > new.fa

That will print the 1st field of each line and, assuming you don't have spaces in your sequence, it should effectively sort this out for you.

I'm not convinced that this is a bug per se -- FASTA headers can use any combination of these specific identifiers. It's parsing it somehow (assuming it's figuring out which of the possible identifiers to use based the parsing the + and - as strings), but it must be getting confused by the addition of three empty columns.

@OrdiNeu what it is after > is left for the user to decide if NCBI decides to add these identifiers or more is their own decision but it should be able to work with other identifiers

@OrdiNeu the FASTA format itself allows anything you care to put there. The identifiers you linked to are just what NCBI tends to use in its FASTA files, but they have nothing to do with the FASTA format. If a program expects a specific format, then the program is buggy. The fact that there's no clear pattern and that the issue occurs both with and without the -parse_seqids flag is further indication of a bug.

Oddly I found on my system that whether or not it 'works' also depends on the sequence lines. For example, take header >| | |151| | | and sequence ATGAAGAG it fails, with ATGAAGAG\nATG it works (in both cases running makeblastdb  -dbtype nucl -in my_fasta.fa)

@Chris_Rands curiouser and curiouser. But that makes it even more likely that this is a bug.

@terdon Yes, can you replicate what I described?

@terdon That's ahem... interesting! The 3 empty columns can be found in the other IMGT sequence headers, which were parsed successfuly- whereas the numerics with +/- are only found in the problem file which does place the blame firmly at their door. It also still doesn't work when leaving out the -parseseq_id flag. Will report it. Thanks for maintaining my sanity.

@Chris_Rands yes and no. The one you say works, kinda works on my system but complains about BLAST Database creation error: Defline lacks a proper ID around line 1. Changing the defline to >a| | |151| | |, however, reproduces what you describe perfectly: sequence ATGAAGAG fails but sequence ATGAAGAG\nATG works. Weird.

@TW93 you're welcome. I know how you feel. I'd have wondered about my own sanity too had I encountered this in the wild.


text:
What is the aim of insertion codes in the pdb file format?
In PDB file we always see insertion codes (ic) in a few proteins. Could anybody please tell me why it has been added in structure and what is the aim of that ???

Hi @user7355 there isn't enough information to answer this question. An example protein would useful together with information about the protein e.g. resolution (Amstrongs), method of structural determination

@MichaelG. the insertion code has a well defined meaning, it's part of the sequence ID.

When I wondered it myself years ago and I found a good explanation in a mailing list archive. I think it was this one, by Frances C. Bernstein:

The problem with numbering started when people wanted to compare the
  'same' proteins from different species.  They found that there were
  the following possibilities that gave rise to differences:

More or fewer residues at either end.
Extra residues at various places within the chain.
Fewer residues at various places within the chain.
Different amino acids at the same place.

Now imagine that residue PHE 195 is very important for the activity of
  the protein in species A.  But in species B it is residue PHE 197 and
  in species C it is PHE 212, because species B and C are not the same
  length as species A.
Because people felt it was important to preserve the amino acid numbering for 'important' residues and to be able to readily discuss
  and compare the structures from different species, various people
  decided to try to number the proteins from species B, C, etc. to match
  the numbering used for species A.  In doing this, one must have gaps
  (missing numbers) where a sequence is shorter.  But what should one do
  when a sequence is longer?  This is the case where it is necessary to
  insert extra numbers and this is done by using insertion codes.
Thus the insertion code is an integral part of the residue number and it is improper to ignore that field when using a PDB entry.  You
  must also allow for 'missing' numbers when using a PDB entry.

In the PDB format the insertion code is in column 27, following the sequence number (columns 23-26):
                      nnnnI
ATOM  11918  CZ  PHE D 100      -6.852  76.356 -23.289  1.00107.94           C
ATOM  11919  N   ARG D 100A     -9.676  74.726 -19.958  1.00105.71           N
...
ATOM  11970  CE  MET D 100H     -8.264  83.348 -19.494  1.00107.93           C
ATOM  11971  N   ASP D 101     -11.329  81.237 -14.804  1.00107.41           N

Together, columns 23-27 make a sequence ID, which is a mostly-numeric ID of the
residue position.
Note: It's not a full ID of the residue. When you have a point
mutation (a.k.a. microheterogeneity) you have 2+ residues with partial occupancy at the same position. To fully identify a residue both the sequence ID and residue name are needed.
As of 2020, the PDB (organization) uses PDBx/mmCIF as the primary format. The PDB format is only a legacy one. The presentation of the insertion codes (and many other things) in the mmCIF format are rather confusing, so it's worth explaining.
The mmCIF format introduced double identifiers for chains, residues and atoms: author's ID and "label" ID. Both look the same and are easy to confuse. For example: two chains named B and A by the author are renamed as, respectively, A and B by the wwPDB; both names are used in parallel: one chain is B(A), the other A(B). Poor design choice at every step. The same with the sequence ID. Originally, we had two tags:

_atom_site.auth_seq_id -- author's sequence ID, possibly with an insertion code (100A),
_atom_site.label_seq_id -- numeric-only sequence ID, numbered with respect to the sequence of this molecule (SEQRES).

Then, as I heard, numeric insertion codes were spotted in the wild. (In a database big enough any mistake that is not automatically checked for will happen). This made a problem with translation from the PDB format to mmCIF. The sequence number 100 with IC 1 would be interpreted as sequence number 1001. The numeric ICs must have been remediated in the meantime, because they absent as of 2020. But at that time (10-20 years ago?), instead of correcting the entries, the decision was made to move the insertion code to a separate column (pdbx_PDB_ins_code).
The documentation was never updated, so the auth_seq_id is still documented to store the sequence ID, not just the sequence number.
Here is an example from 3B9F:
                             !!  I                                      nn
ATOM   770  C CD2 . LEU B 2  46  ? -23.416 21.041  -27.401 1.00 19.94 ? 60   LEU H CD2 1
ATOM   771  N N   . TYR B 2  47  A -21.191 16.717  -28.119 1.00 18.73 ? 60   TYR H N   1

The sequence IDs here are 60 and 60A. Since the meaning of the insertion code is not widely known, it landed next to label_seq_id which has nothing to do with it. So again: the insertion code is just a part of the author's sequence ID (originally, in auth_seq_id), and the letter A above (10th field) is simply an extension to the number 60 (17th field).


text:
Profile of conservation between bacterial sequences and human protein
My goal is to make a conservation plot between bacterial sequences and a human protein.
So far, I have a FASTA file of the protein, and a FASTA file with the sequences of the proteins from the BLAST results.
In my initial attempts, I have tried to do this using msa software:
I have downloaded clustalx and run a profile alignment using the single protein sequence as Profile 1 and the sequences from the BLAST results as Profile 2, and selected the option "Align Sequences to Profile 1".
I'm having difficulty interpreting and dealing with the results. Since the majority of the sequences aren't aligned, the results are a total mess, with a huge number of dashes added.
My goal is to visualize the number of BLAST hits per amino acid in my protein of interest.
ie:
I would like to have a graph with my protein on the x-axis and a plot like the regions of high conservation for the multiple alignments, except with the spikes corresponding to a high number of BLAST hits. This would allow me to identify regions of my protein that have higher sequence similarity to bacteria than others.
Is there a better way to achieve this or to salvage the results from the profile alignment?
Thanks.

I curious why you want to use a profile alignment. Wouldn't it be easier to just add your sequence to the sequences from the blast results and run a standard multiple sequence alignment?

Also, could you clarify what you're trying to achieve? You say you want the number of blast hits per amino acid, but the figure you show is plotting amino acid conservation. Do you care if the amino acids are the same at a given position or just if the high-scoring segment pairs overlap the position?

We could help much more if you gave us the alignment in text form or even a few of your sequences so we can align them here. The image is informative but we can't really work with it. That said, please do explain what you're trying to do here. I am having trouble understanding what the number of blast hits per amino acid would convey. Why would you need that value? And are you sure you want all of those sequences? If the majority are not aligned, are you sure they're relevant? What is the actual question you are trying to answer using this alignment?

A profile alignment generates a new alignment from two sets of aligned sequences (the "families"), where one particular case is between one family and one sequence. Therefore to run clustalx, you need to first align the proteins from BLAST results, before using them as your Profile 2. The other questions remain, about number of blast hits or final goal...maybe there is some confusion with the position-specific score matrix from PSI-BLAST, which also generates a score?

Is there a way to visualize the position specific scoring matrix generated from a PSI-BLAST with a line plot that would resemble the graph from the msa in my original post? My final goal is to produce a plot visualizing regions of high bacterial sequence similarity to my human protein of interest.

@bluescholar1212 have you looked at Jalview? It is designed as an MSA viewer, plus has many ways for viewing conservation which you can extend yourself.

@ithinkiam I have taken a look at Jalview and it looks like it produces some great visualizations. Would a BLAST result be a valid Jalview input?

AFAIK, no it wouldn't be a valid input, but you could use the blast ids to make the alignment and provide the multiple alignment to Jalview. Also please consider [edit]ing the question, you want to have a profile of conservation between bacterial sequences and the human protein, isn't it?

@Llopis thanks, I went ahead and edited my original question, and am going to add a better description as well shortly.
Regarding how best to achieve this, am I better off using a different approach, or is that a reasonable workflow to blast, align the results, and visualize with Jalview? I think potentially a line plot showing consensus percentages of the results against the human protein would be an appropriate metric.

I'm not sure what you mean by it being a mess? That looks like a pretty good alignment to me, and it most certainly has aligned all the sequences. You are nearly always going to have gaps (see my MSA below, which are all paralogs from the same genome so they couldn't really be more closely related, and yet, there are gaps). That comes with its own set of problems mind you, as you need to decide how you're going to deal with them.
You also don't necessarily have to do a profile alignment. Those sequences don't look massively divergent to me, so straight forward sequence alignment would probably give you more or less the same result.

My goal is to visualize the number of BLAST hits per amino acid in my protein of interest.

This doesn't really make sense. You should really be searching for the number of BLAST hits with a particular domain, if you want to infer homology.
If you want a per-position visualisation of the conservation, you could look at the shannon entropy for each column and plot that. I wrote a script to do just that a little while ago: https://github.com/jrjhealey/bioinfo-tools/blob/master/Shannon.py
Just beware it's not super well tested yet. Feed an MSA in with as many sequences as you want to analyse, but you'll have to have identified the sequences and done the alignment first.
For example, given this MSA:
    16    149
PAU_02775  MSTTPEQIAV EYPIPTYRFV VSLGDEQIPF NSVSGLDISH DVIEYKDGTG
PLT_01696  MSTTPEQIAV EYPIPTYRFV VSIGDEQIPF NSVSGLDISH DVIEYKDGTG
PAK_02606  MSTTPEQIAV EYPIPTYRFV VSIGDEQVPF NSVSGLDISH DVIEYKDGTG
PLT_01736  MSTTPEQIAV EYPIPTYRFV VSIGDEKVPF NSVSGLDISH DVIEYKDGTG
PAK_01896  MTTTT----V DYPIPAYRFV VSVGDEQIPF NNVSGLDITY DVIEYKDGTG
PAU_02074  MATTT----V DYPIPAYRFV VSVGDEQIPF NSVSGLDITY DVIEYKDGTG
PLT_02424  MSVTTEQIAV DYPIPTYRFV VSVGDEQIPF NNVSGLDITY DVIEYKDGTG
PLT_01716  MTITPEQIAV DYPIPAYRFV VSVGDEKIPF NNVSGLDVHY DVIEYKDGTG
PLT_01758  MAITPEQIAV EYPIPTYRFV VSVGDEQIPF NNVSGLDVHY DVIEYKDGIG
PAK_03203  MSTSTSQIAV EYPIPVYRFI VSIGDDQIPF NSVSGLDINY DTIEYRDGVG
PAU_03392  MSTSTSQIAV EYPIPVYRFI VSVGDEKIPF NSVSGLDISY DTIEYRDGVG
PAK_02014  MSITQEQIAA EYPIPSYRFM VSIGDVQVPF NSVSGLDRKY EVIEYKDGIG
PAU_02206  MSITQEQIAA EYPIPSYRFM VSIGDVQVPF NSVSGLDRKY EVIEYKDGIG
PAK_01787  MSTTADQIAV QYPIPTYRFV VTIGDEQMCF QSVSGLDISY DTIEYRDGVG
PAU_01961  MSTTADQIAV QYPIPTYRFV VTIGDEQMCF QSVSGLDISY DTIEYRDGVG
PLT_02568  MSTTVDQIAV QYPIPTYRFV VTVGDEQMSF QSVSGLDISY DTIEYRDGIG

           NYYKMPGQRQ AINISLRKGV FSGDTKLFDW INSIQLNQVE KKDISISLTN
           NYYKMPGQRQ AINISLRKGV FSGDTKLFDW INSIQLNQVE KKDISISLTN
           NYYKMPGQRQ AINISLRKGV FSGDTKLFDW INSIQLNQVE KKDISISLTN
           NYYKMPGQRQ AINITLRKGV FSGDTKLFDW LNSIQLNQVE KKDISISLTN
           NYYKMPGQRQ LINITLRKGV FPGDTKLFDW LNSIQLNQVE KKDVSISLTN
           NYYKMPGQRQ LINITLRKGV FPGDTKLFDW LNSIQLNQVE KKDVSISLTN
           NHYKMPGQRQ LINITLRKGV FPGDTKLFDW LNSIQLNQVE KKDVSISLTN
           NYYKMPGQRQ SINITLRKGV FPGDTKLFDW INSIQLNQVE KKDIAISLTN
           NYYKMPGQRQ SINITLRKGV FPGDTKLFDW INSIQLNQVE KKDIAISLTN
           NWFKMPGQSQ LVNITLRKGV FPGKTELFDW INSIQLNQVE KKDITISLTN
           NWFKMPGQSQ STNITLRKGV FPGKTELFDW INSIQLNQVE KKDITISLTN
           NYYKMPGQIQ RVDITLRKGI FSGKNDLFNW INSIELNRVE KKDITISLTN
           NYYKMPGQIQ RVDITLRKGI FSGKNDLFNW INSIELNRVE KKDITISLTN
           NWLQMPGQRQ RPTITLKRGI FKGQSKLYDW INSISLNQIE KKDISISLTD
           NWLQMPGQRQ RPTITLKRGI FKGQSKLYDW INSISLNQIE KKDISISLTD
           NWLQMPGQRQ RPSITLKRGI FKGQSKLYDW INSISLNQIE KKDISISLTD

           EAGTEILMTW SVANAFPTSL TSPSFDATSN EVAVQEITLT ADRVTIQAA
           EAGTEILMTW SVANAFPTSL ISPSFDATSN EVAVQEITLT ADRVTIQAA
           EAGTEILMTW SVANAFPTSL TSPSFDATSN EVAVQEITLT ADRVTIQAA
           EAGTEILMTW SVANAFPTSL TAPAFDATSN EVAVQEISLT ADRVTIQAA
           ETGTEILMSW SVANAFPTSL TSPSFDATSN DIAVQEIKLT ADRVTIQAA
           EVGTEILMTW SVANAFPTSL TSPSFDATSN DIAVQEIKLT ADRVTIQAA
           EAGTEILMSW SVANAFPTSL TSPSFDATSN DIAVQEIKLT ADRVMIQAA
           ETGSQILMTW NVANAFPTSF TSPSFDAASN DIAIQEIALV ADRVTIQAP
           EAGTEILMTW NVANAFPTSF TSPSFDATSN EIAVQEIALT ADRVTIQAA
           DAGTELLMTW NVSNAFPTSL TSPSFDATSN DIAVQEITLT ADRVIMQAV
           DAGTELLMTW NVSNAFPTSL TSPSFDATSN DIAVQEITLM ADRVIMQAV
           DTGSEVLMSW VVSNAFPSSL TAPSFDASSN EIAVQEISLV ADRVTIQVP
           DTGSKVLMSW VVSNAFPSSL TAPSFDASSN EIAVQEISLV ADRVTIQVP
           ETGSNLLITW NIANAFPEKL TAPSFDATSN EVAVQEMSLK ADRVTVEFH
           ETGSNLLITW NIANAFPEKL TAPSFDATSN EVAVQEISLK ADRVTVEFH
           ETGSNLLITW NIANAFPEKL TAPSFDATSN EVAVQEISLK ADRVTVEFH

You'd get this plot:

I would like to have a graph with my protein on the x-axis and a plot
  like the regions of high conservation for the multiple alignments,
  except with the spikes corresponding to a high number of BLAST hits.
  This would allow me to identify regions of my protein that have higher
  sequence similarity to bacteria than others.

I'm not sure your logic is quite right here though. Blast won't give you hits depending on a particular position. It's a local aligner, so it'll just return you hits where at least some part of your query matches at least some part of another.
What you could do is take the logic in the script above, and just use a different metric. For example, perhaps you could count the proportion of sequences which have the most common amino acid at a given position within your MSA. That would be fairly crude though.
As you say in the comments,

My final goal is to produce a plot visualizing regions of high bacterial sequence similarity to my human protein of interest.

your original alignment will show you this intrinsically, if only you include the sequences of all the BLAST hits in the first place. Thus your work flow will be:

Blast sequence of interest.
Download all/as many hits as you want (bear in mind the E-value/bitscore and the number of hits you get. It might only be a few dozen, in which case you can use the lot, but if not, just take all the hits below a certain cut-off.)
Align all the sequences.
Look at the column scores for the whole MSA. You can use whatever metric of conservation you like really. Might be as simple as proportion of sequences with the most common residue, or something more complex like Shannon entropy (though as you can see in the graph above Shannon entropy can be kinda noisy) etc.

this is an excellent answer. I think that this is the work flow that I will adopt. While researching this I also found this blog post link with a script that produces a pretty smooth conservation plot using the Savitzky-Golay technique. The script seems to follow a similar workflow to what you proposed. Now I will need to work on which blast parameters and metrics of conservation will be appropriate for this project. Thanks for your help.

If you don't mind abstracting the data slightly, the script can calculate a moving average/running mean which smooths the data a bit.

this may be worth posting a new question to answer, but in the link that I posted, the author mentions you could easily plug in a new function for the remote_blast function depending on usage. I have a list of parameters I would like to use for my blast, but ideally I would like to access the remote database. Do you know how I could change the function in that script to specify more specific parameters that I need?

I'm not sure what you mean exactly? With the python implementations of remote blast you should have most, if not all, the options of a standard blast. What sort of parameters are we talking about? It's easy enough to add additional arguments to the argparse framework if you want to specify blast thresholds etc. It's pretty much just a case of adding in another function or two to do the heavy lifting as far as generating the MSA for you goes, and then that can be passed directly to the other functions of the file.


text:
How to identify the reference genome assembly used to create a .hic file?
I have a bunch of .hic files of different species. But, I don't know which reference genome assembly (e.g. GRCh37 or GRCh38 in the case of human genome) the genomic co-ordinates of the .hic file correspond to. So I wonder if there is any way to figure out the version of the reference genome assembly from a .hic file.
Therefore, my question is: Does .hic file contain metadata such as the version of reference genome assembly? If yes, then how could I access it?
(I could obtain the information I need by digging into the methods sections of each paper corresponding to each .hic file. But it would be much better if I could get that right-away from the .hic file.).

Turns out HiC file does store such metadata.
straw tool from Aiden Lab contains a script (read_hic_header.py) that prints out the metadata from a given HiC file. Usage: read_hic_header.py <hic file or URL> [verbose]
In order to obtain the metadata in the form of a dictionary, one may consider using a function of the same tool that is currently available in my fork.
from straw import straw_module
straw_module.read_metadata(infile,verbose=False)

Args
    infile: str, path to the HiC file
    verbose: bool, print the metadata if True else return the metadata
Returns
    metadata: dict, the metadata is returned if verbose is False.
                Keys of the metadata:
                HiC version,
                Master index,
                Genome ID (str),
                Attribute dictionary (dict),
                Chromosomes (dict),
                chromosomes2length (dict),
                Base pair-delimited resolutions (list),
                Fragment-delimited resolutions (list).


text:
Get a single number representing the contact probability between a pair of genomic loci using HiC data
Hello to the experts in analyzing chromosome structure data in HiC format,
I have a very basic question. I have a HiC file (specifics are mentioned below), using which I would like to obtain contact probability between a pair of genomic loci.
I obtained contact probabilities between the two loci at different resolutions available in the HiC file. However, rather than having a vector of contact probabilities, I would like to have a single number that appropriately represents the contact probability between the loci.
My guess is that there could be two ways to get a single number as the contact probability between two given genomic loci.

Should I take the highest contact probability?
or,
Depending on the lengths of the genomics loci of interest, may be the resolution that is close to that could be selected as a single number that represents the contact probability between the loci. Eg. with this strategy, if genomic loci of interest are ~5kb in length the contact probability at 5kb resolution is selected as a single number. Then, the issue is that the length of the genomic loci of interest could be lower than the lowest length of resolution available in the HiC file.
or,
Contact probabilities at different resolutions could be aggregated. Then I wonder what would be the appropriate way to aggregate contact probabilities, e.g. by multiplication or median or something else.

Just to reiterate, above these are my guess and I am not sure whats the best way.
So basically my question is: how can I get a single number that would represent the contact probability between a pair of genomic loci using HiC data?

Some specifics that might be useful in answering my question.

source of the HiC file: https://ftp.ncbi.nlm.nih.gov/geo/series/GSE63nnn/GSE63525/suppl/GSE63525_K562_combined_30.hic
organism: human (K562 cell line)
Resolutions in HiC data (bp): [2500000, 1000000, 500000, 250000, 100000, 50000, 25000, 10000, 5000]
Type of genomic loci: upstream non-coding regions of the genes
lengths of genomic loci of interest: average ~1000bp (range: 500-5000bp).
software I am using to "extract" data from HiC file: https://github.com/aidenlab/straw

The fundamental unit of a Hi-C contact is the nearest restriction enzyme cut site.
This is where each mate in a read pair from a Hi-C library comes from and is mapped to when converting sequencing data into a contact matrix.
Because the distance to the nearest cut site is variable at different genomic positions, and depth of coverage of any single locus is typically low, contacts are aggregated into bins and "normalized" to account for differences in cut site densities within those bins.
This normalization transforms the raw contact matrix into something resembling a stochastic transition matrix, where you can think of the elements of that matrix as "what is the probability that at any moment in time bin $i$ is in contact with bin $j$".
So the probability you're looking for is the $(i, j)$th element of the normalized contact matrix, where bins $i$ and $j$ contain your regions of interest.
If your regions straddle the bin boundaries, you may want to average those bins.
You get to pick the resolution of the contact matrix, but if you select a resolution that's too small, you run the risk of making your matrix too sparse and having no contact between those loci.
If you pick one that's too large, your matrix elements aren't really informative for your question at hand.
Given your 500-5000 bp range and median of 1 kbp, it may make sense to pick a resolution between 1 kbp and 5 kbp and use that.


text:
Python API for working with ENSEMBL genomes
Dear bioinformaticians,
I could find two python APIs to mine ENSEMBL genome databases. There is a pyensembl which is more popular and then there is genomepy which is less popular. Following is my experience with both.
1. pyensembl:
pyensembl is a fast and powerful API to mine ENSEMBL genome databases.
I have been using it for a while and in my experience it works very well. The name of the package may make you believe that it is developed/maintained by ENSEMBL. That is not true. Still, being an open source package with a realistic potential to be the default API for mining ENSEMBL genome databases.
However, there are some issues related to it's development that users may find disconcerting.
For example, documentation. I could only find a comprehensive api documentation for release 0.8.10. Current release is 1.1. Not even release notes are available. Therefore users have to rely on few lines of documentation that is available in the README.md.
Also, as far as I know, it is not peer reviewed even though it is an excellent analysis resource.
2. genomepy
genomepy loads all of the genome files into the RAM, making it immediately unusable, especially if you are working with larger genome such as human genome.
So my question is: Is there any other better python API for working with ENSEMBL genomes?

"That makes you slightly worried about it's reliability."  You have the source code, you can see what it does and you can copy and modify as you wish. Why do you think it is unreliable? You can improve the documentation and suggest changes, but it looks like it is documented well with docstrings although the documentation has not been rebuild in a pdf format in a couple of years.

@Llopis I agree with both points. (1) I just wonder why it is not published. It is is a very useful resource. I will remove the reliabilitypart from my question. Thanks.
(2) Conventionally, for contributing (making pull requests), a potential contributor has to understand the code and to understand the code, the he/she has to go through the  documentation (preferably the latest one). So that's why I feel documentation is an important issue. May be I will build the documentation myself and contribute.

If we wrote papers on every tool we create we'd never have time to do anything productive. That and the limited benefit of actually submitting a paper on a tool/package is the reason that most never get published in that way.

The correct API for Ensembl is the Ensembl REST API which is updated and maintained by Ensembl, and language agnostic.

Hi @Emily_Ensembl, I wonder if there is a easy-to-learn kind of tutorial on Ensembl REST API. I want to use it to carry out general tasks such as (1) converting gene ids to symbol and other way around and (2) fetching a genome sequences of a particular species (assembly and ensembl release) within specific coordinates.

We've got an online course here.

Now there's another option, the ensembl_rest module, a thin wrapper around the Ensembl REST API to simplify its usage and make it more pythonic. You can find the documentation here.
To clarify things, I'm the creator and maintainer, but still think it's a legitimate alternative.

This looks pretty neat and actually a lot more cleaner than the default python option in the EnsEMBL page. Great work!


text:
How to link GDC ids to CCLE cell line names?
Hello bioinformaticians,
I've recently downloaded few RNAseq data files from Genomic Data Commons Data (GDC) portal. These files belong to Broad institutes CCLE project.
Now the problem is that GDC portal do not provide CCLE cell line names:

So my question is how could I link the GDC ids (prefixes of File Names) to CCLE cell line names?
I would imagine because it is an important step in the analysis, there should be some way to do that.

Where have you searched? Is there any information page about how to identify the cell lines or the files?

It's an issue with GDC portal. So I looked around on their website. There is not even an email address to query about this issue.

The filenames of those bams contain the cell line names. They are constructed:
SAMPLEID.CELL_LINE_NAME.REPLICATE.bam
Thus the first 5 cell lines in your screen capture are SK-N-F1, TE-9, CJM, HCI-H1915 and CAL-33.
All the information that is need is also present in the experimental metadata xml file attached to each sample record. This is accessible by downloading it from the record on the website, or by querying the GDC metadata API (https://gdc.cancer.gov/developers/gdc-application-programming-interface-api).

Thanks. There is indeed some information in the filenames. However, the CCLE cell line names  are generally in the format {CELL_LINE_NAME}_{TISSUE_NAME}. For example, CJM_SKIN and CAL33_UPPER_AERODIGESTIVE_TRACT.
I'm not sure where can I get the tissue names now.

If you look at the "project" column in the above table, it gives the cancer type code from which the cell lines were extracted.  Now how to map those abbreviations is less obious, but if you click a primary site filter in the left hand column, then you can work it out.

This information is also present in the experimental metadata xml file attached to each bam file. This probably means its also accessible from the GDC metadata API.

Ok. I see. There are XML files along with the BAM files.

GDC API (https://gdc.cancer.gov/developers/gdc-application-programming-interface-api) should be the best way to dig into those. Thanks.

Could you add the xml file related comment in the answer so I can accept it as the correct answer?

I would recommend using UCSC Xena instead of GDC for CCLE data. They don't have everything. For example, BAM files are not available. However, the data that is available is organized in standard text files with consistent identifiers.


text:
Convert Entrez gene ids and HUGO symbols, whose genome assembly is unknown, to GRCh38 Ensembl gene ids
I want to use a file available on a database called DepMap (https://depmap.org/portal/download/). This file called CRISPR_gene_effect.csv contains gene ids in the format HUGO (Entrez) eg. TP53 (7157), without having any complementary information regarding which genome assembly (GRCh37/38) the Entrez ids (and also HUGO symbols) belong to.
I want to convert these gene ids in the HUGO (Entrez) format to  GRCh38 Ensembl ids (specifically to Ensembl release 100), so that I can map the data to Ensembl's annotations. I would like to get suggestions on this.
I feel that figuring out the genome assembly version (GRCh37/38) of the Entrez ids would be critical in this conversion. Although, I am not aware if Entrez ids change from one genome assembly to other.

Hugo names don't really have an associated assembly. They are just the official name of the gene in human. You will find them in all assemblies (except if they hadn't been annotated/identified in one of the older ones). Entrez entries are also not tied to a specific genome assembly. They usually have mappings to multiple ones. For example, Hugo gene name A1BG, with Entrez ID 1 has:

Now, if you want to get Ensembl gene IDs for a list of Hugo gene names, you can use BioMart to convert. First, get the list of target IDs. You can use either Hugo or Entrez. I found the Entrez simpler to parse out, so I went for those (note that this assumes GNU grep):
grep -oP '(?<=\()\d+(?=\))' CRISPR_gene_effect.csv > entrez.ids

I then went to BioMart and selected "Ensembl genes 104" (note that I could only get Ensembl version 104, I assume you're not adamant about getting version 100 given that it doesn't seem to be current) and "Human genes GRCh38.p13" :

Next, I set the filters to "NCBI gene (formerly Entrezgene) ID(s)" and uploaded my file:

And, finally, selected the Attributes I want to see, namely Ensembl gene ID, HGNC (Hugo) symbol and Entrez ID so I can map everything to each other:

I then clicked on "Results" and got:
[
You can now download this to a file and have your mapping. You might even be able to use the one I created: https://www.ensembl.org/biomart/martview/352aae8e0b5c774cbbee73d67634d577.

Your answer is very helpful. Many thanks!