text:
Odds ratio and enrichment of SNPs in gene regions?
I did a QTL analysis with a panel of 7M SNPs, and want to analyze the enrichment of the significant qtl-SNPs in different genic regions (promoters, gene bodies, TFBS, etc.).
A straightforward way to do it would be with an odds ratio (OR), using as 'background' the full 7M SNPs, but I prefer to compare with a background with a similar MAF distribution to that of the qtl-SNPs. For that, I'm thinking in two alternatives:

Take a single random sample of all SNPs with size 10n (n is the number of qtl-SNPs), keeping the same MAF distribution of the qtl-SNPs results, as well as the same proportion of qtl-SNPs in the sample compared to the full 7M panel. Then, test enrichments using ONLY the random sample (with OR and Fisher's exact test).
Do a bootstrapping, taking multiple random samples (excluding qtl-SNPs), each of size n, and keeping as well the same MAF distribution. For each iteration, get the OR between the full qtl-SNPs and the random sample. Then, the mean of the OR will be the point estimate and the quantiles the confidence intervals. For this approach, the problem would be to compute a p-value.

Which approach seems to fit better? If the second one, how to quantify the significance? Can you suggest a better approach?
Thanks.

also posted in https://stats.stackexchange.com/questions/485997/odds-ratio-and-enrichment-of-snps-in-gene-regions. Avoid posting the same question twice. I would suggest thinking about the nature of  your question which forum is more suitable

@StupidWolf the MAF distribution of the full 7M and the qtl-SNPs is different. That's why I don't want to consider the full panel and instead use a more comparable 'background' of SNPs to have a meaningful estimation of the enrichment.

The first approach only address the question of how likely are you to end up with the observed over-representation given the MAF distribution.
My suggestion is to use the second approach, but I am not sure if you would call it bootstrap. Bootstrap in general means sampling with replacement to estimate the uncertainty of a parameter, in this case, OR. So even if you resample your full dataset, you get an estimate of how variable your enrichment is, not how likely you are to observe it by chance. so just name it permutation or sampling.
This is what I think might work

Calculate your true OR, qtl snps versus all others for say TFBS
Divide your snps according to MAF into 10 bins / stratas
Within each bin, we permute the label whether the snp is a qtl
Calculate OR of selection / others for TFBS
Repeat 1-3 N times to collect the permutated OR
p value would be number of (permutated OR > observed OR+1)/(N+1)

Below is an example implementation in R, first we have 1e5 signficant snps and they are 5 times as likely to be in TFBS:
snp_is_qtl = rep(c(1,0),c(0.1e6,6.9e6))
snp_is_TFBS = c(sample(0:1,0.1e6,replace=TRUE,prob=c(0.95,0.05)),
sample(0:1,6.9e6,replace=TRUE,prob=c(0.99,0.01)))

We have the MAF for all snps, with those in qtl being different from others and put them in bins.
snp_MAF = c(runif(0.1e6,0.4,0.8),runif(6.9e6,0.1,0.9))

We sort the snps according to MAF to make it easier:
o = order(snp_MAF)
snp_is_qtl = snp_is_qtl[o]
snp_is_TFBS = snp_is_TFBS[o]
snp_MAF = snp_MAF[o]

And we create the bin:
bin_MAF = cut(snp_MAF,10)

Lastly a function to calculate OR:
calOR = function(snps,TBFS){
tab = table(snps,TBFS)
log(tab[1,1])+log(tab[2,2]) -log(tab[1,2]) -log(tab[2,1])
}

So our observed OR (in log) is:
obs_OR = calOR(snp_is_qtl,snp_is_TFBS)
obs_OR
[1] 1.656859

Now we collect the OR under permutation:
N = 100
permutated_OR = sapply(1:100,function(i){
print(i)
perm_qtl = unlist(tapply(snp_is_qtl,bin_MAF,sample))
calOR(perm_qtl,snp_is_TFBS)
})

hist(permutated_OR)

Your observed OR way exceeds the permuted, even after controlling for MAF. Note the computation is a bit heavy because of the large number of snps. So it makes sense to calculate all the ORs (i.e TFBS, promoters..) with each permutation


text:
How can I get or create a reference genome for Bacteria?
I am a computer engineer and nowadays trying to grasp some concepts of Bioinformatics particularly, reference genomes and genomic variants. My aim is to find the effect of sequence features on variant types that occur in bacteria genomes. I chose Escherichia coli as the bacteria for my study since there seem to be more publication and data about it and its genomic sequences.
When I made a research about variants, I found out that we need to have a reference genome for the comparison with other alleles in order to find variants. And, I even read that we can use genomic sequence of one bacteria as the reference but it won't represent the cohort in an accurate manner. I got the genomic sequence data for several strains of Escherichia coli bacteria that been sequenced before. However, most of these sequences are not complete ones. I filtered out the complete DNA genomes. When I checked the reports on these genomes, I learned that these are complete genomes but not reference genomes. Only the ones which have RefSeq category: Reference Genome in their report are reference genomes for Escherichia coli.
Learning that most of the data I got for different strains of Escherichia coli were not reference genomes, I tried to seek to find reference genomes and I got this website of NCBI ref sequences for prokaryotes. This website shows only two reference sequences which match with what I have find in all sequencing data I got before. So, I want to ask two questions
First: How can I get a reference sequence for other strains of E. coli?
Second: If I cannot get a reference sequence for other strains, can I use complete genome data as the reference one? When I check the complete genome data, I cannot understand their order because there are other data for same strains that show loci of sequences but in the complete genome, I cannot see any order of the loci. Even, I cannot see data of some loci in the complete genome which blows my mind more.
I know it may not be appropriate to ask two questions in a single post, however, these two questions seem to support each other. I will appreciate any of your help.

How can I get a reference sequence for other strains of E. coli?

You can use any sequence you want as reference for a variation analysis. It just depends on the question you want to answer.
So to keep it simple, your Question is:
What are the variant region from a fresh sequenced E. coli K12 sample?
To know the differences from your known strain E. coli K12. A reason for this could be that you work with the strain K12 through out several years. During your studies you changed the conditions for the bacteria, like adding antibiotics or change the growth media. So some of the E. coli would adapt to the antibiotic and become resistant to it or they adapt better to the new media. This happens because some nucleotide changed. With each new experiment you used the best adapted colonies from the previouse one.
So you sequence your bacteria and compare it with the E. coli K12 strain. In this case this is your reference. If you use another strain, like SE15 or SE11, as reference, your result would not give you any answer for your question.
If your used data and publications gives you the information which strain were used, then use the given strain as reference.

in genomics, organisms have a single reference genome which is updated( or not) over a period of time. E.coli has different strains and this means there might be reference genomes for each of the strain. If there is no reference genome available the its likely that in constructing the sequence of those strains they might have used any of those two as their reference. This information you can find by reading the paper that published those strains. Therefore I will advice you pick any of those two available reference genome as your reference and work with and make sure you indicate which one you used in your report or document.

source: Stack Exchange

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

sort: votes


text:
Standard Way to Preprocess Gene Expression?
I am trying to collect gene expression data for the point of fitting gene regulatory networks. My background is primarily in computer science and so I am finding the biological literature a bit difficult to penetrate.
Specifically, I want to fit two regulatory networks for several cancers in TCGA from tumor and normal samples. I have queries using the R package TCGABiolinks and was able to follow the tutorial here up to the normalization part.
So I have three questions;

Is there a standard way to normalize gene expression counts? If so, what is it, and is it available as a library in R, python or some other language?
If not, can I just center and scale (z-transform) the samples? I am not interested in differential expression, only quantifying expression from counts.

Question is, what do you mean by 'quantifying expression from counts'? TCGA has a standard normalization available, which is FPKM. It normalizes expression by fragment-size and library depth. Some people prefer TPM to FPKM though (it alters the order of operations), so you may want to download the count data and do a TPM normalization using edgeR.  For visualization, for TCGA, I personally tend to log2(FPKM) and row-center for heatmaps.

okay so is the FPKM normalization what is given in the scaled_estimate design matrix?

Not sure what that matrix is... in TCGAbiolinks or any other way to download from the GDC, you'll have FPKM. It is not scaled though, so probably no.

First off - do not use FPKM for across-sample normalization. This is such a common issue that I have a prepared text for it:

Why not to use FPKM

They have a poor sensitivity when looking at low expressed genes [1, 2]. Short transcripts/genes may have falsely over-estimated expression values [5].
They may not normalise adequately across samples [3]. The reason for this is: FPKM are normalised across libraries using a single variable, namely the “total number of reads”. This library size factor estimation can easily be biased by highly expressed, highly variable genes, unlike the method used by the other approaches, which model the library size factor based on summary statistics extracted from all genes. In other words as stated in [4]: “The median ratio method is usually quite robust. The biggest problem for library size correction is not actually the genes with low counts, but the genes with very high counts, which have high variance and disproportionately influence estimators like the total sum. This is why the total sum is a notoriously bad estimator for library size.” Or the same by Lior’s Pachter himself [6].
They have a high False Positive Rate in differential expression analysis [1, 2]. A lot of the genes/transcripts deemed to be differentially expressed are not, mostly due to the effects of point 1), i.e. high artefactual variability in the expression and point 2) i.e. inadequate library size estimation.

On the other hand, approaches such as DESeq, DESeq2, vst or voom have adequate between-sample normalisation, have a controlled FDR and can reliably determine differential expression even for lowly expressed genes if there is enough replication or if an adequate log2 fold change cutoff is selected [7]. If you really insist on using cufflinks, consider performing upper quartile normalisation, which corrects for the point 1) issue.

References
[1] Dillies M.A., Briefings in Bioinf. 2012
[2] Soneson & Delorenzi, BMC Bioinf. 2013
[3] Lior Pachter’s lecture at CSHL (Cold Spring Harbour Laboratory) from minutes 32 on: https://www.youtube.com/watch?v=5NiFibnbE8o
[4] From Answer: "Over-correction" in the size-factors of the DESeq2 package:
[5] http://www.biomedcentral.com/1471-2105/14/370
[6] https://liorpachter.wordpress.com/2014/04/30/estimating-number-of-transcripts-from-rna-seq-measurements-and-why-i-believe-in-paywall/
[7] http://rnajournal.cshlp.org/content/22/6/839

What I recommend, is that you get the raw count data from TCGA, create a count matrix using the DESeq2 package in R, use the function varianceStabilizingTransformation() from DESeq2 to further process your counts and model the networks from there. Your pipeline already mentions voom, which is another good alternative, as is limma. All of these are available on Bioconductor.

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


text:
Defining methylation level of CpG island
Differentially methylated site is technically measured by each probe. However, CpG island within gene promoter is comprised of numerous probes. Do we have a consensus of how to determine the methylation level of CpG island using methylation level of probes within corresponding CpG island?

In my experience, the most common approach to find DMR is to use the DSS R package to call and merge DML into regions.
DMC/DML (Differential Methylation Cytosine/Loci) is often calculated by comparing the counts of methylated and unmethylated cytosine at one position between samples.
Steps are to first call DMLs with the callDML() function. Then call DMRs from these first results with callDMR(). Finally annotate the DMLs or DMRs using a GFF file. From there you can look for DMR falling into CpG island. Some of the parameters that can be tuned are the minimum region length, minimum number of differential CpG sites, percentage of CpG sites required to be significant. However I have never heard of what would be the standard values for these parameters. It may not be clearly established yet.
So because the process of calling DMLs and DMRs need robust statistical testing, it is suggested to first call the differential methylated states and then annotate it.

Yes, most current studies use DMR to measure methylation level. However, DNA methylation impose its impact by suppressing gene expression. I think it is more reasonable to measure the magnitude of methylation impact on gene expression than merely defining statistically methylated region. To do this, we may need better understanding on the correlation between methylation and mRNA expression (things like whether different CpG site have equal contribution on the expression of corresponding gene)

Well as you may know the effect of methylation on genes expression is far from being clear, most studies showing poor negative correlation between differential gene expression and methylation. An interesting approach was the one of Vanderkraats in https://doi.org/10.1093/nar/gkt482. The main idea was to look for pattern of differential methylation around genes promotor that could explain genes expression.

Thanks @PaulEndymion a detailed adn useful response


text:
Somatic mutations for normal WES samples
I am trying to identify somatic mutations on healthy subjects (specifically in WES samples). As of now, I have found a couple of papers who have done similar studies such as this paper where somatic mutations in healthy subjects are identified (using RNA-SEQ data). But in my case, I want to find somatic mutations from normal WES samples. For somatic mutations, we need both tumors and matched normal samples to identify somatic mutations but how to identify somatic mutations from normal samples only. Please suggest. Thanks.

You can follow the GATK pipeline

@csgroen: As I have mentioned in my question, I have to find somatic mutations in normal samples only. The link you have shared contains the GATK pipeline for somatic mutations with tumor and matched normal samples. I am already aware of this. Thanks.

Hm, I see... The issue is that you need your 'baseline' germline genome to be able to call somatic mutations. Otherwise, how can you know if they are germline or somatic?

Yes. Precisely. Or I just found one way in this paper (Whole-genome DNA sequencing provides an atlas of somatic mutagenesis in healthy
human cells and identifies a tumor-prone cell type) in which conditions are not clear to me.

From a quick read, my interpretation is that they aligned against a reference genome (hg19), then removed annotated known germline SNPs (using dbSNP as reference)

check out this one:https://gatk.broadinstitute.org/hc/en-us/articles/360035535932-Germline-short-variant-discovery-SNPs-Indels-, this is for germline mutations

@PhoenixMu he wants somatic, not germline mutations. That's actually the challenge, as it is normal tissue and doesn't have a proper reference. But personally, I think what I described from my interpretation of the paper you shared seems to be at least on the correct path. Align against reference genome and remove common SNPs. It's not perfect, there will be false postiives and negatives, but it should work overall.

@csgroen: Yes. I agree with you. But my problem is the result's accuracy. The pipeline used in the above-mentioned paper considered 3 variant callers and allowing only variants with low allele fraction. How low value of allele fraction gives a better solution?? And 3rd condition "variants that occurred during the in vitro culture of
the clone" is also unclear to me.

Paper quote: "To eliminate all other germline variants, including rare ones, we utilized the low confidence germline variants called by GTEx. [...] Our goal was to exclude as many germline variants as possible, so we reasoned that using all germline variants called by GTEx—including the low-quality ones—was the safest option to minimize germline mutation contamination in our somatic mutation calls." I believe they removed all annotated variants, even if they were low confidence. I don't understand your in vitro comment though, makes no sense for this study.


text:
GTF parser for Python
I found this GTF parser for Python. However, it does not work with Scallop's GTF file:
chr01_pilon_pilon       scallop transcript      168145  169166  1000    .       .       gene_id "gene.1.0"; transcript_id "gene.1.0.0"; RPKM "9.8016"; cov "27.2554";
chr01_pilon_pilon       scallop exon    168145  169166  1000    .       .       gene_id "gene.1.0"; transcript_id "gene.1.0.0"; exon "1";
chr01_pilon_pilon       scallop transcript      410672  444190  1000    +       .       gene_id "gene.3.0"; transcript_id "gene.3.0.1"; RPKM "0.5936"; cov "1.6506";
chr01_pilon_pilon       scallop exon    410672  411205  1000    +       .       gene_id "gene.3.0"; transcript_id "gene.3.0.1"; exon "1";
chr01_pilon_pilon       scallop exon    411326  411482  1000    +       .       gene_id "gene.3.0"; transcript_id "gene.3.0.1"; exon "2";
chr01_pilon_pilon       scallop exon    444085  444190  1000    +       .       gene_id "gene.3.0"; transcript_id "gene.3.0.1"; exon "3";
chr01_pilon_pilon       scallop transcript      678431  681781  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.1"; RPKM "0.4733"; cov "1.3162";
chr01_pilon_pilon       scallop exon    678431  678822  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.1"; exon "1";
chr01_pilon_pilon       scallop exon    678917  679068  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.1"; exon "2";
chr01_pilon_pilon       scallop exon    680616  681781  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.1"; exon "3";
chr01_pilon_pilon       scallop transcript      678431  685851  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.2"; RPKM "7.9200"; cov "22.0230";
chr01_pilon_pilon       scallop exon    678431  678822  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.2"; exon "1";
chr01_pilon_pilon       scallop exon    678917  679068  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.2"; exon "2";
chr01_pilon_pilon       scallop exon    680616  681627  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.2"; exon "3";
chr01_pilon_pilon       scallop exon    681874  681965  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.2"; exon "4";
chr01_pilon_pilon       scallop exon    683200  683274  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.2"; exon "5";
chr01_pilon_pilon       scallop exon    683744  683805  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.2"; exon "6";
chr01_pilon_pilon       scallop exon    684640  684709  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.2"; exon "7";
chr01_pilon_pilon       scallop exon    684842  685851  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.2"; exon "8";
chr01_pilon_pilon       scallop transcript      681834  685851  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.3"; RPKM "1.1177"; cov "3.1080";
chr01_pilon_pilon       scallop exon    681834  681965  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.3"; exon "1";
chr01_pilon_pilon       scallop exon    683200  683274  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.3"; exon "2";
chr01_pilon_pilon       scallop exon    683744  683805  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.3"; exon "3";
chr01_pilon_pilon       scallop exon    684640  684709  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.3"; exon "4";
chr01_pilon_pilon       scallop exon    684842  685851  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.3"; exon "5";
chr01_pilon_pilon       scallop transcript      683032  685851  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.0"; RPKM "0.6001"; cov "1.6686";
chr01_pilon_pilon       scallop exon    683032  683274  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.0"; exon "1";
chr01_pilon_pilon       scallop exon    683744  683805  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.0"; exon "2";
chr01_pilon_pilon       scallop exon    684640  684709  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.0"; exon "3";
chr01_pilon_pilon       scallop exon    684842  685851  1000    +       .       gene_id "gene.6.0"; transcript_id "gene.6.0.0"; exon "4";

However, the below code caused the following errors:
Traceback (most recent call last):
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/gtfparse/read_gtf.py", line 103, in parse_gtf
    for df in chunk_iterator:
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/pandas/io/parsers.py", line 1128, in __next__
    return self.get_chunk()
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/pandas/io/parsers.py", line 1188, in get_chunk
    return self.read(nrows=size)
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/pandas/io/parsers.py", line 1154, in read
    ret = self._engine.read(nrows)
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/pandas/io/parsers.py", line 2059, in read
    data = self._reader.read(nrows)
  File "pandas/_libs/parsers.pyx", line 881, in pandas._libs.parsers.TextReader.read
  File "pandas/_libs/parsers.pyx", line 908, in pandas._libs.parsers.TextReader._read_low_memory
  File "pandas/_libs/parsers.pyx", line 973, in pandas._libs.parsers.TextReader._read_rows
  File "pandas/_libs/parsers.pyx", line 1105, in pandas._libs.parsers.TextReader._convert_column_data
  File "pandas/_libs/parsers.pyx", line 1136, in pandas._libs.parsers.TextReader._convert_tokens
  File "pandas/_libs/parsers.pyx", line 1232, in pandas._libs.parsers.TextReader._convert_with_dtype
ValueError: Integer column has NA values in column 3

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/Users/lorencm/projects/bioinf-scripts/reduceIsoforms.py", line 24, in <module>
    df = read_gtf("data/hybrid.gtf")
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/gtfparse/read_gtf.py", line 211, in read_gtf
    restrict_attribute_columns=usecols)
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/gtfparse/read_gtf.py", line 154, in parse_gtf_and_expand_attributes
    features=features)
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/gtfparse/read_gtf.py", line 122, in parse_gtf
    raise ParsingError(str(e))
gtfparse.parsing_error.ParsingError: Integer column has NA values in column 3

This is the code:
df = read_gtf("data/hybrid.gtf")
df_genes = df[df["feature"] == "gene"]
print(df_genes)

Thank you in advance,

Your GTF file contains NA values in the 3rd column. You can use a subset of the GTF file, for example what you have copy pasted above (seems like there are no NAs) and re-run your script to see if it is the GTF file or the script that is causing the error. I can't think of a reason for NAs in the coordinate columns though!

Thank you, I tried a different GTF file and this error disappeared. However, I got a new one

Do you have a link to the gtf-file? import pyranges as pr; pr.read_gtf("file.gtf") might work.


text:
Optitype for Singularity
I am trying to use optitype on an HPC without docker support or sudo privileges. A search brought me https://github.com/researchapps/optitype and https://github.com/nf-core/hlatyping. I tried bootstrapping the singularity image as follows:
singularity build optitype.simg docker://fred2/optitype

However, this does not seem to run as expected. Could someone with experience bootstrapping docker files onto singularity suggest how I could change this? How do the commands change when doing a bootstrap? I would appreciate your help!

'Not as expected' is a little vague. Alternatively try to build from https://hub.docker.com/r/nfcore/hlatyping

Why don't you talk to the HPC admins and ask for advice. Usually there are best practices on how to deal with Docker images. The admins know the system better than any of us strangers on the internet.

Cross-posted: https://www.biostars.org/p/460428/

Have you just tried the singularity container? https://biocontainers.pro/#/tools/optitype

Bit of an old post sorry, but I think you might be getting a little confused here. The nf-core pipelines come with their own software already, so you don't need to install or bootstrap optitype yourself. If it's used in the pipeline, it comes in the pipeline container(s) (or conda environment or whatever packaging system you use). So you should never have to touch docker://fred2/optitype.
When you run the pipeline with the -profile singularity argument, it will automatically pull the image from Docker Hub and convert to Singularity for you (some of the newer pipelines pull pre-build Singularity images directly). You can also use the nf-core download hlatyping --singularity command from the helper tools to pull the image ahead of time.

Docker images for this pipeline: https://hub.docker.com/r/nfcore/hlatyping
Documentation for this pipeline: https://nf-co.re/hlatyping/
nf-core/tools download docs: https://nf-co.re/tools/#downloading-pipelines-for-offline-use

I hope this helps! Come and join the nf-core Slack organisation if you need more help. We have a channel for each pipeline, including #hlatyping.


text:
Does NCBI's blast API block my IP?
I am trying to run a blatn command from both my laptop and within a google colab notebook. I am not sure why, but this command runs properly once and then, on the second execution, it gets stuck for hours without a reply.
Is it possible that NCBI blocks my IP? I've tried also to run web_blast.pl and this didn't work either.
biopython:
from Bio.Blast import NCBIWWW

fasta = '>seq0\n'+'AACCACCTGACGGGCTCTGT'
blast_res = NCBIWWW.qblast('blastn', 'nt', fasta)
print(blast_res.read())

web_blast.pl command:
web_blast.pl blastn nt input.fa > blast.out
Note, if this is due to throttling by NCBI, I am not sure how do they expect people to query their web service if they block you after one call.

If you have a high usage department might explain it. Try via your own ip

You should check the User Guide:

The NCBI BLAST servers are a shared resource. We give priority to
interactive users. In order to ensure availability of the service to
the entire community, we may limit searches for some high volume
users. Interactive users of the NCBI webpages through a web browser
should not encounter problems. We will move searches of users who
submit more than 100 searches in a 24 hour period to a slower queue,
or, in extreme cases, will block the requests. To avoid problems, API
users should comply with the following guidelines:

Do not contact the server more often than once every 10 seconds.

Do not poll for any single RID more often than once a minute.

Use the URL parameter email and tool, so that the NCBI can contact you
if there is a problem.

Run scripts weekends or between 9 pm and 5 am Eastern time on weekdays
if more than 50 searches will be submitted.


text:
Help importing a biom file into R
How can I import a biom file?
When I run the following code:
library(phyloseq)
import_biom(
    "/Users/name/Dropbox/*****.biom",
    parseFunction=parse_taxonomy_greengenes,
    parallel=F)

I get the following error message:

biom either as JSON (BIOM-v1) or HDF5 (BIOM-v2). Check file path, file
name, file itself, then try agai

what happens when you do file.exists("/Users/name/Dropbox/*****.bio")

you can check with file *****.biom as well, to check if your biom file is text or HDF5-like file, to distinguish between versions 1 and 2.

The relevant package for import_biom should be added to your code (e.g. library(phyloseq)).

Thank you.    file.exists("/Users/name/Dropbox/*****.bio")  [1] TRUE

to check if your biom file is text or HDF5-like file, to distinguish between versions 1 and 2.
How can I check? If it is version 1, what should I do and if it is version 2, what should I do?

Biological Observation Matrix 1.0.0, so I think this is JSON file.


text:
Inheritance annotation for nsSNPs
I am relatively new to SNP analysis. Is there a database to find annotation about the Mendelian inheritance of SNPs?
I have a small list of nsSNPs (non-synonymous SNPs) and I need to find how they are they passed to from parents to offspring.
Thank you in advance.

If you only have a small number of genes, OMIM is a very easy way to do this.
For example:

If you have too many to do manually, I suspect there is an easy way to query lists, although I'm not sure exactly the way myself. Perhaps someone can add a new answer if so.


text:
Error using htseq-count: Could not retrieve index file
I have a total RNAseq dataset that I aligned using STAR producing BAM files (sorted by coordinates). I am now trying to get counts for the lncRNA sequences using htseq-count, with the command:
htseq-count -m intersection-nonempty -f bam -i gene_id -r pos -s reverse ./OVCA429_Aligned.sortedByCoord.out.bam ./GRCh38_latest_genomic.gtf > ./OVCA429.counts.txt

This is the output:
[E::idx_find_and_load] Could not retrieve index file for './OVCA429_587680_1_1.rnaseq.fastqAligned.sortedByCoord.out.bam'
100000 GFF lines processed.
200000 GFF lines processed.
300000 GFF lines processed.
400000 GFF lines processed.
...(truncated output here)...
3775001 GFF lines processed.
[E::idx_find_and_load] Could not retrieve index file for './OVCA429_587680_1_1.rnaseq.fastqAligned.sortedByCoord.out.bam'
100000 alignment records processed.
200000 alignment records processed.
300000 alignment records processed.
...(truncated output here)...

Output text file:
__no_feature    20368875
__ambiguous 0
__too_low_aQual 0
__not_aligned   0
__alignment_not_unique  2634439

I've read through the htseq-count literature but it never asks to include an index file, therefore i'm not sure how to deal with this problem. Any help/ suggestions would be welcome! Thank you, in advance.

The problem seems to be missing index files of the bam files (therefore the error message that the index file is missing). You should be able to generate those with samtools, e.g.
samtools index ./OVCA429_587680_1_1.rnaseq.fastqAligned.sortedByCoord.out.bam

As it happens I did have the index file as well, I just ended up naming it differently and so the application probably didn't recognise it. Once I changed the name to match the bam file, it worked out perfectly file. Nevertheless, thank you for your help! :)

Ah, I did not even think this could be a problem as samtools do not give you an option to give it a different name. GLad you resolved it!

I have created the index using what Kamil suggested but I am not able to generate the counts file. This is what I use: for i in *fastqAligned.sortedByCoord.out.bam; do htseq-count -r pos -f bam $i /gpfs/scratch/mice/GCA-000001635.90-GRCm39-genomic.gtf > $i_output_basename.counts ; done

@christosnoutsos 1. that was not an answer. You should open a new question, link this one and explain why the solution did not work for you. 2. I would put quotation marks around $i in $i_output_basename.counts, i.e. "$i"_output_basename.counts. You probably created a bunch of files .count that replaced each other.


text:
Pointing Cromwell to local EBS mounted storage on AWS Batch
I'm trying to run the WDL tasks from the GTEx RNA-Seq pipeline with Cromwell using AWS Batch as a backend.
I store the STAR alignment index on an Elastic Block Store since my Cromwell server instance lacks sufficient disk space. The EBS Block is mounted to my Cromwell server at /data. The STAR index is located at /data/genomic-references/rnaseq/star_index_oh50.
Try as I might, I cannot get Cromwell to run the STAR task correctly.
My JSON inputs file:
 {
   "rnaseq_pipeline_fastq_workflow.prefix": "s3://mybucket/test/CL100148068",
   "rnaseq_pipeline_fastq_workflow.star.fastq1": "s3://anotherbucket/genotypedata/sample1.fq.gz",
   "rnaseq_pipeline_fastq_workflow.star.star_index": "/data/genomic-references/rnaseq/star_index_oh50/",
   "rnaseq_pipeline_fastq_workflow.star.memory": "32",
   "rnaseq_pipeline_fastq_workflow.star.num_threads": "4",
   "rnaseq_pipeline_fastq_workflow.star.disk_space": "200",
   "rnaseq_pipeline_fastq_workflow.star.num_preempt": "0"
 }

Error message:
tar: /data/genomic-references/rnaseq/star_index_oh50/: Cannot open: No such file or directory

Note that /data/genomic-references/rnaseq/star_index_oh50/ does exist and is readable.
The relevant command in the GTEx container is here:
tar -xvvf ${star_index} -C star_index --strip-components=1

If the files are stored locally, then why does tar fail to wrap the directory?

Somebody helped me with this one through other channels, but I wanted to make it publicly available:
Make a tarball of the STAR index. The option --strip-components=1 strips the first level of the path from the untarred items; thus, the expected tarball was made with something like
tar -cvf ${star_index} the_star_index_directory/

Then point Cromwell/WDL to the tarball, e.g. on S3. The GTEx container will unwrap the tarball on the fly during execution.


text:
median normalization for proteomics
I am using the data from a proteomics study were the data was log2 transformed and then a median normalization was applied. The data was normalized by groups of conditions (normal, mutant), not for all the arrays at the same time, so the median for the normal group is slightly higher than for the mutant group.
I am wondering if this is correct or should I normalized the data en general, I mean taking into account all the arrays and normalizing all at the same time so the 2 groups have the same median.
Regards

The kind of normalization depends on what you what to explore. So, there is no absolute answer to this. Different normalizations highlight different aspects of your dataset. There is a nice paper for a higher throughput quantitative study that I would recommend reading:
Nusinow, David P., John Szpyt, Mahmoud Ghandi, Christopher M. Rose, E. Robert McDonald 3rd, Marian Kalocsay, Judit Jané-Valbuena, et al. 2020. “Quantitative Proteomics of the Cancer Cell Line Encyclopedia.” Cell 180 (2): 387–402.e16.
What you want to aim for is a normalization that removes batch effects and that clusters together your quality control samples that you know are very similar to each other.


text:
TopHat2 versus HISAT2 inner workings
In my intro to bioinformatics course, we mentioned that TopHat2 and HISAT2 will both try to align as many reads as possible to the reference genome (TopHat2 has been superseded by HISAT2). For the reads that were not mapped, these probably cross exon boundaries. To solve this issue, both TopHat2 and HISAT2 chop up the unmapped reads into around 25bp fragments. Then:
"TopHat2 tries to align many fragments from each read using all possible exons as the search space. HISAT2 focuses on one fragment from each read, and once that fragment has been anchored it only searches locally."
My question is about that last part "once that fragment has been anchored, [HISAT2] only searches locally." Isn't it possible that a given fragment will map to many different regions of the genome? So how is it possible that it will be "anchored"? Are we going to anchor it multiple times and see which anchor makes most sense or what exactly?

As far as I know, TopHat2 hasn't been supported for a while now.  It shouldn't be mentioned any more or, at least, you should say that it has been superseded by Hisat2.

@Ray thanks for pointing that out I didn't know this

@Ray but even though TopHat2 is no longer supported I still think it's worthwhile to understand how it works since it would give a better understanding of how this technology developed over time

Oh!  I agree with you about understanding it as there are many papers that have relied on TopHat2.  However (I don't have a reference for you; it's from my memory), there were bugs found in TopHat2 and once Hisat2 was released, I don't think they are actively fixing it.  (Which is totally understandable, given the software maintenance cycle of any program.)  Oh -- it's actually the main site where I saw it!  It says, "Please note that TopHat has entered a low maintenance, low support stage as ... superseded by HISAT2...".  This was 2016...

I certainly don't think you should not mention TopHat2, of course.  But if you're teaching a course, you should mention that as I still see papers published using TopHat2.  No one reads the main page...they just go to the download link, I guess!

I don't believe this statement is true:

To solve this issue, both TopHat2 and HISAT2 chop up the unmapped reads into around 25bp fragments.

TopHat2 did not support spliced alignment directly.  This is why it followed this general approach of chopping the read up into small parts, followed by "contiguous" alignment of the sub-reads.  Hisat2 / Hisat-genotype support spliced alignment directly insofar as they don't pre-partition the reads into small parts, but allow dynamic extensions of alignments across splice junctions (in a way similar to how STAR works).  This is generally a fundamentally better approach and is one reason, among many, why TopHat/TopHat2 have been put into maintenance-only mode and have been officially superseded by Hisat2/genotype.


text:
GCF VS GCA or combination of both for pangenomic studies
I have just started to learn bioinformatics and pangenomics. So if this question seems to you pretty basic then I apologize in advance.
As we know, NCBI for the genome database, there are two kinds of sequences available RefSeq and Genbank. I have read the differences in both RefSeq and Genbank. But I was just curious which one would be preferred for pangenome studies? For example, phylogeny built by RefSeq sequences would be different from GenBank sequences? As I read, RefSeq is already curated, annotated and contaminants are removed. Will it affect the phylogeny that we will build after using any pangenome pipeline?
I am particularly interested in non-synonymous mutations (frameshifts or stop codons inside coding sequences). Do RefSeq curators also remove this kind of unusual errors or mutations? If yes then it would be difficult to study such mutations in RefSeq sequences.
Somebody may say a combination of both Genbank and RefSeq would be good, after removing the duplicate one. Which one should I remove RefSeq or Genbank?.
Regards


text:
Advantages of paired-end sequencing compared to single end
One of the advantages of paired end sequencing over single end is that it doubles the amount of data. Another supposed advantage is that it leads to more accurate reads because if say Read 1 (see picture below) maps to two different regions of the genome, Read 2 can be used to help determine which one of the two regions makes more sense. However, I don't understand how this second advantage pertains only to paired end sequencing. With single-end sequencing, couldn't you just use the read that comes right after read 1 to essentially get a longer net read?

An explanation for the downvote would be appreciated

You seem to think you can put the reads in order and they somehow join up, the reality is there is not order to the reads and it's only by assembling them like a jigsaw puzzle or aligning them to a reference genome that you can infer their location. I didn't downvote

@Chris_Rands thank you for your answer. Hmm, could you elaborate more as to why there is no order? I mean, when we're obtaining the reads, aren't we starting with a primer at one end of the fragment and working our way to the other end to generate the multiple reads that make up that fragment?

You have a picture of an Illumina read, which clearly does not work by 'working your way' from one end all the way through to the other.  Why technology are you talking about that works as you suggest, and why did you post a picture of an Illumina read when it's obviously not what you are talking about?  Why are you talking about multiple reads when you are asserting that there is an inherent ordering of single end reads?

There isn't any "order" that could be assumed from the FASTQ file.  That's why after you align reads, you often have to sort the generated SAM/BAM file.

@Ray so if there is no order to the FASTQ file, then how is it possible to infer which reads are paired together (in the case of paired end sequencing)?

Sorry for the confusion!  My "there is no order in the FASTQ file" comment refers to single-ended sequencing -- a single file.  For paired-end sequencing, the reads are usually in two separate files:  one for forward, the other for reverse.  Each read has a name, and the names are unique within a file.  Between files, they are matched up to form a pair by name.  So, yes, they are ordered in the sense that read called "ABC" is the 100th read in the forward file and it should be the 100th read in the reverse file.  A good program would also double-check the names match, too!

But within either file, there is no meaning between the 99th read and the 100th read (i.e., why it is the 100th read).  Thus, returning to single-ended sequencing (i.e., one file), there is no order that can be assumed.  I think some programs accept paired-end sequencing as a single file so that reads are interleaved.  Since they alternate, you can create two files from that easily, if you wanted to.

@Ray okay thanks this makes sense now

@AnIgnorantWanderer You're welcome!  Just to emphasize what someone else hinted, this refers to Illumina reads and how it (BaseSpace or equivalent programs) changes the images to FASTQ files.  Another sequencing technology might choose another format that is still FASTQ, but maybe order of reads exists.  (I don't know of one off hand, but I won't rule that out that it exists now or in the future.)  Good luck!

Like you mentioned, having read pairs can help with alignment.
Your FASTQ for a single-end sequencing run will look something like this:
@SEQ_ID_0001
AGCTAGCGCGGTTGGCTTAGCGACT
+
!''*((((***+))%%%++)(%%%%
@SEQ_ID_0002
AGGTGGTTGTAGGGAAAAAAGTCTC
+
!''*((((***+))%%%++)(%%%%
...

There is not necessarily a relationship between @SEQ_ID_0001 and @SEQ_ID_0002.
They could come from neigbouring regions in the genome or they could come from different species if that's what your experiment involves.
So aligning @SEQ_ID_0001 doesn't exactly give you information about where SEQ_ID_0002 goes.
@SEQ_ID_0001 could equally map to multiple locations in the genome with the same probability.
How do we decide where it "should" go?
We can randomly distribute multi-mapped reads or do some other fancy stuff, but that usually relies on facts about the data as a whole, and not something specific to the original DNA molecule that sequencing read came from.
If we have paired-end sequencing, then like your diagram shows, you know that read 1 (R1) and read 2 (R2; these have the same @SEQ_ID) will relate to each other in some way.
@SEQ_ID_0001 R1 may be multi-mapped, but if @SEQ_ID_0001 R2 isn't, then we know that R1 should probably be relatively close to R2.
Because of how you prepare your DNA samples before sequencing, you'll have some information about the insert length distribution.
You can use this information with sequence aligners to better identify where R1 should be placed.
You can also do this if both R1 and R2 are multi-mapped; certain pairwise alignments will be more likely than others, even if the marginal alignment probabilities for each read are the same.
Finally, there are other unexpected benefits that you can get from paired-end sequencing.
All of the above applies to just the mapping step, but applications like structural variant detection and chromosome conformation capture really benefit from paired-end sequencing.
If you have single-end reads for structural variant detection, then you will only detect breakpoints if the read overlaps a breakpoint.
But if you have paired-end, you not only have twice the amount of DNA sequenced to look for breakpoints, you can also estimate if there is a breakpoint between the two reads (in the inner distance from your diagram) if the two reads align really far apart from each other, or on separate chromosomes.
The idea is similar for chromosome conformation capture and the ligation junction sites.

In paired-end sequencing, you have an expected fragment length. So if Read 1 maps to two different places, whereas Read 2 only maps to one place. Then you can determine which one is more likely for Read 1 to map to. It cannot be too far away from Read 2 given the expected fragment length.
It is not necessarily true that paired-end sequencing gets twice the amount of data. The total amount of data can be determined when running sequencing. If single-end has 20M reads, I can just stop sequencing at 20M as well for paired-end.


text:
How do I Read a Sam File into a String Using ifstream in C++?
I am trying to read a .sam file into a single string in C++. My end goal is to read this string into a vector, where tabs indicate separations between elements of the vector. After making a vector containing all of the data from the sam file, I will separate the relevant data types (chromosome name, locus, etc.) into their own vectors based on position in the master vector.
The sam file was pre-fetched from the SRA database, split into paired-end reads, and aligned to the reference genome using the following Linux commands:
prefetch --output-directory ./ SRR
fastq-dump --defline-seq '@$sn[_$rn]/$ri' --split-files SRR.sra
bwa mem -M reference_genome.fasta SRR_1.fastq SRR_2.fastq > SRR.sam

To make the file easier to read in C++, I altered the .sam file in Linux by removing the header, replacing newline characters with tabs, and replacing gaps between "columns" with true tabs rather than spaces. Since the number of spaces between each column was variable, I accomplished this last part by replacing successively smaller numbers of spaces (12, 11, 10, etc) with tabs. Here is the Linux command line code for the three steps listed above:
grep -v '^@' SRR.sam > headerless_SRR.sam

# Condense sam files into a single line
sed 's/\n/\t/g' headerless_SRR.sam > condensed_SRR.sam

# Replace all spaces in sam files with tabs. Since some columns have multiple spaces, the spaces must be replaced in progressively smaller chunks
sed 's/            /\t/g' condensed_SRR791885.sam > tch_SRR791885.sam     # "tcd" is tabbed, condensed, headerless
sed 's/           /\t/g' tch_SRR.sam > tch2_SRR791885.sam     # Repeat this line with progressively smaller numbers of spaces until all spaces are replaced with single tabs

To open the resulting file and read it into a string, I created the following C++ script, compiled it using the g++ command in Linux, and ran the resulting .out file using "./"
#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <string>

int main(){
    // Creates an input file stream opbject and fills it with data from the sam file
    std::ifstream sam_file;
    sam_file.open("Sample_Sam.txt", std::ifstream::in);
    // Reads the first character in the sam file and stores it as a "character" object
    char sam_character = sam_file.get();
    // Creates a string into which characters can be pushed
    std::string sam_string;
    // While characters are still detected in the sam file, add the characters to a string and move to the next character in the file
    while(sam_file.good()){
        sam_string.push_back(sam_character);
        sam_character = sam_file.get();
    }
    sam_file.close();
}

When I run the C++ script on a dummy sam file containing random blocks of data, it runs fine and I can even output the resulting string to "std::cout" or another file. However, when I run it on a "tch_SRR" sam file, the "./" file hangs without creating an output or anything. Does anyone know how I could fix this?

Are you sure it hangs, or are you instead noticing how incredibly inefficient it is (at least read a line at a time and don't modify the file. SAM files are formatted how they are for a reason)? 2. This seems more appropriate for the codereview stackexchange.


text:
Setting Contrast for DESeq2 results
This question was also asked on BioStars
I am trying to recreate this heatmap. How would I compare the four variables: Ly49+/- and MOG/MOGSP to get a single heatmap? Thank you for helping me through this.
I need to compare both Ly49+ versus Ly49− and MOG versus MOG plus SP. Column data in at the bottom. The function I used for DESeq2 is,
dds = DESeqDataSetFromMatrix(countData=countData,

colData=colData,

design= ~ cell_treatment)

I am not sure if the design is right. For the results, how should I write the contrasts? I am very confused and would really appreciate any help!
DataFrame with 11 rows and 2 columns

cell treatment

<character> <character>

MOGSP_2 SPtetramer+CD8+TCell MOGSP

MOGSP_3 SPtetramer+CD8+TCell MOGSP

MOG_3 SPtetramer+CD8+TCell MOG

MOG_4 SPtetramer+CD8+TCell MOG

MOG_5 SPtetramer+CD8+TCell MOG

Ly49N_1 Ly49-CD8+TCell MOGSP

Ly49P_1 Ly49+CD8+TCell MOGSP

Ly49N_2 Ly49-CD8+TCell MOGSP

Ly49P_2 Ly49+CD8+TCell MOGSP

Ly49N_3 Ly49-CD8+TCell MOGSP

Ly49P_3 Ly49+CD8+TCell MOGSP

Thank you!!

Cross-posted: https://www.biostars.org/p/461091/

I closed the cross-post at Biostars because @StupidWolf basically asked for the same details as I did there, and there is no need to split information across two communities.

Your table is not very clear. If I understood it correctly, there are three groups of t-cells, SPtetramer+, Ly49- and Ly49+ and within it there is a MOGSP vs MOG comparison. I would set colData like this:
                   name             Ly treatment
1  SPtetramer+CD8+TCell SPtetramerplus     MOGSP
2  SPtetramer+CD8+TCell SPtetramerplus     MOGSP
3  SPtetramer+CD8+TCell SPtetramerplus       MOG
4  SPtetramer+CD8+TCell SPtetramerplus       MOG
5  SPtetramer+CD8+TCell SPtetramerplus       MOG
6        Ly49-CD8+TCell      Ly49minus     MOGSP
7        Ly49+CD8+TCell       Ly49plus     MOGSP
8        Ly49-CD8+TCell      Ly49minus     MOGSP
9        Ly49+CD8+TCell       Ly49plus     MOGSP
10       Ly49-CD8+TCell      Ly49minus     MOGSP
11       Ly49+CD8+TCell       Ly49plus     MOGSP

I generate an example dataset:
counts = matrix(rnbinom(11000,mu=100,size=1),ncol=11)

And in the design you specify these two factors, meaning gene expression can be affected by these two conditions:
dds = DESeqDataSetFromMatrix(counts,colData=colData,~Ly +treatment)
dds = DESeq(dds)

You can compare them like this:
head(results(dds,c("Ly","Ly49plus","Ly49minus")))
log2 fold change (MLE): Ly Ly49plus vs Ly49minus
Wald test p-value: Ly Ly49plus vs Ly49minus
DataFrame with 6 rows and 6 columns
          baseMean     log2FoldChange            lfcSE               stat
         <numeric>          <numeric>        <numeric>          <numeric>
1 91.4458114779915  0.712987165797688 1.28291188260395  0.555756927241586
2 81.7945234843954 -0.418266530066077 1.14378290088875 -0.365686993345568
3 94.4351681904941 -0.775592496878288 1.32965355129901 -0.583304196886903
4 107.484335203298  -1.36694885973712 1.38595911743023 -0.986283680770933
5 61.2205696982841 -0.760219864564472 1.16779888688769 -0.650985262188887
6 95.5908168437137  0.497427543096335 1.12768919469662   0.44110340458672
             pvalue              padj
          <numeric>         <numeric>
1 0.578377034355333 0.977537706056193
2    0.714598652706 0.992321035598301
3 0.559688537757548 0.977537706056193
4 0.323993925580938 0.954772990635094
5 0.515056000552832 0.977537706056193
6 0.659138138867294 0.992321035598301

To compare groups within another factor:
results(dds,c("treatment","MOGSP","MOG"))
log2 fold change (MLE): treatment MOGSP vs MOG
Wald test p-value: treatment MOGSP vs MOG
DataFrame with 6 rows and 6 columns
          baseMean     log2FoldChange            lfcSE               stat
         <numeric>          <numeric>        <numeric>          <numeric>
1 91.4458114779915 -0.875905107403943 1.43007567177694 -0.612488642867122
2 81.7945234843954  0.595588650280833 1.28109080878515  0.464907441530728
3 94.4351681904941 0.0255578400583986 1.48592401497315 0.0171999643325371
4 107.484335203298  -1.46652057431159 1.55285829334496 -0.944400774105796
5 61.2205696982841  -1.27442596494316 1.30489822606005  -0.97664778715433
6 95.5908168437137  -2.35402796630456 1.26608624796739  -1.85929510733078
              pvalue              padj
           <numeric>         <numeric>
1  0.540214509839105 0.944881849859735
2  0.641997741235223 0.973651054534338
3  0.986277090644156 0.994944323499595
4  0.344964886056295  0.87421063696937
5  0.328743552427222                NA
6 0.0629853201285643 0.525889480109579

The colData i used:
structure(list(name = structure(c(3L, 3L, 3L, 3L, 3L, 1L, 2L,
1L, 2L, 1L, 2L), .Label = c("Ly49-CD8+TCell", "Ly49+CD8+TCell",
"SPtetramer+CD8+TCell"), class = "factor"), Ly = c("SPtetramerplus",
"SPtetramerplus", "SPtetramerplus", "SPtetramerplus", "SPtetramerplus",
"Ly49minus", "Ly49plus", "Ly49minus", "Ly49plus", "Ly49minus",
"Ly49plus"), treatment = structure(c(2L, 2L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L), .Label = c("MOG", "MOGSP"), class = "factor")), class = "data.frame", row.names = c(NA,
-11L))

so it can used for multiple variable or cell type comparisons instead of many pairwise comparison which i do normally

it depends on the design but in general yes. In this example, you can see it's necessary because the groups have some overlap

conceptual doubt so if i have to perform multiple pairwise comparison do i need to set reference level for every comparison or i can simply use "results(dds,c("treatment","MOGSP","MOG"))" or the contrast function?

using results(dds,c("treatment","MOGSP","MOG")) takes care of everything. The reference level of the factor does not matter here

will try this..


text:
What programs account for structural alignment of different parts of distant homologs which have significant structural differences?
If there is a need to perform structural alignment of different parts of distant homologs, which program one should use?
Since distant homologs often have significant structural changes, meaning the assumption of a similar physico-chemical environment for aligned residues would be violated.
Any program/tools or suggested reference will be a great help.
Thanks.

You might want to check out TCoffee's Espresso as it does a structure based alignment (requires structures to be given/exist in PDB). However, in my hands it has always performed less than ideally.

TM-align has some options to do this, see for example https://zhanglab.ccmb.med.umich.edu/TM-align/help.html.
TM-align is similar to TM-score but the alignment is purely structural.

You can use progressiveMauve (http://darlinglab.org/mauve/user-guide/progressivemauve.html), or bubbz (https://www.cell.com/iscience/fulltext/S2589-0042%2820%2930409-0)


text:
Details of DESeq2 modeling a batch effect
When correcting my data for a batch effect using removeBatchEffect, some of the gene expression values become negative.
When searching for differentially expressed genes, I do not use the data above, but rather model the batch using deseq2 (design=~Batch + Condition).
However, I started worrying. When DESeq2 introduces a batch into the model, does it allow for negative values "behind the scenes"?
If it does, I do not understand how that can make sense, in the context of RNA-Seq data.

DESeq2 accounts for batch effect by treating it as a covariate when you use a formula that includes the batch variable. removeBatchEffect attempts to remove batch effects, which is why you get data that may not make sense - removing batch effects is not an exact process, IMO. Rest assured, accounting for batch as a covariate wiil not have DESeq2 operating on batch-effect-removed data in the background.

DESeq2 uses the batch information (and everything else in the design) to produce offsets for its GLM. For a background on that please check how linear models work, e.g. using the StatQuest series of statistics videos over at YouTube.
It still operates on the raw counts. The same goes for the normalization factors.
removeBatchEffect fits a linear model to the data including the batch information and then subtracts the batch component from the counts (that is basically the baseline difference).
If you are interested in preserving the integer nature of counts and preserving zeros as smallest values after explicit batch correction you may want to check ComBat-seq() from the sva package. It operates on raw counts and returns batch-corrected raw counts which you could then normalize and calculate CPMs from (if you need batch corrected CPMs with no negative values). I find it useful and prefer it over removeBatchEffect as it avoids the unfortunate negative counts which sometimes messes up plotting scripts that expect zero as the smallest possible value.


text:
How is the GT field in a VCF file defined?
As my question in SO was closed and asked to be posted in this forum, I am posting it here.
I am not from the bioinformatics domain. However, for the sake of analysis, I am trying to pick up certain basics related to the GT field in the VCF file.
I know we have a field called Alleles. May I know under what circumstances GT takes a certain value and how they are called? Can you confirm my understanding?
Ref Alt    GT   Name
A    A     0/0  Homozygous
A    G     0/1  Heterozygous (does 1/0 also mean the same?) What's this 0 and 1 actually?
A   [C,CA] ??    ??
??   ??    1/1  HOM_ALT? How and why?

Can the experts here help me to fill the question marks and also help me understand with the Ref and Alt combinations when a genotype can take a value of 0/0 or 1/0 or 1/1 or 1/2 etc. and what are the names for those values? Like when it is called HOM_ALT etc.?
Any simple explanation for beginner like me (with no background in bioinformatics/biology) can be helpful.

You can get most of the info from this paper. See Fig. 1 and the surrounding text. Quoting from there, "GT, genotype, encodes alleles as numbers: 0 for the reference allele, 1 for the first allele listed in ALT column, 2 for the second allele listed in ALT and so on."
In your case, the reference allele, here a single nucleotide, is A. When the alternate allele is also A, the genotype GT is reference, or 0. There are 2 copies of each allele in the human genome in non-sex chromosomes chr1-chr22, hence 0/0, or homozygous reference or HOM_REF.
When ALT=G, and the GT column is 0/1, this means that you have 1 reference allele (0), and 1 alternate allele (1). This means that you have A on one copy of this locus, and G on another. The convention is write GT field in ascending order, so 0/1 rather than 1/0. This is called heterozygous, or HET.
When ALT=C,CA, the GT is probably 1/2, because there are 2 alternate alleles, and I assume we continue with the same chromosome present in 2 copies. This means there are no reference alleles here at all, only alternate alleles. It is a heterozygous genotype composed of two different ALT alleles, or HET_ALT. Note that it is not enclosed in square brackets in the vcf file format: A <tab> C,CA <tab> 1/2 ....
Finally, these are some examples of HOM_ALT:
A    C     1/1
A    G     1/1
A    CA    1/1

This means that the same ALT allele (either C, or G, or CA) is present in 2 copies. There is no reference allele present. This is called homozygous alternate genotype.
In general, the name homo means the same, and hetero means different, in the context of genotypes.
REFERENCES:
Danecek P, Auton A, Abecasis G, et al. The variant call format and VCFtools. Bioinformatics. 2011;27(15):2156-2158. doi:10.1093/bioinformatics/btr330 : https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3137218/
SEE ALSO:

VCF - Variant Call Format
What Does Genotype ("0/0", "0/1" Or "1/1") In *.Vcf File Represent?
Difference between 0/0 and ./. for genotype in VCF
1/2 in VCF genotype field?

0 and 1 is just a way of coding the reference (0) and alternate (1) allele. These could be A/G, C/T etc. It's just a simplified way of expressing the different alleles.
0/1 and 1/0 functionally mean the same thing (that the individual is a heterozygote) - since the genotype is unphased, the alleles aren't ordered. The / symbol tells you the genotype is unphased. However, in practice, a heterozygous genotype is always written as 0/1 as a matter of convention. 0/0 is referred to as homozygous reference, and 1/1 as homozygous alternative.
0|1 and 1|0 do mean something different, since the pipe symbol | tell us the order of the alleles matters.
Occasionally you will have multi allelic positions, where you have more than one alternate allele, and therefore the field will look like:
Ref   Alt
A   [C,GT]

The extra alt allele is given the number 2, so in this case if an individual had the genotype CT then their genotype code would be 0/2.
Full detail is given in the official VCF specification page.


text:
Installing HiC-Pro issues
I am setting up a HiChIPseq pipeline, but installing software necessary for the HiC side of things has been unbecoming.
Both my PI and I have tried to install Hi-C Pro onto our Mac Pro that we use as our compute server (though in a couple of weeks we hope to migrate to an EPYC 64 Core computer server of our own!).
There are three approaches that I have taken to install Hi-C Pro.

From Source

To install from source, you need to supply the paths to its dependencies. Believing we did a good job, we run it and find problems accessing bx-python.
No matter how many times we update pip and try to install bx-python and specify that we want it installed for python2.

From singularity

This is the one I tried. From what I understand this creates a VM with the necessary dependencies installed, and allows you to use the program without building. However, it appears to abort when conda tries to install numpy.
+ conda install -y -c anaconda numpy=1.16.3
Collecting package metadata (current_repodata.json): done
Solving environment: failed with initial frozen solve. Retrying with flexible solve.
Collecting package metadata (repodata.json): / Killed
ABORT: Aborting with RETVAL=255
Cleaning up...

I have seen the failed solving environment bug many times, but have no idea how to solve it. I have tried updating conda, but that doesn't solve it.

Install Someone's Premade Conda Environment

Running conda.sh give me this:
Solving environment: failed

ResolvePackageNotFound:
<Lots of packages...>

Sorry if this is sloppy.

Is the conda environment running under python2 or 3?

@ATpoint python2

If you're reading this, I can only assume that's because you had trouble installing HiC-Pro too. I was finally able to install from source and run the software, so here is how I did it, giving detail to the error that I was hit with and the clues I used to solve the puzzles and get through installing this program:

From Source To install from source, you need to supply the paths to
its dependencies. Believing we did a good job, we run it and find
problems accessing bx-python.
No matter how many times we update pip and try to install bx-python
and specify that we want it installed for p

Setup
I downloaded from source once again using git clone.
The I cated the config file to find the programs I needed to find the path of, which I did using which. I ensured that I was using python2.
I also used pip2 to install pysam, numpy, and scipy. I needed to use sudo to install pysam, though. I also installed something called bx-python.
Make Configure
Here is the error I get when I run make configure
Checking Python libraries ...
Traceback (most recent call last):
File "../scripts/install/check_pythonlib.py", line 24, in <module>
      import bx.intervals
  ModuleNotFoundError: No module named 'bx'

  Traceback (most recent call last):
    File "../scripts/install/check_pythonlib.py", line 28, in <module>
      raise ImportError('Error - bx-python cannot be imported')
ImportError: Error - bx-python cannot be imported
Can not proceed without the required Python libraries, please install them and re-run
make[1]: *** [configure] Error 1
make: *** [configure] Error 2

Looking at this error for the fifth time, I realize that the first part, and not the last part, is the important one. It's yelling at me that I don't have the bx module installed. It also gives me the path to the file that threw the exception, so I read it and saw that this was what was being run:
## bx
try:
 import bx.intervals
 if vcmp(bx.__version__, '0.5.0') == -1:
   raise ValueError('bx-python '+bx.__version__+' detected. Version >= 0.5.0 required')
 except ImportError:
   raise ImportError('Error - bx-python cannot be imported')

I notice that there's nothing fancy, the thing just needs to be imported, and its actual name is bx. Wonderful, I pip2 it. But... I hit another error in pip:
pyrsistent requires Python '>=3.5' but the running Python is 2.7.16

Apparently, pyrsistent dropped support for python2. Wonderful. I get around it by:
sudo pip2 install awsebcli pyrsistent==0.16.0

In the middle of this installation, Im reading the output and see:
bx-python 0.8.9 has requirement six>=1.13.0, but you'll have six 1.11.0 which is incompatible.
So:
sudo pip2 install six==1.13.0
At this point, I, without any more problems, sudo pip2 install bx.
I run make configure again, and I run into...
Checking Python libraries ...
Traceback (most recent call last):
File "../scripts/install/check_pythonlib.py", line 24, in <module>
      import bx.intervals
  ModuleNotFoundError: No module named 'bx'

  Traceback (most recent call last):
    File "../scripts/install/check_pythonlib.py", line 28, in <module>
      raise ImportError('Error - bx-python cannot be imported')
ImportError: Error - bx-python cannot be imported
Can not proceed without the required Python libraries, please install them and re-run
make[1]: *** [configure] Error 1
make: *** [configure] Error 2

huh? Didn't I fix this?? Something seems fishy. For whatever reason, my mind then wanders to the Makefile to see if anything there pops out, cause... what else would be calling that python file? In the main Makefile I see:
make -f ./scripts/install/Makefile CONFIG_SYS=$(CONFIG_SYS)

So I head on over to ./scripts/install/Makefile and see:
./scripts/install/install_dependencies.sh -c $(CONFIG_SYS) -o  $(realpath $(PREFIX))/HiC-Proo    _$(VNUM) -q

Nothing like a goose chase. So I head on over to install_dependencies.sh, where I see:
# python
which python > /dev/null;
if [ $? != "0" ]; then
     echo -e "$RED""Can not proceed without Python, please install and re-run""$NORMAL"
     exit 1;
else
     pver=`python --version 2>&1 | cut -d" " -f2` #<- THIS LINE
     vercomp $pver "2.7.0"
     if [[ $? == 2 ]]; then
         echo -e "$RED""Python v2.7.0 or higher is needed [$pver detected].""$NORMAL"
         exit 1;
     fi
fi

I check that lines output, AND IT'S PYTHON3!! For some reason, it didn't give me that error message, it just kept going from that line.
So, I changed all instances of python with python2 in that file.
I reran make configure and got The required Python libraries appear to be already installed. Bliss.
Make Install
Next step was make install, and guess what? A roadblock:
Traceback (most recent call last):
File "setup.py", line 2, in <module>
from numpy.distutils.core import setup
  ModuleNotFoundError: No module named 'numpy'
make: *** [iced] Error 1

Hmm... well, it's probably the same thing. python3 being called where python2 should be. So I double checked the main Makefile and changed the setup.py call to use python2.
I run make install and it runs for awhile, chugging along using the g++ compiler. Then I hit an odd error:
cp: /Users/bosslab/Downloads/HiC-Pro-master/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4: name too long (not copied)
I'm guessing that the Makefile has a recursive path somewhere. I check the appropriate section and find:
cp -Ri $(MK_PATH) $(INSTALL_PATH)
I simply replace that line with the actual paths, with the install_path being adjacent instead of inside the mk_path.
Then I get hit with:
HiC-Pro installed in /Users/bosslab/Downloads/HiC-Pro-master/HiC-Pro_2.11.4 !
So running $install_dir/HiC-Pro --help returns the help. YAY!
Leftovers
I installed it, but how about when I actually run it?
Error: The 'ice' command is not in your path. Please check where the 'iced' python package has been installed and update your PATH !
Ok:
$ sudo pip2 install iced
Requirement already satisfied: iced in /Users/bosslab/Library/Python/2.7/lib/python/site-packages (0.4.2)

...so, I set up a test to see whether it actually is installed:
$ python2
Python 2.7.16 (default, Nov  9 2019, 05:55:08)
[GCC 4.2.1 Compatible Apple LLVM 11.0.0 (clang-1100.0.32.4) (-macos10.15-objc-s on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> iimport iced
>>> quit()

It's already installed. But then I had an epiphany:
$ python3
Python 3.8.5 (default, Jul 21 2020, 10:48:26)
[Clang 11.0.3 (clang-1103.0.32.62)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> ipomport iced
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
ModuleNotFoundError: No module named 'iced'
>>> quit()

Aha! Somewhere (i don't know where though) python3 is being called. So I just installed that module for python3 with sudo pip3 install iced.
Then I had some trouble with pointing to the right bowtie index (just give the directory where the index exists) and then I hit my next error dealing with HiC-Pro itself:
Pairing of R1 and R2 tags ...
Logs: logs/dixon_2M/mergeSAM.log
make: *** [bowtie_pairing] Error 1

Checking the log file named in the error message:
/usr/local/bin/python /Users/bosslab/Downloads/HiC-Pro-master/HiC-Pro_2.11.4/scripts/mergeSAM. py -q 0 -t -v -f bowtie_results/bwt2/dixon_2M/SRR400264_00_R1_hg19.bwt2merged.bam -r bowtie_re sults/bwt2/dixon_2M/SRR400264_00_R2_hg19.bwt2merged.bam -o bowtie_results/bwt2/dixon_2M/SRR400 264_00_hg19.bwt2pairs.bam
   File "/Users/bosslab/Downloads/HiC-Pro-master/HiC-Pro_2.11.4/scripts/mergeSAM.py", line 26
   print "Usage : python mergeSAM.py"
    ^
SyntaxError: Missing parentheses in call to 'print'. Did you mean print("Usage : python mergeS AM.py")?

Missing parentheses? Oh my god, it's using python3 for some reason again. I grep the installation directory for a file that calls merge_sam.py, and find config-system.txt. I go in there and just need to change the python path from /usr/bin/local/ where python3 is installed to /usr/bin/ where python2 is installed.
Lastly, I get an error during the ICE normalization step that numpy is not installed. Using my big fat brain, I isntall numpy to my python3 and call it a day. Now it all works as intended.

Yikes, this sounds like an awful time troubleshooting. Glad you were able to solve your issue and post a solution here

I could install HiC-Pro using conda when creating a separate environment (which I always do for side-projects), and then pull the conda build of HiC-Pro:
conda create --name HiCPro python=2.7
conda activate HiCPro
conda install -c davebx hicpro

Not sure what conda.sh is that you describe, sometimes tricky to get random code from GitHub running. Give it a go with a standard conda install as I described.


text:
BioJava - Protein to DNA
Using Biojava 5, I am trying to convert a sequence of Amino Acids to a Nucleotide sequence. Any idea how to do this?

Welcome to SE.bioinformatics. What have you tried so far?

Just looked in BioJava documentation - I don't think they have anything to do this. There is http://www.bioinformatics.org/sms2/rev_trans.html, but I need to write code for this - so can't use that tool

"convert a sequence of Amino Acids to a Nucleotide sequence". Like you've mentioned in a comment, this is called "reverse translation". Avoid the word "convert" as much as you can - it is as generic as a word can get in our context. It doesn't look like BioJava has a reverse translate provision - you might need to write some code yourself.

Does this answer your question? How to translate amino acid sequences to Nucleotide sequences

Thanks Chris for that link. Interesting. I ended up using files from NCBI, that worked.

Note that unlike the reverse (nucleotide sequence to protein sequence), which is uniquely defined given a reading frame, the mapping from protein sequence to nucleotide sequence is not unique (it is not a function in the mathematical sense). The method you use to select which of the many (exponential in input sequence length) possible sequences to output depends critically on what the question is.
If the goal is to design a synthetic gene that will produce a desired protein, then you need to optimize it for expression in the organism into which the gene will be inserted, subject to the constraint of the encoded amino acid sequence.
If the goal is to search for nucleotide sequences in one or more genomes that could give rise to a certain polypeptide, there is a much better way than to do reverse translation at all. In particular, most putative protein-coding open reading frames in most sequenced genomes have already been run through forward translation code and the expected amino acid sequence encoded in them can be found in databases. So, this problem is reduced to a substring search in these databases.

Very interesting! Thanks. That is what I ended up doing. Matching up FASTA files to get sequences.


text:
Find intersection among multiple columns and find the intersected elements
file input file
a <- read.csv("Studies_studies.tsv",sep = "\t")

#https://stackoverflow.com/questions/43658982/cross-comparison-of-columns-of-the-same-data-frame
bb <- tcrossprod(table(melt(as.matrix(a))[-1])) * !diag(37)

dd<- setNames(cbind.data.frame(t(combn(names(a), 2)),
                               combn(names(a), 2, function(x) length(intersect(a[, x[1]],a[, x[2]])))),
              c("col1", "col2", "count"))

So far my code.
I can get the count between each pairwise intersection happening between the columns but Im not able to get the elements that is common or intersected between two columns.
For example
This C0026998 C1879321  1310 I do find that there are 1310 elements common between those to ID but I'm not able to find the name of the ID that are column.
So that way if I extrapolate my problem what i would like to see is the pairwise intersection among all the elements and i would like to store the output dataframe.
Any suggestion or help would be really appreciated.
input file dput
structure(list(C0026998 = c(32621177L, 32586363L, 32504186L,
32391628L, 32333156L, 32171069L), C1879321 = c(27276561L, 29286103L,
29988143L, 30076173L, 30024784L, 30892988L), C0023470 = c(27276561L,
29988143L, 30076173L, 30024784L, 26837842L, 29286103L), C0280449 = c(30024784L,
30076173L, 30892988L, 31688884L, 30651561L, 31516032L), C0522631 = c(26055299L,
12801839L, 10516750L, 28825596L, 12652469L, 16203817L), C1292775 = c(28883080L,
10400416L, 22581002L, 1503929L, 9499663L, NA), C0023472 = c(26837842L,
29091516L, 29498925L, 20525995L, 27217448L, 20525993L), C1292773 = c(30135184L,
24440648L, 26273061L, 21508125L, NA, NA), C4545381 = c(30093401L,
29933073L, 29437791L, 31227358L, 30541745L, 28555084L), C4707228 = c(26492932L,
21791467L, 28179273L, 24077845L, 30287390L, 22362038L), C2939461 = c(28588020L,
29988143L, 30076173L, 30510081L, 30024784L, 30963592L), C1336735 = c(28645776L,
28844816L, 31400961L, 30545576L, 29241450L, 28588020L), C2827362 = c(29119847L,
29486661L, 29567772L, 31775485L, 30670826L, NA), C3839741 = c(28494506L,
26492932L, 30923103L, 25715404L, 28556489L, 25006131L), C3839868 = c(28473620L,
25512507L, 31015209L, 29933069L, 29321554L, 21791467L), C4520840 = c(28246192L,
23935018L, 21750082L, 10741700L, 9407724L, NA), C4528668 = c(28645776L,
28588020L, 30510081L, 30409776L, 30093401L, 28844816L), C2825139 = c(29241450L,
30541745L, 31152020L, 24663049L, NA, NA), C1332153 = c(24440648L,
26273061L, NA, NA, NA, NA), C4721505 = c(18623376L, 21618422L,
21685468L, 25196856L, 16985182L, 18353632L), C2826177 = c(27276561L,
29286103L, 25512507L, 29079128L, 28473620L, 26789727L), C4324477 = c(29286103L,
22375971L, 30153096L, 30368038L, 25006131L, 30054307L), C2826176 = c(27276561L,
29988143L, 27959731L, 29286103L, 25512507L, 28588019L), C2827356 = c(27716285L,
29177434L, 30952642L, 30418178L, 27126994L, 29656442L), C3711382 = c(30510081L,
26492932L, 29991558L, 27927766L, 28490572L, 29860938L), C4704767 = c(30555165L,
30055822L, 25540937L, 30135184L, 30545923L, 30361063L), C0027013 = c(22375971L,
28962635L, 22375970L, 29515114L, 30185431L, 23216616L), C0153886 = c(30510081L,
28644114L, 28588020L, 30555165L, 30024784L, 29304833L), C0279094 = c(28644114L,
28588020L, 30510081L, 30555165L, 30076173L, 30024784L), C0279623 = c(27276561L,
30555165L, 28588020L, 30076173L, 30024784L, 30510081L), C1142169 = c(28588020L,
30510081L, 25512507L, 25540937L, 30409776L, 30668890L), C1275661 = c(29660836L,
30963592L, 22579233L, 21880637L, 26772158L, 27121471L), C1292776 = c(30055822L,
30115541L, 10025899L, 31152020L, 26577691L, 23129738L), C1320638 = c(27959731L,
18489989L, 25425683L, 28179273L, 30510079L, 30655375L), C1737261 = c(29988143L,
28645776L, 29933073L, 30545923L, 30361063L, 28336527L), C2826178 = c(25512507L,
29286103L, 28473620L, 25611491L, 21242187L, 19059939L), C4744720 = c(30555165L,
28644114L, 29597002L, 30541745L, 30118897L, 29660836L), relapse = c(32786187L,
32601798L, 32564196L, 32391628L, 32333156L, 32062741L)), row.names = c(NA,
6L), class = "data.frame")

output what i have
dput(head(dd))
structure(list(col1 = structure(c(3L, 3L, 3L, 3L, 3L, 3L), .Label = c("C0023470",
"C0023472", "C0026998", "C0027013", "C0153886", "C0279094", "C0279623",
"C0280449", "C0522631", "C1142169", "C1275661", "C1292773", "C1292775",
"C1292776", "C1320638", "C1332153", "C1336735", "C1737261", "C1879321",
"C2825139", "C2826176", "C2826177", "C2826178", "C2827356", "C2827362",
"C2939461", "C3711382", "C3839741", "C3839868", "C4324477", "C4520840",
"C4528668", "C4545381", "C4704767", "C4707228", "C4721505", "C4744720"
), class = "factor"), col2 = structure(c(18L, 1L, 7L, 8L, 12L,
2L), .Label = c("C0023470", "C0023472", "C0027013", "C0153886",
"C0279094", "C0279623", "C0280449", "C0522631", "C1142169", "C1275661",
"C1292773", "C1292775", "C1292776", "C1320638", "C1332153", "C1336735",
"C1737261", "C1879321", "C2825139", "C2826176", "C2826177", "C2826178",
"C2827356", "C2827362", "C2939461", "C3711382", "C3839741", "C3839868",
"C4324477", "C4520840", "C4528668", "C4545381", "C4704767", "C4707228",
"C4721505", "C4744720", "relapse"), class = "factor"), count = structure(c(1310L,
1310L, 164L, 9L, 6L, 23L), .Dim = 6L)), row.names = c(NA, 6L), class = "data.frame")

from the above output i do know that  C0026998 C1879321  1310 between C0026998 and C1879321 there are 1310 elements are common. Now the next thing i would like to find out is what are the elements that are common between those two IDs same goes for rest of the column which are to be compared pairwise.
The output I'm expecting is this
structure(list(`intersect(a$C0026998, a$C1879321)` = c(32621177L,
32586363L, 32504186L, 32391628L, 32333156L, 32171069L)), row.names = c(NA,
6L), class = "data.frame")

structure(list(`intersect(a$C0026998, a$C0280449)` = c(32171069L,
31841594L, 31770437L, 31753773L, 31688884L, 31516032L)), row.names = c(NA,
6L), class = "data.frame")

So in each column contain one Mesh term where there PMIDS which i want to find if these studies are overlapping or not.So in the above example i have given intersection between first four columns so i can see the IDs which are common. So i want to do the same with all the columns and save the results.

It'll probably help if you give a reproducible example and use dput(head(data)) or dput(data[100,]) or something, since most people won't want to download a random file from google drive.

Since StackExchange is not a code writing service you should try to ask a generic question that people can actually can benefit from. Please provide a representative input with dput and make it clear how the output is expected to look. I personally will never download random stuff from any GoogleDrive just to understand the question. It is the OPs duty to make this clear at first sight.

sure i will do it dput normal i do that but i thought file would be helpful i will update my question thank you for the suggestion

updated the question with dput.

Can you add an example of the output that you want?

okay will add the output what im expecting

I have updated the question

I’m voting to close this question because SE is not a code-writing service.
There is a tendency recently that users post their data and then expect code to solve it.
Since these questions are often highly custom the community does not benefit from it.

A tidyverse answer.
library("tidyverse")

results <- df %>%
  {gtools::combinations(ncol(.), 2, colnames(.))} %>%
  as_tibble(.name_repair="unique") %>%
  rename(col1=1, col2=2) %>%
  rowwise %>%
  mutate(intersection=list(intersect({{df}}[[col1]], {{df}}[[col2]]))) %>%
  unnest_wider(col=intersection, names_repair="unique") %>%
  pivot_longer(starts_with("..."), names_to="col", values_to="val") %>%
  pivot_wider(names_from=c(col1, col2), names_sep="_", values_from=val) %>%
  select(!col)

The first few results.
> results[1:5, 1:5]
# A tibble: 5 x 5
  C0023470_C00234~ C0023470_C00269~ C0023470_C00270~ C0023470_C01538~
             <int>            <int>            <int>            <int>
1         26837842               NA               NA         30024784
2               NA               NA               NA               NA
3               NA               NA               NA               NA
4               NA               NA               NA               NA
5               NA               NA               NA               NA
# ... with 1 more variable: C0023470_C0279094 <int>

let me try and i tired i do end up with this .."New names:

`` -> ...1
`` -> ...2
Error in rename(., col1 = 1, col2 = 2) :
unused arguments (col1 = 1, col2 = 2)"

I think the problem comes this part the renaming part of the column df %>%   {gtools::combinations(ncol(.), 2, colnames(.))} %>% as_tibble(.name_repair="unique") %>%   rename(col1=1, col2=2)

i think the issue is with my R version perhaps i was getting error in my Cent OS where R version 3.6.0 and when I ran the same code in Ubuntu which run R 3.6.3 i get the result as you have shown ..i think it is a version issue.

Thanks for the elegant solution it works..accepted the answer

@ rpolicastro i did a session clean up your code works perfectly fine in both then environment i think there was some issue with the namespace..

There are some other libraries that have a 'rename' function, so from what you said there was probably a namespace conflict. Another way to fix that is to explicitly state which library a function comes from, so for that case dplyr::rename(...).

yes somehow i was able to figure it out...and thanks for the answer and code greatly helpful


text:
get sequences begining with TA
I have a fasta file and the sequences in them are arranged like this:
>chr8:106028475-106028517()
GTaccaggcctagggttgcactcctgtagtctcagcattggg
>chr8:106040430-106040472()
TAACCGGCTCTCCGCTGCCCTCATCTGTCTAGAATTCAGGTA
>chr9:100245605-100245648()
TACCACCTCCCTCACTCTTTCCCAAATTTATGAGTTAGGCTAT

I want to extract only those which begin with TA. so my final output should look like
>chr8:106040430-106040472()
TAACCGGCTCTCCGCTGCCCTCATCTGTCTAGAATTCAGGTA
>chr9:100245605-100245648()
TACCACCTCCCTCACTCTTTCCCAAATTTATGAGTTAGGCTAT

Kindly help.

Hi. People won't write your code for your (especially as a bit of googling leads to some very simple one-liners which do what you need), but they will help you if you show what code you have already tried to write. Thanks.

seqmagick has a handy --seq-pattern-include feature that works similarly to what Kamil did with a regular grep.
$ cat sequences.fasta
>chr8:106028475-106028517()
GTaccaggcctagggttgcactcctgtagtctcagcattggg
>chr8:106040430-106040472()
TAACCGGCTCTCCGCTGCCCTCATCTGTCTAGAATTCAGGTA
>chr9:100245605-100245648()
TACCACCTCCCTCACTCTTTCCCAAATTTATGAGTTAGGCTAT

$ seqmagick convert --seq-pattern-include "^TA" sequences.fasta -
>chr8:106040430-106040472()
TAACCGGCTCTCCGCTGCCCTCATCTGTCTAGAATTCAGGTA
>chr9:100245605-100245648()
TACCACCTCCCTCACTCTTTCCCAAATTTATGAGTTAGGCTAT

It's a great little tool for this sort of thing built around Biopython, though it's been falling behind in compatibility with the latest versions as you can see from some of the open GitHub issues.

Assuming your sequences are always on a single line
grep -B 1 "^TA" seqences.fasta | grep -v "\-\-" > filtered_seqences.fasta

should do the job. Argument -B specifies how many lines before the pattern match should be printed. But there grep separates individual entries with -- that's why I used a second grep to get rid of those.
Ha, apparently there is an undocumented --no-group-separator flag, so this should do as well
grep -B 1 "^TA" seqences.fasta --no-group-separator > filtered_seqences.fasta

GNU grep actually has an undocumented argument --no-group-separator that can do the trick too.


text:
Finding original papers for gene annotations
I am wondering if there's a database to look up the annotation history of an arbitrary gene/ORF, with links to the primary literature for each step of the annotation. I have found that neither NCBI Gene nor NCBI Protein makes this easy.
As a case in point, I was looking at the "TMEM79" gene, which was recently reported as an inhibitor of Wnt signaling: https://pubmed.ncbi.nlm.nih.gov/32924931/. Searching Pubmed for "TMEM79" gave references establishing that this gene has the synonym "Mattin" and has roles in skin allergy and secretion from the Golgi. However, none of the papers I found reference a paper for the origin of the name "TMEM79" itself. Looking at both the gene and the protein in NCBI only gives these same "late" papers, that were evidently published after the annotation "TMEM79" had already been applied to this ORF for some time, and all acted as if this name had simply fallen from the sky without quoting a reference.
Presumably the ORF was first detected by some automated gene scanning program, then subsequently someone (possibly the same group or possibly a completely different one) ran it through a transmembrane prediction program that detected that it was likely an integral membrane protein. However, strangely nobody provides a citation for either of these steps, which I would certainly want to do if I were writing a paper on this gene (as this protein certainly seems never to have been cloned and biochemically characterized in even the most rudimentary manner, computational inferences at the sequence level are the only means to contextualize a role in any biological process).
So I'm wondering, is there a standard database that makes this easy, i.e. given an ORF and an identifier, it returns the literature citation where that identifier was first applied to that ORF? In this case that would clearly point to the transmembrane prediction step because nobody would name an ORF "TMEMxx" unless they thought it was a transmembrane protein. And that would hopefully cite the first identification of the ORF itself, if those were not one and the same paper.


text:
Get Dockerfile of a bioconda package
Bioconda packages have associated containers images. E.g. bioconductor-csaw has the image quay.io/biocontainers/bioconductor-csaw.
I would like to customize the image, in particular, I would like to add rstudio.
It is possible to get the Dockerfile associated to the image of a bioconda package?
How do you suggest to proceed to customize the image?

The Dockerfile is generated by galaxy-utils, so I'm not entirely sure what it looks like. However you should be able to generate such a "mulled" container as documented here. In short, something of the form:
mulled-build build 'bioconductor-csaw=someversion--somebuildnumber' 'someotherpackage=someversion--somebuildnumber' ...

That probably requires setting the CONDA_IMAGE environment variable. Bioconda currently sets that to quay.io/dpryan79/mulled_container:latest.


text:
Choosing Fisher's Exact or Binomial test for overrepresentation in PANTHER
The PANTHER website offers a tool to obtain the GO-based overrepresentation of a gene list (the analyzed list) versus a reference gene list. After entering both these lists, the tool asks to choose between the Fisher's Exact or Binomial tests.
I guess that the convenience for one or another test would depend on the number of genes to analyze. If that is the case, which would be the ideal or advisable rank of gene numbers in both analyzed and reference lists (or their ratio) for each type of test?
(In particular I am working with a list of around 150 versus a reference list of around 2000 genes).

For a guide on using the tests in PANTHER, see https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6519457/.
In answer to this specific question, the defaults are the best choice (Fisher's Exact Test, with FDR correction). Fisher's Exact test gives more accurate statistics for relatively small counts (e.g. when only a few genes in your list are associated with a particular GO term). The FDR is a better correction to use for these types of analysis in general, especially given the large overlap between classes.
The other options are there mostly to support replication of older analyses, e.g. to replicate a published analysis that used the PANTHER over-representation tool prior to 2018, when the binomial with Bonferroni correction was the default.


text:
How to get the uniparc ids for all members of uniref cluster in batch request?
I currently use the 'UniProt website REST API' (example) to collect the UniParc IDs of all the members of the given cluster. Fetching each entry individually is very slow though, so I wanted to ask if there is a way that I can do this request in batches.
What I have already tried:
There are pages mentioning a batch retrieval of entries. As I understood it, all the examples only adress the ID mapping service. Here I have tried adjusting the given example to the following
import urllib.parse
import urllib.request

url = 'https://www.uniprot.org/uploadlists/'

params = {
'from': 'NF100',
'to': 'UPARC',
'format': 'tab',
'query': 'UniRef100_Q8RKD8'
}

data = urllib.parse.urlencode(params)
data = data.encode('utf-8')
req = urllib.request.Request(url, data)
with urllib.request.urlopen(req) as f:
   response = f.read()
print(response.decode('utf-8'))

which resulted in the output
From    To
UniRef100_Q8RKD8        null

I also checked the interface of this service on the website and there I can't seem to map NF100 to UPARC, so I'm assuming that is why I'm not getting a result via the API.
I also know that the UniRef Id contains either a UniParc ID or the UniProtKB accession, so I could use the UniParc ID directly or fetch it for the UniProtKB via the proteins api with the uniparc endpoint, but from what I can see I only get the UniParc ID/entry of one of the members.

Update 17.10.22
It seems that the api changed and the previous solution doesn't work anymore. I updated my solution by using the Uniprot id mapping service. All the following instructions are copied from the documentation page but added here because link answers are discouraged. The documentation page also includes an implementation in python. The basic example is as follows (using commandline tools)
$ curl
    --request POST 'https://rest.uniprot.org/idmapping/run'
    --form 'ids="P99999,P31946"'
    --form 'from="UniProtKB_AC-ID"'
    --form 'to="UniRef100"'

This returns a job id as a response
{"jobId":"27a020f6334184c4eb382111fbcad0e848f40300"}

You then have to poll the api for the status of the job
$ curl -i 'https://rest.uniprot.org/idmapping/status/{jobId}'

Which returns you a status of FINISHED and the url where you can fetch the result in the HTTP Header. Note that since the job might take some time to finish you need to poll this url until you get the status FINISHED.
HTTP/1.1 303
Server: nginx/1.17.7
Location: https://rest.uniprot.org/idmapping/uniref/results/27a020f6334184c4eb382111fbcad0e848f40300
Content-Type: application/json
Access-Control-Allow-Origin: *
...
{"jobStatus":"FINISHED"}

Access the result url with curl
$ curl -s "https://rest.uniprot.org/idmapping/uniref/results/27a020f6334184c4eb382111fbcad0e848f40300"

Which gives you a json dictionary
{ "results": [...], "failedIds": []}

Each entry in the result list is a dictionary similar to the following example
{
    'from': 'P99999'\n\t
    , 'to': {
        'id': 'UniRef100_P99999'
        , 'memberCount': 17
        , 'representativeMember': {
            'uniparcId': 'UPI000013EAA0'
            , ...
        }
        , 'memberIdTypes': [
            'UniProtKB Unreviewed (TrEMBL)'
            , 'UniProtKB Reviewed (Swiss-Prot)'
        ]
        , 'members': [
            'P99999,0'
            , 'Q6WUX8,0'
            , 'P99998,0'
            , 'Q5RFH4,0'
            , 'G4XXL9,1'
            , 'A0A2J8VDX6,1'
            , 'K7CAI8,1'
            , 'A0A663D840,1'
            , 'G4XXM0,1'
            , 'G4XXM2,1'
            , 'G4XXM4,1'
            , 'G1RY04,1'
            , 'G4XXM3,1'
            , 'G4XXM6,1'
            , 'C9JFR7,1'
            , 'A0A2J8VE06,1'
            , 'A0A2J8KZC0,1'
        ]
        , ...
    }
}

The uniparc id of the cluster is given in the representative Member.
If you want the group members directly they are given as a list of strings that contain the accession and the memberIdType index separated with a comma. You would use the index to look up the type of the accession in the memberIdTypes array.
The result also contains information such as the sequence of the representative Member, or the number of organisms in the cluster.
Note: From what I can tell in the UniRef100 Clusters there will only be one Uniparc Id which you will find with the representative Member. If you want to map to UniRef90 or UniRef50 there will be other Uniparc ids in the member list, however the representative Uniparc Id will still only be in the representative Member, so you it's probably safest to parse both.

Old solution
I found a way to batch the mapping from uniref cluster ids to the ids of all members. The url
https://www.uniprot.org/uniref/?query=UniRef100_A0A003&format=tab&columns=id,members
can be used to get the link from one uniref Cluster to it's member ids. Multiple uniref cluster can be fetched at once by connecting the ids with a +OR+ like so https://www.uniprot.org/uniref/?query=UniRef100_P99999+OR+UniRef100_P31946&format=tab&columns=id,members
Because the uniref cluster ids are not stable I ended up using the following syntax, where instead of querying the UniRef cluster by it's id, I query by one of it's members, namely the one that was used to create the current uniref cluster id: https://www.uniprot.org/uniref/?query=*&format=tab&limit=10&columns=id,members&fil=identity:1.0+AND+member:(P99999+OR+P31946)
In the notes section of the 'Retrieve / ID mapping' documentation of the uniprot they talk about fetching chunks up to 20.000 identifiers at a time. I am not sure if I use the service they are talking about though, so I'm keeping my requests at 100 identifiers per request and that seems to work.
One thing I noticed was that fetching results via the uniref endpoint is faster for small sets of identifiers (<20) since the fetch time of the batch request is similar for smaller and larger (<100) sets. So batching only makes sense for more than 20 identifiers at a time. (Disclaimer: I have only done very basic tests on this)
Once I have all the members of a cluster I can use the proteins/uniparc endpoint of the proteins api to fetch the uniparc data. For uniparc accessions (starting with the string 'UPI') I use the GET /uniparc/upi/{upi} endpoint. For the other accession, which should all be Uniprot accessions, I use the GET /uniparc/accession/{accession} endpoint.
This doesn't answer the question of a direct mapping between the uniref ids and uniparc member ids, but it is a workable solution.