text:
Use of "baseMean" in DeSEQ2
The basemean is described as the "mean of normalized counts of all samples". My question is, how would we interpret differences in basemean between genes, pertaining to reliability and such? For example, if two genes that have similar p-values and l2fchanges, but one has a much higher basemean, would that gene's data be more dependable?
Sorry if this is a trivial question!


text:
Where to find rat genome for STAR alignment?
I'm trying to find the rat genome to carry out STAR alignment for several FASTQ files.  Where would I find these files?

source: Stack Exchange

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

sort: votes


text:
Build protein-protein interaction prediction deep learning model
I'm an undergraduate biology student and my thesis is on designing a deep learning architecture to predict whether two proteins interact or not given their primary sequences. I have read some papers with different approaches to the problem including the paper [https://academic.oup.com/bioinformatics/article/35/14/i305/5529260] uses PIPR model which is the state of the art model for this problem in 2019 and recently in 2021 the paper [https://www.biorxiv.org/content/10.1101/2021.01.22.427866v1] uses D-SCRIPT which surpass PIPR with the aspect of the small dataset and has better generalization feature for many model organisms.
For my problem, I have two questions

can you suggest to me the approach to build a model for my thesis. Should I follow these papers to build a model, how my model will be different from the models in the papers?
Which architecture or model recently that the literature considers the most proper for the protein-protein interaction prediction approach?

I very much appreciate any help from everyone!

Albeit to make a full structure... Check out https://github.com/sokrypton/ColabFold it's a mod of the AlphaFold2 notebook. Albeit not mentioned by AlphaFold2, the Rosettafold paper talks about complexes. It's not benchmarked. If two proteins don't interact the replicates will differ.

Here is the post mentioned by @ATPoint: https://www.biostars.org/p/9485815/

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


text:
Extract read names and the associated nucleotides on specific positions from a BAM file (in R)
Let's assume I have a BAM file and several positions that I would like to examine more closely in this alignment. My goal is to find out whether these positions are on the same reads and which nucleotides each read has in these two positions.
So my output should be a table of all reads that occur in the positions and list the associated nucleotides.
The output should look like this (with 980, 1350 and 1400 the first three positions; I made up the read names):

read name
980
1350
1400
...

ReAdN4m1L03
A
T

...

ReAdN4m2L02
A
T
G
...

ReAdN4m3L01

T
G
...

...
...
...
...
...

At the end I would like to have a table that links the nucleotides via the reads.
So far, I have tried to read in the Bam file at the desired positions:
#Positions of interest
pos <- c(980, 1350, 1400)

#path to BAM file
path <- pathToBamFile

#name of reference used for mapping
refname <- "myreference"

#loading the bam file
bamFile <- BamFile(path)
gr <- GRanges(seqnames = refname, ranges = IRanges(start = pos, end = pos))
params <- ScanBamParam(which = gr, what = scanBamWhat())
aln <- scanBam(bamFile, param = params)

With the following code I can access the reads, but how do I get the corresponding nucleotides in the individual positions?
pos.no <- c(1:length(pos))

i <- 1
aln[[pos.no[i]]]$qname

I have already tried to work with the CIGAR strings aln[[i]]$cigar with the sequences aln[[i]]$seq in order to extract to the nucleotides manually, but that led to more and more problems.
Is there any easy way to extract the nucleotides?
I have searched for a long time now on the internet for an easy solution but without success. Hopefully someone will help me here! Many thanks for your time in advance.

Do you have to use R? This is more convenient if you use the pileup interface from pysam (perhaps Rsamtools provides a similar functionality).

@DevonRyan but does pileup give you an output where the nucleotides are linked to the reads?

Yes, it does, at least in non-R interfaces (i.e., python and C).

This answer does it in Python, perhaps there are some ideas there: https://bioinformatics.stackexchange.com/questions/7401/access-base-aligned-to-particular-reference-position


text:
Transcriptomics vs Proteomics for understanding molecular responses
I am planning some experiments to investigate how mono cultures of microbes respond to different nutrient concentrations. I want to use an -omics based approach to understand how, on a molecular level, the microbes utilise these nutrients by looking for up regulation of associated genes. With the idea looking out for upregulation in mind, I was going to go for transcriptomics. But now I am beginning to realise that regulatory changes in genes does not necessarily translate to the molecular utilisation of a nutrient by a cell. To that end, perhaps I would be better using quantitative proteomics to capture changes in the concentration of functional proteins which are actually interacting with the nutrients. However, I worry that by relying on proteomics I might miss important regulatory changes for microbe species where the proteome is less well defined.
I should say that it's an either/or situation, there is only budget for one.
I'd appreciate any thoughts on the best approach.
Thanks.

My gut feeling is that you will be figuring out a lot more causality using the transcriptomic approach, but that might be just because I am a lot more familiar with that type of data. This is a really open-ended question I think, a lot more suited to discussion forums than Q&A sites.

I agree that this is a discussion question. Proteomics is less available but also less sensitive, albeit requiring less corrections analysis-wise. For transcription factors and other low copy number protein transcriptomics is the best option. Bacteria do not go around sequestering mRNA as much as Eukaryotes, which is nice, but some phyla, produce compounds that make them problematic for RNA work (e.g. the osmolyte in Firmicutes spores (something-picolinate)), so proteomics > transcriptomics for Geobacillus etc. Also... there is metabolomics if enzymes are involved (but way less sensitive).

Also, before starting any -omics, check the quality of your reference genome if it is a non-standard species. This sounds unimportant, but if your have an auto-annotated reference genome, you will spend more time fixing it than anything else —e.g. pathway X is missing enzyme Y, but in the operon with WYZ there is a mystery enzyme, which guilty by association would suggest is a candidate, but you have to get more evidence.

@MatteoFerla Thanks for this. When you say quality, are you referring to coverage? The species I'm interested in has a reference genome with 39% coverage, is that likely to be sufficient?

No, I meant the quality of the annotations: the data is meaningless if too many of the genes are genes of unknown function or worse not what the automatic annotation claimed. I am assuming you have full coverage. reference genome with 39% coverage is unusable.


text:
Creating Experimental Components on SBOL
I am developing software to store experiments or studies using SBOL3. In the regular synthetic biology workflow we will find the experimental components Strain and Media. A safe way to implement them is using tyto.SBO.functional_entity as type for both.
Other attempts to implement media (https://github.com/Bioprotocols/paml/blob/main/examples/growth_curve.py)looks like this:
SC_media = sbol3.Component('SC_Media', 'TBD', name='Synthetic Complete Media')

Therefore there is a need for standardise these types on SBOL3 to provide more accurate representations. The Experimental Factor Ontology (EFO) have ontology terms for Strain http://www.ebi.ac.uk/efo/EFO_0005135 and Media http://www.ebi.ac.uk/efo/EFO_0000579 that can be used as type. Then draft SBOL components looks as follow:
strain = sbol3.Component(fj_strain.name(), http://www.ebi.ac.uk/efo/EFO_0005135)
strain.name =fj_strain.name
strain.description = fj_strain.description()

media = sbol3.Component(fj_media.name(), http://www.ebi.ac.uk/efo/EFO_0000579)
media.name =fj_media.name
media.description = fj_media.description()

The problem with this is the lack of readability on the code given the obscure ontology URIs. This problem has been solved for other ontologies using tyto, but EFO is not supported out-of-the-box and the configuration of new ontologies is not working. On tyto repository there is an example:
from tyto import EBIOntologyLookupService, Ontology
KISAO = Ontology(uri='http://www.biomodels.net/kisao/KISAO_FULL#', endpoints=[EBIOntologyLookupService])
KISAO.Gillespie_direct_algorithm
'http://www.biomodels.net/kisao/KISAO#KISAO_0000029'

But I copied and runed that code but it raise this error:
[ERROR] tyto.py 70: Ontology http://www.biomodels.net/kisao/KISAO_FULL# is not available at EBI Ontology Lookup Service
What I would like to do is something like this:
from tyto import EBIOntologyLookupService, Ontology
efo = Ontology(uri='http://www.ebi.ac.uk/efo', endpoints=[EBIOntologyLookupService])
efo.media

But this raises the same error.
How can I add EFO to tyto? What is the URI that should be provided? Should I raise an issue in tyto to ask for EFO out-of-the-box support?
Do you agree with the ontology terms for Strain and Media? Any feedback?

Section 7.9 of the SBOL3 specification already provides recommended terms for experimental elements.  Among these, Section 7.9.2 specifically addresses strain and media, reading in pat:

The roles of reagents may vary in context: for example, arabinose may serve as an inducer or as a media carbon source. As such, contextual role SHOULD be indicated by an NCI Thesaurus (NCIT) term in a role property of the Feature. Examples include:

NCIT:C64356: PositiveControl
NCIT:C48694: Cell
NCIT:C85504: Media
NCIT:C14419: Strain
NCIT:C120268: Inducer

NCIT is one of the ontologies that is built-in for tyto:
tyto.NCIT.Strain
Out[4]: 'https://identifiers.org/ncit:C95501'
tyto.NCIT.Media
Out[5]: 'https://identifiers.org/ncit:C19442'

One thing that makes noise to me is that these terms are exemplified as roles not as types. These terms make sense as types to me.  Should I use another type like SBO_FUNCTIONAL_ENTITY and keep the NCIT terms as roles?

Type is used to indicate the nature of the material. Most of these are complex entities that cannot be categorized as being a single class of molecule, and thus SBO_FUNCTIONAL_ENTITY is the appropriate type.


text:
datasets or biobanks with both (WGS or genotype intensity) and methylation data
My aim is to associate CNVs to DNAm, but I don't know where to find datasets which contain both (sequencing + methylation). Do you know any?


text:
consortium, biobank or browser with CNV eQTLs
I found that GTEx contains several SNPs eQTLs, but I'm looking for previously reported CNV-eQTLs, where can I find such data (freely available)?


text:
Can multiplexing in Sequel II SMRTcells reduce the coverage?
I would like to have high depth of coverage of the transcriptome, but I need to analyze several tissues. I was suggested to pool all 5 tissue samples in same SMRT 8M cell. Will this result in a smaller coverage when compared to non-multiplex version?

Well, necessarily multiplexing reduces coverage of each sample. 1M reads / 1 sample > 1M reads / 5 samples (assuming even coverage). Is this what you are asking?

Yes, all right makes sense. Just another question. Are these reads per base pair or per insert?

A little unsure what that's asking, but if I understand you to mean library inserts, it will be reads per insert (species of DNA in the library). That should be strongly correlated to genome/transcriptome bases covered, so unless you have some rather wacky insert sizes it will imply reads per base as well.


text:
How does one distinguish nuclear DNA from mitochondrial DNA when doing WGS?
I'm interested in doing de-novo sequencing but also phylogenetic analysis. In particular, after de-novo sequencing and annotating the genome, I need to align the CO1 gene and the nuclear 28S rRNA gene of several species. When extracting DNA and sequencing with PacBio for example, how does the assembler know what corresponds to nuclear and what to mitochondrial DNA?

Does it need to? I mean, and I stress that I have never done de-novo sequencing so I may be completely wrong here, but won't the assembler simply try to assemble the longest contiguous contigs it can? And, if so, won't that naturally result in the mitochondrial sequences getting assembled together in one set of contigs and the non-mito sequences in others? I would expect this to be even easier with long read technology like PacBio. Finally, if you're only interested in a couple of genes, can't you directly target the highly conserved parts of them to amplify and sequence?

Cannot comment on CO1, but rrRNA should be different enough between the mitochondria and the nucleus as to not be confused. For one, mitochondria (like bacteria) should have 23S subunits instead of 28S. Many programs (e.g. barrnap, whatever is used by Prokka) automatically distinguish them.

Could you kindly review the answer provided? Please do consider upvoting or accepting an answer.

According to this pdf,  from pacb.com
isolation with mitoISO MITOISO
The kit will include three steps: 1) mechanical
rupturing of tissues or cells 2) removing cellular debris and nuclei by low speed centrifugation and 3) harvesting
mitochondria by high speed centrifugation.
Then, you can follow the DNA extraction with mitochondrial DNA, and then sequenced on the PacBio Sequel.
If you are doing whole genome sequencing, the mitogenome will have a higher read depth according to this article stating that lower read depth of the nuclear genome separating the two different origins:

those of mitochondrial origin & those of nuclear origin


text:
How can Iso-Seq reverse transcriptase artifacts be avoided?
My end goal is annotate a de-novo assembled genome. When trying to select the best method for transcriptome assembly I read Iso-seq was the preferred method. However other people suggested Nanopore as reverse transcriptase artifacts can be avoided. Is there any other way one can avoid these artifacts without the need for Nanopore?
Also, what is the recommended depth? Is using more SMRTcells worth the cost?

It is possible to filter out iso-seq reads with RT artifacts, though I don't know how effective these informatics approaches are. With nanopore, you will get more false junctions in alignment or transcript assembly due to sequencing errors. I more often see iso-seq used for new species.

Do you think polishing with RNA-seq reads is worth it? Or is Iso-Seq alone good enough?

@Caterina My other answer is also applied here: polishing iso-seq data with short reads is likely to reduce the accuracy.


text:
WGCNA hub genes error
I am trying to find hub genes for modules made using WGCNA:
hub = chooseTopHubInEachModule(datExpr, moduleColors)

But I am getting error:
Error in hubs[m] <- colnames(adj)[hub] : replacement has length zero

Thank you for the help1

Cross-posted on biostars: https://www.biostars.org/p/9486378/ Please don't shotgun post on multiple fora. See this discussion: https://bioinformatics.meta.stackexchange.com/questions/265/about-posting-queries-in-stack-bioinformatics

Can you please share the code, how you fixed this problem? I am facing the same problem and my data has colnames.
Any help will be appreciated.
Thanks

The most likely culprit here is that your expression data (the datExpr input to chooseTopHubInEachModule) has no colnames. Check that your data do have appropriate column names. Also check that the module color vector you supply has at least actual module (type table(moduleColors) and check that you have at least one module apart from "grey" or 0).

thanks as it worked. I have another question for extracting genes within modules of interest. I used unsigned network and interested in two modules for a trait. Earlier I made scatter plot and extracted the genes in all the modules. But now I only want to extract genes in these two modules for further analysis. I am not able to find a tutorial with more than one module. Can you please suggest something? Thank you1

Well, there are a few sections of the tutorial that may help. You can export the entire network analysis results (which you should do anyway) as illustrated in section 3d and then restrict the table (or just the gene identifiers) to the two modules you want. Alternatively, section 4a shows how to export gene identifiers in selected modules.


text:
no expression for hub genes identified in wgcna
I tried to identify intramodular hub genes in WGCNA analysis. I used:
hub = chooseTopHubInEachModule(datExpr, moduleColors)

But when I checked the hub genes for the modules, I found that the hub genes are expressed only in less than 5 out of 300 samples. I did not filter the gene expression data for any low expressed genes in certain number of samples. How does this make biological relevance if the hub genes are not expressed in most of the samples?
Thank you!

It could be that your modules group together genes (over)expressed in only certain small groups of samples that may be biologically relevant. You would need to check whether the samples in which the hubs (and presumably the modules) are up share some sort of a characteristic.
In general though, there are good reasons why authors of WGCNA recommend filtering out non- and low-expressed genes - correlations tend to not be a good co-expression measure when counts are very low, and technical effects can be much more pronounced. It is also recommended to spend a bit of time on preprocessing to make sure the data set does not have strong global drivers that are technical or uninteresting.

thank you! I applied some filters on RNA-seq data input, gsg = goodSamplesGenes(GE.adjusted,  minFraction = 1/2, minNSamples = 150, verbose = 3). I still have the genes which only have expression values for even one sample out of 300. I have a lot of zero's for no gene expression in the input file. Does WGCNA treat zero's as no or low gene expression? Should I replace zero's with NA so that it drop off the genes that have gene expression for less than say 100 samples? Thank you!


text:
GROMACS RMSF Result Analysis Problem
This is RMSF for protein from GROMACS after 100ns simulation-

I have few questions-

Residue numbers are starting from 1650 to 1850, has it ignore LIG
and HOH as residues?
There are three lines, so are there three chains and their residue number is same?
Are the straight lines connecting these chains or this is kind of loop something?
Since I am seeing 2nm or 20A fluctuation, is the protein quite disturbed?

The residue numbering in a PDB file mainly depend on the person who deposited it. There are usually reasons for how things are numbered, but ultimately you should look at PDB file, and check the details on how the structure was obtained. To give you an example, during my PhD I worked on a protein that had some missing residues at the N-terminus. That was because the N-term was very unstable and they didn't manage to crystallise it properly.
You need to check the PDB, but it is possible that the 3 chains have the same residue numbering (A1-200, B1-200). This is common in structures where you have multiple identical sub-units (e.g. dimers).
The lines are merely an artefact. Those lines are the attempt of the engine to draw a single line for all the chains. You go from residue 1850 of chain 1 to residue 1650 of chain 2, and so on. Some tools handle this behaviour, others don't. I would suggest you to parse the PDB yourself and plot (e.g. ggplot) each chain as a layer of the same figures .
It is indeed, but only you can tell if this is good or bad. It really depends on how you calculate the RMSF. Was it on the backbone or on the side-chains? How is your dynamics? High temperature? Does you protein partially unfold? are the chains separating?

Thanks for your response and clarification to understand the result.


text:
BWA-MEM2 alignment-Snakemake
I have started using snakemake 6.5.2 to align fastq files with reference file. I have pasted the error below in this question. How to allocate memory in the snakefile and read the header from samfile, '-'. This is the snakefile (wrapper for running alignment):
 rule bwa_mem2_mem:
        input:
            reads=["/scicore/home/cichon/GROUP/test_workflow/samples/{sample}.1.fq", "/scicore/home/cichon/GROUP/test_workflow/samples/{sample}.2.fq"]
        output:
            "/scicore/home/cichon/GROUP/test_workflow/output/{sample}.bam"
        log:
            "/scicore/home/cichon/GROUP/test_workflow/output/{sample}.log"
        params:
            index="/scicore/home/cichon/GROUP/test_workflow/data/gch38.fa",
            extra=r"-R '@RG\tID:{sample}\tSM:{sample}'",
            sort="none",             # Can be 'none', 'samtools' or 'picard'.
            sort_order="coordinate", # Can be 'coordinate' (default) or 'queryname'.
            sort_extra=""

        threads: 8
        wrapper:
            "0.77.0/bio/bwa-mem2/mem"

I used this command to execute the above snakemake file:
snakemake --use-conda --snakefile Snakefile --cores -p /scicore/home/cichon/GROUP/test_workflow/output/DE98NGSUKBD117612_1.bam /scicore/home/cichon/GROUP/test_workflow/output/DE71NGSUKBD117613_3.bam

Error from log file:
Looking to launch executable "/scicore/home/cichon/GROUP/test_workflow/.snakemake/conda/0ecf6ba401c70d7501c73bb1c733a060/bin/bwa-mem2.avx", simd = .avx
Launching executable "/scicore/home/cichon/GROUP/test_workflow/.snakemake/conda/0ecf6ba401c70d7501c73bb1c733a060/bin/bwa-mem2.avx"
-----------------------------
Executing in AVX mode!!
-----------------------------
* SA compression enabled with xfactor: 8
* Ref file: /scicore/home/cichon/GROUP/test_workflow/data/gch38.fa
* Entering FMI_search
* Index file found. Loading index from /scicore/home/cichon/GROUP/test_workflow/data/gch38.fa.bwt.2bit.64
* Reference seq len for bi-index = 6544178411
* sentinel-index: 2782712944
* Count:
0,  1
1,  1914743810
2,  3272089206
3,  4629434602
4,  6544178411

* Reading other elements of the index from files /scicore/home/cichon/GROUP/test_workflow/data/gch38.fa
* Index prefix: /scicore/home/cichon/GROUP/test_workflow/data/gch38.fa
* Read 0 ALT contigs
[fread] Cannot allocate memory
[main_samview] fail to read the header from "-".

Thanks
Priya

Please don't close the question. I am the one posted the question in Biostars as I didn't get a clear answer, I have posted it here again . Thanks

Clearly, this seems a problem with memory allocation. You should tell us what you tried from the suggestions of Biostars so we can help you in your problem.


text:
Dictionary and index of vcf for base recalibration step
I need suggestions in creating index and dictionary with vcf files. For the base recalibration step, I downloaded Homo_sapiens_assembly38.known_indels.vcf.gz from the given link: https://console.cloud.google.com/storage/browser/genomics-public-data/resources/broad/hg38/v0;tab=objects?prefix=&forceOnObjectsSortingFiltering=false
I am using GATK version 4.2.2.0
When I try to create index and dictoionary of the vcf file, I get the following error: no suitable codecs found. I also tried to check the vcf file using ValidateVariants and I get the same error: no suitable codecs found.
Can anyone suggest which indels vcf file is used for masking sites in the base recalibration step. Has anyone come across this above error and how did they resolve the issue.
Thanks

Can you show which commands are you trying on the file?

I have also encountered this error before in cases where there was something wrong with the VCF file I was supplying.

Double check your command to ensure that you are feeding it the file you believed you were and verify that file is a vcf (gzipped or otherwise)
If you know how, please compare MD5 or SHA checksums of your copy vs. the original to ensure you have an uncorrupted copy.
If this is something you haven't dealt with before, a less effective version would be to check and make sure that the file size on the server is the same as the file size on your machine.
Additionally, you may want to just use the less command to look inside the file and see if anything appears to be corrupted (you may need to decompress it for that purpose).
It may also be worth trying to gunzip or gzip -d commands to decompress the file (assuming you have the drive space).  This will both ensure that decompressing the file on the fly is not causing any grief to GATK (which is unlikely) and that the file decompresses without error (which would also be an indication of corruption).
Lastly, you could try downloading a fresh version.
Please let me know if this helps and you discover the source of the error.


text:
Proper formatting of .ped files, in order to run ADMIXTURE
I'm trying to run ADMIXTURE on some files I've generated using the Stacks pipeline. ADMIXTURE accepts .bed and .ped files, but when I try to pass the .ped file generated by Stacks to ADMIXTURE (v1.3.0), I get a "PLINK Input file error". I've pasted the head of my (tab separated) .ped file below.

head populations.plink.ped -n 1
garden  QUAC_G_SH_Q1485 0       0       0       0       0       0       T       T       C       C       G       G       T       T       C       C       C       C       G       G       T

Alternatively, I could convert my .ped file to .bed using plink, and try to pass this to ADMIXTURE--but Stacks didn't provide me with a map file to do this, so I don't think that will work.
If anyone knows what needs to change with my .ped file (remove sample name? population designation?), I'd appreciate any tips. This is my first time running most of these programs.
Thanks!

Did you obtain a solution to this? I am currently facing the same issue. Thanks

@AbosedeOlarewaju this would be a good question to ask, and include some specifics of your project.

Admixture requires a '12'-coded .ped file, rather than the allele-coding you're using. 'plink --file [file prefix] --recode 12' produced a .ped file that Admixture could use for me.


text:
Approaches for filtering bacterial/fungal contaminant sequences from RADseq results
I'm working with RADseq data (288 compressed FASTQ files, *.fq.gz) from 3 plates of plant tissue. I've done demultiplexing and preliminary QC (using the Stacks process_radtags command), and before performing a reference alignment, I want make sure the sequences in each sample (roughly) correspond to the organism I expect them to be derived from, rather than potential contaminants (bacteria and/or fungi).
Does anyone have a straightforward, "canned" approach for doing this, either in R or a different environment? I believe the blastn algorithm can be used to query the sequences in these files against a database (like NCBI), and is usable on the command line. However, I've not used it before, and I'm not sure if it makes sense to query all the sequences in each file. Additionally, I don't know if the command can process this many files of this size (~600 MB per FASTQ file).
This post from ~4 years ago asks a similar question for a pair of human blood samples, but is concerned with bacterial/viral DNA. The answer linking kmerfinder might have promise, but I'd appreciate any input from folks familiar with this issue.
(Background: I've gotten my Master's degree in Biology and have worked with NGS data for ~4 years now. I'm performing analyses on a 32-core Linux server. I'm less skilled in Python, but am willing to use it if there are no other options.)

Do you have a high quality reference genome of the plant you are working with?

Yes, I do. I've used the same reference genome in other projects, with similarly related taxa.

Thanks for the question, I'm familiar with working on RAD format and particularly don't know whether you've generated mixed genome RAD data. If know this information and I've time I'll put an answer down using either Bash or Python. I don't use R unless I've no other option. Please post a note below the question if you do edit it to notify me.

I ultimately addressed this issue using the kraken2 software. I found this approach to rely on a good deal of my server's memory, but it was sufficient for doing a viral/bacterial classification of sequences.


text:
Confidence intervals around Q-matrix values from multiple STRUCTURE runs
Context: I've run STRUCTURE on a set of 84 individuals (RADseq data) for a range of K values (2--7). I've specified 20 replicate STRUCTURE runs per K value, and these were passed to CLUMPAK to summarize results within and across K values, generating a summary Q-matrix for each K value.
Question: is there a way to calculate confidence intervals surrounding these summary Q-matrices? My preference is to do so using CLUMPAK, but I am open to other options if they're available. The output of each STRUCTURE run generates a text file containing a Q-matrix, but it isn't clear to me how to combine these for multiple runs to generate confidence intervals.
Thanks for any suggestions!

If we consult the structure docs, we observe that they directly address the idea of a $Q$ distribution for each individual (p. 24), in the form of a Bayesian credible interval, which can be obtained by using the ANCESTDIST flag. Example output is shown on p.34.
It is not clear to me how to summarize across runs given the different initializations, but I'd trust the internal methods more.
I am unsure how to integrate this with CLUMPAK. I'd say that they are probably two different ways of getting at a conceptually similar interpretation.

Thanks--setting the ANCESTDIST flag to 1 did generate Bayesian Credible intervals surrounding each value in the Q matrix, for each STRUICTURE run.
I will look into how to combine these results across runs.


text:
Protein-Protein Simulation using Gromacs
Can we do protein-protein simulation through Gromacs? if yes, can anyone teach me the correct procedure to get this done.
Thank you,

This is a very broad question. Could you kindly limit it to a specific technical issue? Secondly a reproducible example is always helpful. FYI the question has received a closed request.

I'd suggest looking at one of the many tutorials on the subject


text:
running fragpipe
I am getting the following error when I run FragPipe on my Ubuntu machine 22.04 how do I fix this?

Thanks for your question. As a general point it's always good to write out the error in asci rather than have an image. It makes it difficult to read and remember at least on my browser. The context of the question is processing mass-spec and I don't know if anyone here knows (maybe someone does). Is there a file format issue per chance?

the file format that I use as input is control.raw and treated.raw for fragpipe and affter running it i got the following log which is below:

Checking database...
Failed in checking /home/ms/MS/Control_1.raw. Will ignore it.
Failed in checking /home/ms/MS/Treated_1.raw. Will ignore it.
Batmass-IO binaries for Thermo support and/or Thermo native libraries not found
Batmass-IO binaries for Thermo support and/or Thermo native libraries not found
There are corrupted files. Please remove those files and re-start the task.
Checking spectral files...
/home/ms/MS/Control_1.raw: Scans = 0
/home/ms/MS/Treated_1.raw: Scans = 0
Process 'MSFragger' finished, exit code: 1
Process returned non-zero exit code, stopping

Edit your post and add the plain text content there with proper formatting.

Thanks for your valuable time, well my error has been resolved through Fragpipe support Here I am adding the discussion link below:
https://github.com/Nesvilab/FragPipe/issues/1207


text:
Accessing expression data in an ExpressionSet
I used rma from the oligo package to normalize CEL files, returning an object of class ExpressionSet. How do I access the actual expression data?
> data = read.celfiles(names)
> data2 = oligo::rma(data)
class(data2)
[1] "ExpressionSet"
attr(,"package")
[1] "Biobase"
> data3 = as.matrix(data2)
Error in as.vector(data) :
  no method for coercing this S4 class to a vector

Have you looked at the oligo::rma documentation? Can you clarify why this doesn't make sense to you?


text:
java error when running CIRI-full.jar Merge module
I have an issue running a module trying to detect circRNA from RNAseq data. The error is a java error and I have not knowledge of java so I don't understand the issues and googling hasn't helped.
the issue arises when I try to run this module of the algorithm
https://ciri-cookbook.readthedocs.io/en/latest/CIRI-full.html#the-merge-module
as I try to run the module I get the following output and error
Merge module start

Loading Annotation...
Annotation Loaded
Loading CIRI_AS output...
CIRI_AS output Loaded
Loading CIRI output...
CIRI output Loaded
Loading CIRI_RO output...
CIRI_RO output Loaded
Combine AS and RO output
Combine completed, 71 reads are used
Exception in thread "main" java.lang.reflect.InvocationTargetException
        at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
        at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:62)
        at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
        at java.lang.reflect.Method.invoke(Method.java:483)
        at org.eclipse.jdt.internal.jarinjarloader.JarRsrcLoader.main(JarRsrcLoader.java:58)
Caused by: java.lang.NullPointerException
        at Exon_gtf.get_exon(Exon_gtf.java:16)
        at Merge3.getmerge(Merge3.java:784)
        at CIRI_Full2.main(CIRI_Full2.java:551)
        ... 5 more

All prior steps preceding this one seem to have run properly and I was able to generate all the files to run the Merge module.
any help or explanation as to what the issue could be would be greatly appreciated. As I mentioned before I have no experience in java.
kind regards

Whatever the reason, it shouldn't crash with a null pointer exception.  You should bring that to the attention of the developers so they can try to debug.  Does it run successfully with their test dataset?  That might help narrow it down.

Hi Jesse, I did bring it up on github but no response so far. I have run the same pipeline of another dataset of mine and it worked. So the issue might be data I am feeding it. Thanks

Makes sense (something about the input data crashing the software).  It might be meaningful that the crash is within a GTF-handling function, too.  Are you giving it the optional GTF annotation file, and if so, does that same file work with the other dataset?  If you can make a small enough reproducible test case to post here that shows the crash, maybe someone can help narrow it down.

thanks for the feedback. I don't use the optional GTF. Comparing the inputs for that piece of code there is a clear difference in size so that means the difference is in there. I discussed with our collaborator who generated the data and they did process it differently compared to the previous one that did work. So the problem is not the code but the input and how the data was generated.


text:
samtools view command not found error
When I tried to use samtools to split a bam file based on different chromosomes, I used this command:
samtools view input.bam -b chr21 | chr21.bam

However, I get error messages like this:
-bash: chr21.bam: command not found
[W::hts_idx_load3] The index file is older than the data file: input.bam.bai

How could I solve this problem? Thanks in advance~~~
The version I am using is:
Program: samtools (Tools for alignments in the SAM format)
Version: 1.10 (using htslib 1.10.2)

Try re-indexing your bam file:
samtools index input.bam

Or update the timestamp via touch index.bam.bai but I think re-indexing is safer. You might want to save any original index file first, e.g. via mv input.bam.bai old_input.bam.bai or similar
Also I think you meant > chr21.bam instead of piping via | chr21.bam if you want to save the output file

Cool! I think the problem is because I use | chr21.bam.


text:
How to see effect of point mutation in PyMOL?
I am trying to model how a specific point mutation would affect the protein structure.  So far, I've figured out how to create the point mutation in PyMOL using the mutagenesis wizard.  I'm not sure how to take the next step toward actually visualizing the consequences of the mutation though (i.e. would it disrupt a helix, etc).  Is there a way to do this in PyMOL, or would I need other software?  Thank you!

I would have said the question lacked detail, but the answers are good. Duet and equivalents are good. Pymol creates mutations, but can't analyse them.

PyMOL just find the rotamer that causes the least clashes, it does not repack other nearby residues, alter the backbone or adopt a conformer that is not part of the rotamer set it knows, this is a problem in some complicated cases core cases, but in in most cases the clashes reflect what the ∆∆G may be so I have to reluctantly agree here...

A common metric used to assess if a mutation destabilises the structure is the difference in Gibbs free energy. PyMOL does not give a ∆∆G (difference in Gibbs free energy, kcal/mol), but a very crude steric clash nor does it repack neighbouring sidechains out of the way, or alter the backbone, so unfortunately is not the way forward. However, if one must there is a way (see SIMBAI below).
Beyond destabilisation
Destabilisation predictors have a ~60% or lower success at predicting if something is pathogenic (ref). This is because the effect a mutation has on a protein is rather complicated as you can have destabilisation, but also aggregation, loss of regulation, mislocalisation, loss of binding to another protein, loss of binding to ligand, loss of activity (catalysis, inability to change conformation) etc etc. Then you can get curious subcases, such as loss of activity that result in the sequestration of a binding partner or homologues in a complex resulting in a worse phenotype than destabilisation or truncation.
Modelling
The behaviour of a protein can be assess by the energy of the system using forcefields, which are multi-term equations that describes the atomic interactions (attraction, repulsion etc.). The result is the Gibbs free energy mentioned above, which is a potential, so negative is good. Kinetic energy is release when the unfolded protein "rolls down" down the potential energy funnel —ProTherm dataset (an empirical dataset of ∆∆Gs) and the tool SDM give this value inverted as the negative can be confusing.
An uphill change in the Gibbs free energy (∆∆G > 0) generally correlates with a destabilisation of a protein —aggregation and solubility being the caveat in the "generally". The destabilisation a single state (static conformation) can be assessed. gmx energy in Gromacs, FoldX or Rosetta can be used for this.
For altered ability to switch conformation steered MD simulations are required, but generally the ∆∆G at different states can be used.
For the non-destabilisation options, https://www.phosphosite.org/ and http://elm.eu.org/ are very useful.
I have been working on a server that gives a ∆∆G value and note if there are any nearby post-translational modifications, gnomAD or changes in potential linear motifs etc.: Venus (unpublished). However, there are many ∆∆G calculation servers out there and are compared in many papers (example), but it should be said that the dataset used in these comparisons is always ultra-biased despite how much is tried. Missense by the Phyre2 group (Michael Sternberg) give detailed information of the neighbourhood changes the mutation causes.
In terms of extra info, Miscast gives Uniprot info.
SIMBA-I via PyMOL
In terms of PyMOL, the SIMBA-I algorithm can be used by getting the RSA from it —I would not recommend it as it may underestimate the severity of unusual cases. But if one must, in PyMOL one can paste the following into the console with the wild type variant:
# defined by you so change accordingly
from_residue = '' # single letter
to_residue = '' # single letter
resi =
chain = ''

## data
maxASA = {'A': 121.0, 'R': 265.0, 'N': 187.0, 'D': 187.0, 'C': 148.0, 'E': 214.0, 'Q': 214.0, 'G': 97.0, 'H': 216.0,\
          'I': 195.0, 'L': 191.0, 'K': 230.0, 'M': 203.0, 'F': 228.0, 'P': 154.0, 'S': 143.0, 'T': 163.0,\
              'W': 264.0, 'Y': 255.0, 'V': 165.0}
volumes = {'A': 88.6, 'R': 173.4, 'N': 114.1, 'D': 111.1, 'C': 108.5, 'Q': 143.8, 'E': 138.4, 'G': 60.1,\
               'H': 153.2, 'I': 166.7, 'L': 166.7, 'K': 168.6, 'M': 162.9, 'F': 189.9, 'P': 112.7, 'S': 89.0,\
               'T': 116.1, 'W': 227.8, 'Y': 193.6, 'V': 140.0}
hydrophobicities = {'A': 8.1, 'R': 10.5, 'N': 11.6, 'D': 13.0, 'C': 5.5, 'Q': 10.5, 'E': 12.3, 'G': 9.0,\
                    'H': 10.4, 'I': 5.2, 'L': 4.9, 'K': 11.3, 'M': 5.7, 'F': 5.2, 'P': 8.0, 'S': 9.2, 'T': 8.6,\
                    'W': 5.4, 'Y': 6.2, 'V': 5.9}
cmd.set('dot_solvent', 1)
cmd.set('dot_density', 3)
sasa = cmd.get_area(f'resi {resi} and chain {chain}')
rsa = sasa / maxASA[from_residue]
offset = -1.6
rsa_factor = +2.2
dV_factor = +0.59
dH_factor = -0.29
rsa_dV_factor = -0.54
rsa_dH_factor = +0.49
# calculate
dV = (volumes[to_residue] - volumes[from_residue]) / 100
dH = (hydrophobicities[to_residue] - hydrophobicities[from_residue])
simba = offset + rsa_factor * rsa + dV_factor * dV + dH_factor * dH + rsa_dV_factor * rsa * dV + rsa_dH_factor * rsa * dH
print -simba  # negative stabilised, positive destabilised

Human reasoning
However, the major problem with all the tools mentioned is that human reasoning is often required to read the literature and make a solid hypothesis.
A resource I cannot stress enough, for human mutations, its utility is gnomAD, which contains phenotypically neutral mutations: for example, one really ought to not claim a homozygous missense is pathogenic because it is causing a loss of function by destabilising the protein, yet there are homozygous truncations roaming around in the human population.
The paper describing the crystal structure used and those of the homologues may point out key features which would missed by automated servers, such as secondary-structure elements that are speculated to perform certain roles.
In fact, without doing any calculations seeing in PyMOL if a mutation causes a clash in the core of a protein, alters the active site, or changes a residue to a smaller one in a region known to alter conformation upon a given condition, is a lot more informative even if not quantitative...

You could run a molecular dynamics simulation for example with Gromacs or Charm. The problem with pymol is that the structure is static and you don't know how much tension that point mutation creates. First of all you would run an energy minimization to remove bad non-physical overlaps between atoms, then add solvent molecules, and then simulate your system first with position restraints for equilibration and then without for 1ns to 100ns+.
Even then you cannot be sure that what comes out is the 'final' result as sometimes an induced fold change can take much longer. It depends on the system and the energy barriers included and the time scales that you are interested in. If the protein is large it may have many different physiologically relevant conformations and the point mutation might only change the equilibrium distribution of those states or totally prevent a transition.
In case you have more than one crystal structure (in different conformations) you may introduce these point mutations in all of them and simulate them for a while. But it is not that trivial to setup such a simulation you should contact an expert. There are many things that can go wrong. Other software you could use are Modeller and Rosetta, if I am not mistaken.

You could give SAAPpred a try.
This tool uses machine learning (in the form of random forest) to predict the local effects of a single point mutation on the protein structure.
I did my PhD in the group and I used SAAPpred to find the most (predicted) damaging mutations of G6PD, to then ran extensive MD simulations to evaluate the final stability of those mutants.
Keep in mind that the tools doesn't return a structure, but just the possible/likely effects of that substitution.

Cool this looks good.

then ran extensive MD simulations to evaluate the final stability of those mutants..... Did You try to mutare the actual protein too ? Any interesting validation result ?


text:
GTF file format causing GeneMark-ES error
I am predicting genes using GeneMark-ES, I managed to generate a gtf file, but when I tried to extract the nt/protein sequence using get_sequence_from_GTF.pl, I am running into an error in the format of the gtf file, which I find puzzing as it was the gtf file was generated using gmes_petap.pl (an script in the GeneMark-ES package).
Here is the command I used and the error message which also prints the first line of the genemark.gtf file:
$ get_sequence_from_GTF.pl genemark.gtf DT_CAJQTU010000037_1.fasta

error in file format GTF: CAJQTU010000037.1 Desoria tigrina genome assembly, contig: utg000040l, whole genome shotgun sequence     GeneMark.hmm3   gene    2099    4039    .       +       .       gene_id "1_g"; transcript_id "1_t";

I think I have found the line in get_sequence_from_GTF.pl that is causing the error, but not being a perl monk, I am having troubles to read figure out what is wrong with my gtf file and more importantly figure out how to fix it.
...
if ( /^(\S+)\t\S+\t(\S+)\t(\d+)\t(\d+)\t.\t([\+-])\t([012.])\tgene_id \"(\S+)\"/ )
    {
    ... (some stuff)
    }
else
{
  if ( /^\s+$/ ) { next };
  if ( /^#/ ) { next };
  die( "error in file format GTF: $_\n" );
}
...

I'm afraid this is a clear bug in the tool you are using. The regular expression requires that all fields in the GTF file contain no whitespace. The \S+ means "one or more non-whitespace characters" but your "source" field, the second field in the GTF file, contains spaces: Desoria tigrina genome assembly, contig: utg000040l, whole genome shotgun sequence.
The thing is, that is a perfectly valid field for a GTF file! GTF fields can contain spaces, they just can't contain tabs. I tried to find a repository for the software to suggest a fix, but found only an empty one. In any case, try changing the regex you show to this:
#if ( /^(\S+)\t\S+\t(\S+)\t(\d+)\t(\d+)\t.\t([\+-])\t([012.])\tgene_id \"(\S+)\"/ )
if ( /^([^\t]+)\t[^\t]+\t([^\t]+)\t(\d+)\t(\d+)\t.\t([\+-])\t([012.])\tgene_id \"(\S+)\"/ )
    {
    ... (some stuff)
    }

I can't guarantee that this will work as expected, please test! The details will depend on what they do with the captured groups (the some stuff), but at least that regex will match and I think I have kept the numbering of the capture groups.

I agree with this 100%. I read the regex and didn't spot the reply. There are quite a few problems, the use of a bare "." just means 'any character' rather '.' which in this example would fit (albeit it might be intentional). The final section '; transcript_id "1_t";' isn't captured. However, the key problem is \S which will trip over the space characters in "CAJQTU010000037.1 Desoria tigrina genome assembly, contig: utg000040l, whole genome shotgun sequence". The use of '[^\t]', i.e. not a tab character will work. Personally I would split on the tab character.

Just to clarify this "\t.\t" means any single character between two tabs. However, this "[012.]" might mean a \. i.e. just a full-stop, I can't remember. If it doesn't [.] is exactly the same as [012.] or even [a-z|0-9] and beyond. Anyway this is the accepted answer.

You might use agat_sp_extract_sequences.pl from AGAT instead.
To extract and translate the coding regions:
agat_sp_extract_sequences.pl -g infile.gxf -f infile.fasta -p
see here for all possibilities with that script: https://agat.readthedocs.io/en/latest/tools/agat_sp_extract_sequences.html
Otherwise with the modification suggested by @terdon, your script should work fine.


text:
Help with MinION sequencing data species identification
Hi I'm new to bioinformatics and have just completed my first run on the MinION (long read sequencing Oxford Nanopore Technologies). I was hoping someone could direct me towards R packages, workflow, tutorials or guides that will help me identify species that are present in my sample mainly for fungi as my amplicons had ITS1 and ITS2 regions.
I have been watching Youtube tutorials on RStudio for over a week and I feel like I'm no closer than I was when I started. I have both fastq and fast5 files and I tried to run the fastq files using EPI2ME, however, the database is lacking in fungi. I have tried using NCBI's BLAST but I can only run <5kb reads and most of my sequencing reads are >6kb. Also, I've tried using Galaxy to convert my fastq files to fasta but I'm not sure how to run my long reads for species identification using Galaxy or if that is even possible.
Any help would be greatly appreciated. Thanks.

What sequencing kit did you use? You mention amplicons... are you amplifying specific regions, or doing a whole-genome metagenomic analysis?

There are a number of useful metagenomics tools out there, including Centrifuge and Kraken2, as @gringer has mentioned.
Also, FWIW, the WIMP tool from ONT, a Centrifuge based metageonomic classifier, does include fungus sequences derived from RefSeq, though it may not include the species you are looking for.
If you want to use NCBI BLAST, you could split your long reads into pieces that are <5kb (there are a number of tools out there for FASTQ manipulation - my current favorite tool for this is seqkit, where you could do something like this):
seqkit sliding -s 4900 -W 4900 test1.fastq
The above will split your FASTQ reads into chunks of 4.9kb. You can adjust the -s parameter if you instead want overlapping chunks.
From a GUI standpoint, you can use Galaxy for the above by making a workflow with the following steps:

Trim Sequences (enter reasonable numbers here, to trim each of your reads to <5kb)
FASTQ to FASTA converter
NCBI BLAST+

Another possibility is QIIME2, which can be used with the UNITE fungal training set. See recent reference (looking at the ITS1/ITS2 regions), UNITE, and QIIME2.
Long term, you may want to look into using the command line, as it is much more powerful and customizable, but if BLAST gives you reasonable results for your reads, the above process in Galaxy may work for you.

I've found Centrifuge and Kraken2 to be useful for species classification. There are pre-prepared databases that include Fungi, e.g. see here for the Kraken 2 "PF" database, which is the Standard database plus protozoa & fungi.
Can't really provide much more help than that with more information about what you've done.
You'll need a computer with around 64 GB of memory in order to be able to run these programs, and it's probably not the best idea to jump right into command line compilations when starting out... but it looks like you're up for a challenge given that your first MinION run has involved jumping into the deep end with a fungi metagenomic project.


text:
ncbi fastq dump error in loop
I have a very basic objective which I want to give list of id to a ncbi-srafastq tool kit which will prefetch the id which is in .sra file extension and then it will run fastq-dump on the file to split the files into fastq files, single file if there is single end file or paired end file if its paired end or two set of files.
This is my code
#!/bin/bash

input=$1

a=$(basename "$input")
id=${a/.id/}

# Download SRA data
prefetch "$id"

# Extract SRA ID from the downloaded data
#sraID=$(ls -1 | grep -v run.sh)
sraID=$(ls -1 | grep -v run.sh | sed 's/\.sra$//')
# Validate the downloaded SRA data
vdb-validate "$sraID"

# Dump fastq data from the SRA file
fastq-dump "$sraID"  --gzip

# Concatenate fastq files for read1, read2, and unpaired reads (if available)
for s in "$sraID"
do
  echo "$s"
  if [ "$s" != "$id" ]
  then
    [ -f "${s}_1.fastq" ] && cat "${s}_1.fastq" >> "${id}_R1.fastq"
    [ -f "${s}_2.fastq" ] && cat "${s}_2.fastq" >> "${id}_R2.fastq"
    [ -f "${s}.fastq" ] && cat "${s}.fastq" >> "${id}.fastq"
  fi
done

# Compress fastq files using pigz
#for i in "${id}"*.fastq
#do
  #echo "$i"
 # pigz "$i" & sleep 1
#done
#wait

# Remove the original uncompressed fastq file if both read1 and read2 files are present
[ -f "${id}_R1.fastq.gz" ] && [ -f "${id}_R2.fastq.gz" ] && [ -f "${id}.fastq.gz" ] && rm "${id}.fastq.gz"

# Test deleting intermediate folder
echo "Space before tmp folder delete"
echo "Removing tmp folder $sraID "
ls *
rm -rf "$sraID"*
ls *
echo 'done'

I run using as small example file as such
./run.sh SRR12332899

The error I get is this
2023-08-05T12:37:32 prefetch.3.0.6: Current preference is set to retrieve SRA Normalized Format files with full base quality scores.
2023-08-05T12:37:33 prefetch.3.0.6: 1) Downloading 'SRR12332899'...
2023-08-05T12:37:33 prefetch.3.0.6: SRA Normalized Format file is being retrieved, if this is different from your preference, it may be due to current file availability.
2023-08-05T12:37:33 prefetch.3.0.6:  Downloading via HTTPS...
2023-08-05T12:40:48 prefetch.3.0.6:  HTTPS download succeed
2023-08-05T12:40:48 prefetch.3.0.6:  'SRR12332899' is valid
2023-08-05T12:40:48 prefetch.3.0.6: 1) 'SRR12332899' was downloaded successfully
2023-08-05T12:40:48 prefetch.3.0.6: 'SRR12332899' has 0 unresolved dependencies
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Database 'SRR12332899.sra' metadata: md5 ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Table 'SEQUENCE' metadata: md5 ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Column 'ALTREAD': checksums ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Column 'QUALITY': checksums ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Column 'READ': checksums ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Column 'SPOT_GROUP': checksums ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Database 'sra/SRR12332899.sra' contains only unaligned reads
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Database 'SRR12332899.sra' is consistent
2023-08-05T12:40:49 fastq-dump.3.0.6 err: item not found while constructing within virtual database module - the path 'sra' cannot be opened as database or table
fastq-dump quit with error code 3
sra
Space before tmp folder delete
Removing tmp folder sra
run.sh

sra:
SRR12332899.sra
run.sh
done

If I do just
prefetch SRR12332899 and then fastq-dump SRR12332899
I get a SRR12332899.fastq  SRR12332899.sra
This is for one file since I want to download list of files using the script which is not working.
Any suggestion would be really helpful

Note that you are catting an already gzipped files since you used --gzip in fastq-dump, so the suffix of the output must be .gz
It would be muchfaster to simply search for the dataset at sra-explorer.info and get fastq download links right away.

The main problem is that your script is inadvertently treating a set of file and directory names as a single SRA accession number, all in one variable.  You should update it to work with each accession separately and make sure they really are accessions.

On this line you're taking all the names of things in your working directory (minus the script), each on a separate line, and trimming .sra off the end of each line:
sraID=$(ls -1 | grep -v run.sh | sed 's/\.sra$//')

In your case that list was apparently just a directory named "sra".  (Even though I don't have that directory on my end, you can see it in your output, from the first of two ls * at the end, where you have sra/SRR12332899.sra.)  So later on it tries to fastq-dump sra, which gives you that "cannot be opened as database or table" error.  (vdb-validate can take a directory as input, so it doesn't complain.)  If I run your script from scratch in an empty directory with SRR12332899, it actually seems to run OK, aside from removing everything it creates at the end.  But if you run it with a second accession you see the problem, since it then tries things like vdb-validate "several accessions and filenames separated by newlines".
I think a proper fix would involve a bunch of changes, but here are some tips, at least:

use subdirectories to keep things organized so you don't have to juggle multiple types of file in one place
set -x at the start of a script will show you each command as it happens, which can be very helpful for troubleshooting
be careful with quotes and whitespace in bash.  If you want to loop over each word in a variable you can do for thing in $var; do ...; done (unquoted) but if you run into filenames with whitespace, it breaks them up.  Another pattern is command | while read thing; do ...; done which will go line-by-line.  (In that case if you were working with a variable storing the command output, you might do something like echo "$output" | while..., quoted that time, since you'd want to preserve the newlines.)
ls * will expand each file/directory name as separate arguments (so bash really ends up calling ls "sra" ...whatever else...), and when ls sees directory names is lists their contents too.  If you just want to show all the files/directories you can leave off the asterisk altogether, or if you wanted to give directory names but not have it list their contents, you can use ls -d.
you don't necessarily have to call ls to store the names as a block of text.  You can loop over each filename, including those with spaces, directly in bash (e.g. for thing in *; do echo "path: '$thing'"; done)
be careful with the rm command, especially combined with wildcards, and especially with -rf.  If sraID happened to be empty somehow, think of what rm -rf "$sraID"* would actually mean.  (And yes, scripts can remove themselves!)

will work on your suggestion and show you what i have done


text:
How to plot the gap distribution of contigs wrt to reference genome
I have a contigs file that I generated using Minia and then I have used Minimap2 to map these contigs again to the reference genome. Now I want to plot the gap distribution i.e. gap distance of the set of contigs with respect to the positions in the reference genome. For example: something like , if there are m contigs obtained from the Minia, then X axis will show 0 to m-1 gap points (in an ideal case where there is no overlap between the contigs and the contigs are aligned back to back without space).
The file from minimap2 looks like this:
0   17679   0   17679   +   NZ_BCML01000091.1   39426   18736   36415   17679   17679   60  NM:i:0  ms:i:35358  AS:i:35358  nn:i:0  tp:A:P  cm:i:3274   s1:i:17675  s2:i:49 de:f:0  rl:i:75 cg:Z:17679M
1   16011   0   16011   +   NZ_BCML01000021.1   105137  86      16097   16011   16011   60  NM:i:0  ms:i:32022  AS:i:32022  nn:i:0  tp:A:P  cm:i:3008   s1:i:15996  s2:i:0  de:f:0  rl:i:15 cg:Z:16011M
2   2779    0   2779    -   NZ_BCML01000048.1   72397   14902   17681   2779    2779    60  NM:i:0  ms:i:5558   AS:i:5558   nn:i:0  tp:A:P  cm:i:519    s1:i:2751   s2:i:0  de:f:0  rl:i:45 cg:Z:2779M
4   177     0   177     +   NZ_BCML01000150.1   16425   15710   15887   177     177     60  NM:i:0  ms:i:354    AS:i:354    nn:i:0  tp:A:P  cm:i:30     s1:i:160    s2:i:0  de:f:0  rl:i:15 cg:Z:177M

where the 1st column is the contig ID and 8th and 9th column show the start and end position of the contig on the reference genome.
I have used the samtools for obtaining a length distribution plot, but not understanding how to do a gap distribution plot. I can use awk on linux or python or R. Anything that is going to work. Thank you for all the help.

a better solution
I wasn't too happy with the previous answer that I gave, so I was fiddling around, and I realized that bedtools genomecov will give you this for free with the following workflow:
minimap2 -ax asm10 ref.fasta contigs.fasta | \
samtools sort --output_fmt BAM > contigs_to_ref.bam
# -bga means output in bedgraph format, including all positions
bedtools genomecov -ibam contigs_to_ref.bam -bga > contig_cov.bg

If you look at the resulting file, it will report all regions of identical coverage:
head contig_cov.bg
# prints the following:
# reference chr/contig, start position, end position, coverage
NODE_1_length_442924_cov_598.199    0   8151    0
NODE_1_length_442924_cov_598.199    8151    8392    1
NODE_1_length_442924_cov_598.199    8392    20203   0
NODE_1_length_442924_cov_598.199    20203   20463   1
NODE_1_length_442924_cov_598.199    20463   44953   0
NODE_1_length_442924_cov_598.199    44953   45216   1
NODE_1_length_442924_cov_598.199    45216   52633   0
NODE_1_length_442924_cov_598.199    52633   52899   1
NODE_1_length_442924_cov_598.199    52899   57698   0
NODE_1_length_442924_cov_598.199    57698   58207   1

You can simply filter this file for all regions wherein the fourth column is equal to zero, and take the first 3 columns, and that is a list of gaps in the assembly relative to the reference:
awk '$4 == "0" { print $0 }' contig_cov.bg | cut -f1-3
# prints
NODE_1_length_442924_cov_598.199    0   8151
NODE_1_length_442924_cov_598.199    8392    20203
NODE_1_length_442924_cov_598.199    20463   44953
NODE_1_length_442924_cov_598.199    45216   52633
NODE_1_length_442924_cov_598.199    52899   57698
NODE_1_length_442924_cov_598.199    58207   73398
NODE_1_length_442924_cov_598.199    73733   78601
NODE_1_length_442924_cov_598.199    78959   80147
NODE_1_length_442924_cov_598.199    80419   86967
NODE_1_length_442924_cov_598.199    87223   93872

You could even do this as a one-liner probably though that's getting a little convoluted.
Original answer
It's a bit circuitous, but one way to do this would be to follow this answer's suggestion and generate a pileup for the alignment of your contigs to the reference.
You can get minimap2 to give SAM/BAM output with the -a flag. You then use the SAM/BAM to generate a pileup or a coverage file, using samtools mpileup or samtools depth -a. You can then use a heuristic to find contiguous regions with zero coverage of contigs. You could for example loop through a depth file, here is some python code that gives the basic idea:
#!/usr/bin/env python

# minimap2 -ax asm10 ref.fasta contigs.fasta | \
# samtools sort --output_fmt BAM > contigs_to_ref.bam
# samtools depth contigs_to_ref.bam > depth.txt
depth_file = "depth.txt"

with open(depth_file) as dep_handle:
  gap_sizes = list()
  start = None
  pos = None
  for line in dep_handle:
    # needs also logic to handle multiple chromosomes in reference
    # but I'm lazy
    fields = line.split()
    if int(fields[2]) == 0:
      pos = int(fields[2])
      if start is None:
        start = int(fields[1])

    elif fields[2] >0:
      if start is None:
        continue
      else:
        gap_sizes.append(pos - start)
        start = None

print(gap_sizes)

Probably there is a simpler way to do this but I don't know it.


text:
Why do we count reverse compliments when counting kmers in a single DNA strand?
When looking for repeat patterns in a DNA sequence, we can look for the pattern and its reverse compliment with up to d mismatches. However, why do we look for reverse compliments if we're analyzing a single DNA strand?
I thought about it and my guess is that this lets us analyze both strands of DNA. For example, if we find the reverse compliment at position x, that means that the other strand of DNA has the pattern at position x - pattern.length.
Is that correct?

When we sequence genomes, the individual reads come from both strands, and that also applies for assembled sequences* (contigs / scaffolds). Therefore for most of the genome assemblies out there, strandeness of individual scaffolds is actually an arbitrary property, and reverse complementing any sequence will generate an equally valid representation (with reverted coordinates).
The coordinates of the reverse complementary k-mer matches depend on the notation. Usually, the notations are still in the coordinates of the reference, for example, blast gives for reverse complimentary strand matches the from subject coordinate greater than the to coordinate. Some others are denoting "the leftmost" position and strand encoded in a different column (e.g. sam).
*perhaps with the exception of chromosomal level assemblies

Do you have any suggested reading for insight on composing a single sequence from both strands? I feel like that's a big mental roadblock for me.

To be completely honest, I don't! I would be looking out for assembly literature (especially of short reads). I will also try to keep my eyes open and let you know if I run into something.


text:
Why must a maximal non-branching path be a contig?
The following is from Bioinformatics Algorithms:

Fortunately, we can derive contigs from the de Bruijn graph. A path in a graph is called non-branching if in(v) = out(v) = 1 for each intermediate node v of this path, i.e., for each node except possibly the starting and ending node of a path. A maximal non-branching path is a non-branching path that cannot be extended into a longer non-branching path. We are interested in these paths because the strings of nucleotides that they spell out must be present in any assembly with a given k-mer composition. For this reason, contigs correspond to strings spelled by maximal non-branching paths in the de Bruijn graph. For example, the de Bruijn graph below, constructed for the 3-mer composition of TAATGCCATGGGATGTT, has nine maximal non-branching paths that spell out the contigs TAAT, TGTT, TGCCAT, ATG, ATG, ATG, TGG, GGG, and GGAT. In practice, biologists have no choice but to break genomes into contigs, even in the case of perfect coverage (like in the figure below), since repeats prevent them from being able to infer a unique Eulerian path.

But this quote doesn't make sense to me:

We are interested in these paths because the strings of nucleotides that they spell out must be present in any assembly with a given k-mer composition.

Is it not possible that we simply have a missing read that would have otherwise branched this path? How do we know that a maximal non-branching path must be a contig?

The definition of contig is derived from assembly algorithms. A contig is not necessarily present in the genome. PS: i.e. "contig" is algorithmically defined, not biologically defined.

Is it not possible that we simply have a missing read that would have otherwise branched this path?

Yes, it is possible. Then you are likely to get a misassembled contig.

By the way, a maximal non-branching path in the assembly graph is more often called a unitig.


text:
Efficient way of mapping UniProt IDs to representative UniRef90 IDs?
I want to map ~500,000 UniProt IDs to their representative UniRef90 IDs.
What is the most efficient way of doing this?
Do we have a table with the mapping?

You can do this directly on UniProt: https://www.uniprot.org/uploadlists/
Just paste or upload your list of UniProt IDs, and select "UniProtKB AC/ID" in the "From" field and "UniParc" in the "To" field

I've also written a script, pasted below, that can do this with some useful options:
$ uniprot_map.pl -h

uniprot_map.pl will connect to UniProt and map IDs from the input file to ehatever ID was requested.

USAGE:
       uniprot_map.pl [options] <ID_LIST>

OPTIONS:
      -h : Print this help and exit.
      -f : From ID type. This is the type of ID found in the input file.
      -F : Format. See http://www.uniprot.org/faq/28#batch_retrieval_of_entries for details.

      -l : When a requested ID maps to many output IDs, print only the first (probably the best) match.
      -L : List all recognised ID types.
      -u : Print unreviewed (trEMBL) prots as well. By default, if an ID maps to SWISSPROT and trEMBL IDs
           only the SWISSPROTs will be printed, use this flag to print trEMBL as well.
      -t : To ID. This is the ID type that the input IDs will be mapped to.
      -v : Verbose. Print progress messages to STDERR.

If you save the script as uniprot_map.pl in a directory in your PATH and make it executable, then you can use a file like this:
$ cat file
P31946
P62258

And feed it to the script with the appropriate options to get the mapping:
$ uniprot_map.pl -f ID -t UniParc file
P31946  UPI000013C714
P62258  UPI0000021A46
2 of 2

The script:
#!/usr/bin/perl
use strict;
use warnings;
use LWP::UserAgent;
use Getopt::Std;
my %opts;
getopts('Lsvuhf:F:t:l', \%opts);
my $base = 'http://www.uniprot.org';
my $verbose=$opts{v}||undef;
my $tool = 'mapping';
my (@NAMES,@names);
my $print_only_first=$opts{l}||undef;
my $allow_trembl=$opts{u}||undef;
list_ids() if $opts{L};
&usage() if $opts{h};

&usage() unless $ARGV[0];

####################################
# Collect the IDs we want to match #
####################################
foreach my $ar(@ARGV){
if (-e $ar){
open(A, $ar);
while(<A>){
    chomp;
    push @names,$_;
}
}
else{push @names,$ar};
}

######################################
# What are we matching FROM and TO ? #
######################################
my $from=$opts{f}||'ID';
my $to = $opts{t}||'ACC';

if($opts{t}){
$to=translate_terms($to);
}
if($opts{f}){
    $from=translate_terms($from);
}
#################
# Output format #
#################
my $format = $opts{F}||'tab';
#my @query=

my $agent = LWP::UserAgent->new;
push @{$agent->requests_redirectable}, 'POST';

my $params= {
'from' => $from,
'to' => $to,
'format' => $format,
    'query' => "@names"
};

print STDERR "Requesting...\t" if $verbose;
my $try_again=1;
my $response;
while($try_again==1){
    $response = get_response($params);
    if(defined($response->is_success)){$try_again=0}
elsif($verbose){print STDERR "Failed, retrying...\n"; }
}

if( $print_only_first)
{
if($response->is_success){
my %k;
my @a=split(/\n/,$response->content);
foreach my $line(@a){
    next if $line=~/From/;
    my @b=split(/\t/,$line);
    next if(defined($k{$b[0]}));
    $k{$b[0]}=$b[1];
}
map{print STDOUT "$_\t$k{$_}\n"}keys(%k);
}
else{
die 'Failed, got ' . $response->status_line . ' for ' . $response->request->uri . "\n";
}
}
else{
if($response->is_success){
my %map;
my @a=split(/\n/,$response->content);
#############################
    # remove "From   To"   line #
    #############################
shift(@a) if $a[0] =~/^From/i;
my %kk;
map{
    ########################################################
        # UniProt sometimes returns the SAME line twice    #
        ########################################################
    $kk{$_}++;
    #print STDOUT "$_\n" unless $kk{$_}>1;
    next if $kk{$_}>1;
    my @b=split(/\t/);
    push @{$map{$b[0]}},$b[1];
    }@a;

    foreach my $query (keys(%map)) {
    my %sw;
    foreach my $target (@{$map{$query}}) {
    #############################################
            # Does this prot have SWISSPROT IDs?        #
            #############################################
    if ($to eq 'ID')
    {
        my @b=split(/_/,$target);
        $sw{$target}++ if length($b[0])!=6;
    }
    }
    foreach my $target (@{$map{$query}}) {
    ################################################
    # If a protein has a least one SWISSPROT ID    #
    # skip printing its trEMBL ones unless -u      #
    ################################################
    if ($to eq 'ID' and not $allow_trembl and scalar(keys(%sw)) > 0) {
        print "$query\t$target\n" if defined($sw{$target});
    }
    else {
        print "$query\t$target\n";
        }
        }
    }

    ## Count how many were actually found
    my %found;
    foreach my $resp (@a){
    #print STDERR "R:$resp\t";
    foreach my $id (@names){
    #print STDERR "I:$id\n";
    $found{$id}++ if $resp=~/^$id\t/;
        }

    }
    print STDERR scalar(keys(%found)) . " of " . scalar @names . "\n";

    }
    else{
        die 'Failed, got ' . $response->status_line .
    ' for ' . $response->request->uri . "\n";
    }
}
################################################
# This will translate $to and $from values     #
# between the names used in uniprot flat files #
# and those used by the uniprot database       #
################################################
sub translate_terms{
    my $id=shift();
my %hash=(
"AC" => "ACC", "ID" => "ID", "UniParc" => "UPARC", "UniRef50" => "NF50", "UniRef90" => "NF90", "UniRef100" => "NF100", "EMBL/GenBank/DDBJ" => "EMBL_ID", "EMBL/GenBank/DDBJ CDS" => "EMBLD", "PIR" => "PIR", "UniGene" => "UNIGENE_ID", "GI number*" => "P_GI", "IPI" => "P_IPI", "RefSeq" => "P_REFSEQ_AC", "PDB" => "PDB_ID", "DisProt" => "DISPROT_ID", "HSSP" => "HSSP_ID", "DIP" => "DIP_ID", "MINT" => "MINT_ID", "MEROPS" => "MEROPS_ID", "PeroxiBase" => "PEROXIBASE_ID", "PptaseDB" => "PPTASEDB_ID", "REBASE" => "REBASE_ID", "TCDB" => "TCDB_ID", "Aarhus/Ghent-2DPAGE" => "AARHUS_GHENT_2DPAGE_ID", "ECO2DBASE" => "ECO2DBASE_ID", "World-2DPAGE" => "WORLD_2DPAGE_ID", "Ensembl" => "ENSEMBL_ID", "Ensembl Protein" => "ENSEMBL_PRO_ID", "Ensembl  Trscrpt" => "ENSEMBL_TRS_ID", "Ensembl  Genomes" => "ENSEMBLGENOME_ID", "Ensembl    Genomes" => "Protein", "Ensembl Genomes" => "Transcript", "GeneID" => "P_ENTREZGENEID", "GenomeReviews" => "GENOMEREVIEWS_ID", "KEGG" => "KEGG_ID", "TIGR" => "TIGR_ID", "UCSC" => "UCSC_ID", "VectorBase" => "VECTORBASE_ID", "AGD" => "AGD_ID", "ArachnoServer" => "ARACHNOSERVER_ID", "CGD" => "CGD", "ConoServer" => "CONOSERVER_ID", "CYGD" => "CYGD_ID", "dictyBase" => "DICTYBASE_ID", "EchoBASE" => "ECHOBASE_ID", "EcoGene" => "ECOGENE_ID", "euHCVdb" => "EUHCVDB_ID", "EuPathDB" => "EUPATHDB_ID", "FlyBase" => "FLYBASE_ID", "GeneCards" => "GENECARDS_ID", "GeneDB_Spombe" => "GENEDB_SPOMBE_ID", "GeneFarm" => "GENEFARM_ID", "GenoList" => "GENOLIST_ID", "H-InvDB" => "H_INVDB_ID", "HGNC" => "HGNC_ID", "HPA" => "HPA_ID", "LegioList" => "LEGIOLIST_ID", "Leproma" => "LEPROMA_ID", "MaizeGDB" => "MAIZEGDB_ID", "MIM" => "MIM_ID", "MGI" => "MGI_ID", "NMPDR" => "NMPDR", "Orphanet" => "ORPHANET_ID", "PharmGKB" => "PHARMGKB_ID", "PseudoCAP" => "PSEUDOCAP_ID", "RGD" => "RGD_ID", "SGD" => "SGD_ID", "TAIR" => "TAIR_ID", "TubercuList" => "TUBERCULIST_ID", "WormBase" => "WORMBASE_ID", "WormBase Trscrpt" => "WORMBASE_TRS_ID", "WormBase Protein" => "WORMBASE_PRO_ID", "Xenbase" => "XENBASE_ID", "ZFIN" => "ZFIN_ID", "eggNOG" => "EGGNOG_ID", "HOGENOM" => "HOGENOM_ID", "HOVERGEN" => "HOVERGEN_ID", "OMA" => "OMA_ID", "OrthoDB" => "ORTHODB_ID", "ProtClustDB" => "PROTCLUSTDB_ID", "BioCyc" => "BIOCYC_ID", "Reactome" => "REACTOME_ID", "CleanEx" => "CLEANEX_ID", "GermOnline" => "GERMONLINE_ID", "DrugBank" => "DRUGBANK_ID", "NextBio" => "NEXTBIO_ID"
);
defined($hash{$id}) ?
return($hash{$id}) :
return($id) ;
}

sub get_response{
    my $params=shift;
my $response = $agent->post("$base/$tool/", $params);
    while (my $wait = $response->header('Retry-After')) {
    print STDERR "Waiting ($wait)...\n";
sleep $wait;
    $response = $agent->get($response->base);
}
return($response);

}

sub usage{
    $0=~s/.+?([^\/]+)$/$1/;

    print STDERR <<EndOfHelp;

$0 will connect to UniProt and map IDs from the input file to ehatever ID was requested.

USAGE:
       $0 [options] <ID_LIST>

OPTIONS:
      -h : Print this help and exit.
      -f : From ID type. This is the type of ID found in the input file.
      -F : Format. See http://www.uniprot.org/faq/28#batch_retrieval_of_entries for details.

      -l : When a requested ID maps to many output IDs, print only the first (probably the best) match.
      -L : List all recognised ID types.
      -u : Print unreviewed (trEMBL) prots as well. By default, if an ID maps to SWISSPROT and trEMBL IDs
           only the SWISSPROTs will be printed, use this flag to print trEMBL as well.
      -t : To ID. This is the ID type that the input IDs will be mapped to.
      -v : Verbose. Print progress messages to STDERR.
EndOfHelp

    exit(0);

}
sub list_ids{
    open(HELP, "| less") ;
    print HELP <<EndOfHelp;

Database            Abbreviation
===================== UniProt =============================
UniProtKB AC/ID :   ACC+ID
UniProtKB AC    :   ACC
UniProtKB ID    :   ID
UniParc         :   UPARC
UniRef50    :   NF50
UniRef90    :   NF90
UniRef100   :   NF100

=============== Other sequence databases ==================
EMBL/GenBank/DDBJ:  EMBL_ID
EMBL/GenBank/DDBJ CDS:  EMBLD
PIR     :   PIR
UniGene     :   UNIGENE_ID
GI number*  :   P_GI
IPI     :   P_IPI
RefSeq      :   P_REFSEQ_AC
HGNC name       :       GENENAME
================== 3D structure databases =================
PDB     :   PDB_ID
DisProt     :   DISPROT_ID
HSSP        :   HSSP_ID

============ Protein-protein interaction databases ========
DIP     :   DIP_ID
MINT        :   MINT_ID

============ Protein family/group databases ===============
MEROPS      :   MEROPS_ID
PeroxiBase  :   PEROXIBASE_ID
PptaseDB    :   PPTASEDB_ID
REBASE      :   REBASE_ID
TCDB        :   TCDB_ID

================== 2D gel databases =======================
Aarhus/Ghent-2DPAGE:  AARHUS_GHENT_2DPAGE_ID
ECO2DBASE   :   ECO2DBASE_ID
World-2DPAGE    :   WORLD_2DPAGE_ID

============= Genome annotation databases =================
Ensembl         :   ENSEMBL_ID
Ensembl Protein :   ENSEMBL_PRO_ID
Ensembl Trscrpt :   ENSEMBL_TRS_ID
Ensembl Genomes :   ENSEMBLGENOME_ID
Ensembl Genomes :   Protein
Ensembl Genomes :   Transcript
GeneID          :   P_ENTREZGENEID
GenomeReviews   :   GENOMEREVIEWS_ID
KEGG        :   KEGG_ID
TIGR        :   TIGR_ID
UCSC        :   UCSC_ID
VectorBase  :   VECTORBASE_ID

============ Organism-specific gene databases =============
AGD     :   AGD_ID
ArachnoServer   :   ARACHNOSERVER_ID
CGD     :   CGD
ConoServer  :   CONOSERVER_ID
CYGD        :   CYGD_ID
dictyBase   :   DICTYBASE_ID
EchoBASE    :   ECHOBASE_ID
EcoGene     :   ECOGENE_ID
euHCVdb     :   EUHCVDB_ID
EuPathDB    :   EUPATHDB_ID
FlyBase     :   FLYBASE_ID
GeneCards   :   GENECARDS_ID
GeneDB_Spombe   :   GENEDB_SPOMBE_ID
GeneFarm    :   GENEFARM_ID
GenoList    :   GENOLIST_ID
H-InvDB     :   H_INVDB_ID
HGNC        :   HGNC_ID
HPA     :   HPA_ID
LegioList   :   LEGIOLIST_ID
Leproma     :   LEPROMA_ID
MaizeGDB    :   MAIZEGDB_ID
MIM     :   MIM_ID
MGI     :   MGI_ID
NMPDR       :   NMPDR
Orphanet    :   ORPHANET_ID
PharmGKB    :   PHARMGKB_ID
PseudoCAP   :   PSEUDOCAP_ID
RGD     :   RGD_ID
SGD     :   SGD_ID
TAIR        :   TAIR_ID
TubercuList :   TUBERCULIST_ID
WormBase    :   WORMBASE_ID
WormBase Trscrpt:   WORMBASE_TRS_ID
WormBase Protein:   WORMBASE_PRO_ID
Xenbase     :   XENBASE_ID
ZFIN        :   ZFIN_ID

================= Phylogenomic databases ==================
eggNOG      :   EGGNOG_ID
HOGENOM     :   HOGENOM_ID
HOVERGEN    :   HOVERGEN_ID
OMA     :   OMA_ID
OrthoDB     :   ORTHODB_ID
ProtClustDB :   PROTCLUSTDB_ID

=================== Enzyme and pathway ====================
BioCyc      :   BIOCYC_ID
Reactome    :   REACTOME_ID

================= Gene expression databases ===============
CleanEx     :   CLEANEX_ID
GermOnline  :   GERMONLINE_ID

========================== Other ==========================
DrugBank    :   DRUGBANK_ID
NextBio     :   NEXTBIO_ID

EndOfHelp
close(HELP);
    exit(0);
}

If you have too many IDs to do in a single batch, you could either simply split the file into smaller ones and run the script on each file in a loop or modify the script to do it.
Split the input:
$ wc -l names
$ wc -l names
500000 names
## split into 500 files of 1K lines each
$ split -l 1000 names chunk
## iterate over each chunk of names
$ for i in chunk*; do uniprot_map.pl -f ID -t UniParc "$i" >> mapped; done

Thank you @terdon for the detailed answer and your script! Let me try it out!

Seems both of the solutions require ID list < 2Mb. I will try to look into the XML file of UniProt. Please check out: https://blog.liang2.tw/posts/2018/01/read-uniprotkb-xml/

@JohnnyTam just split your file and run the script in a loop. Or change the script to run them in batches.

@JohnnyTam see updated answer for how to split the input. Alternatively, if you can give me your list of names, I can have a go at modifying the script to handle them.

Ops, I forgot I can split. It worked perfectly! Thank you very much!


text:
Where can we download a curated dataframe of the full archive of UniProt?
I want to download a curated dataframe of the full archive of UniProt.
Any format as long as it is a table (.csv, .tsv, SQL etc.).
Solution 1: Download the .fasta and do the formating myself -> very slow + information incomplete.
Solution 2: Use API like BioServices -> still very slow but more information.
Do we have an already-curated dataframe containing as complete as possible the entry + metadata of UniProt that can be downloaded right away?

Answer from @matteo-ferla, converted from comment:
There are many tables in the downloads page such tables enumerating each protein in a species proteome with basic stats, but not all the data. The uniprot records are very structured so would convert poorly to a table (even filling cells with lists/arrays and dictionaries/objects)… so a fuller table would be a mess.


text:
mpiblast alternatives?
I am looking for a sequence alignment software capable of pairwise alignment of a large number of protein sequences (10^7).
I tried to go for mpiblast: https://wiki.canterbury.ac.nz/display/UCHPC/Installing+mpiBLAST
but seems their website was down.
Do you know any alternatives to mpiblast which is also highly efficient in protein sequence alignment?
Thanks!

The Diamond protein aligner is a very popular tool for pairwise protein alignment which is much faster than blastp while having comparable sensitivity. Diamond also supports the use of using multiple compute nodes (similar to mpiBLAST). You can view the details of the new Diamond 2.0 release in this recent Nature Methods paper (source of the image as well).

...we aligned all 281 million protein sequences from the NCBI nr database against the UniRef50 database, which consists of 39 million sequences, using DIAMOND (v2.0.7) in ultra-sensitive mode...

Throckmorton thanks a lot for the suggestion! I'll have a look at this aligner now.


text:
Mapping PDB ID + chain ID to UniProt ID
This question was addressed here.
Unfortunately, including bioservices, the solutions pypdb and map_pdb_to_uniprot do not work for now.
Besides, I do not need residue level mapping as one of the suggested solution.
Does the failure of known methods have anything to deal with the recent update of the UniProt database?
Please let me know if you know a workaround.
Thank you!

https://www.ebi.ac.uk/pdbe/docs/sifts/quick.html

Please find the solution here:
https://www.biostars.org/p/9540519/#9540826


text:
Site Mutations in a 3D protein structure to locate the core of mutations
I have a pdb file for a protein (3D structure predicted by a tool) and around hundreds of known point mutations in the protein. I want to visualize the protein with all the mutations to see the hub of or core where the mutations are mostly. Is there are tool that does this? I would give the structure in a pdb file and a list of point mutations myself.

Mutations can have different effects and severity. "Visualising" them is very powerful for discussion purposes, but quantifying them may be better, see this past Q. You can open your protein in PyMOL and select the residues in the command line part (select interested, resi 1+2+3+5) and then show them (show sticks, interested)... but the problem then the format you have the mutations in influences how you make the concatenated indices...

Also this would be showing the wild type residue, but showing all mutations in a single structure makes no sense in most cases.

@MatteoFerla Clustering mutations could reveal lots about the protein providing one has the functional impact of each mutation (this cluster are benign, this cluster causes misfolding etc). Indeed a region of no mutations could indicate constraint.

Clustering mutations could reveal lots about the protein providing one has the functional impact of each mutation (this cluster are benign, this cluster causes misfolding etc). Many tools are emerging to explore this and there are plenty of established tools that can handle this.
VarMap is a tool that is able to map genomic variant coordinates to the best available structure.
PyMol can open PDB files and you can programmatically or manually mutate each position using the mutagenesis tool.


text:
efetch multiple records at once
If I want to efetch two records, I can do the following:
efetch -db nuccore -id AB610939,AB610940 -format fasta

I have hundreds of records to fetch. Can I put all the IDs in a file and efetch. If so, how?

Yes, efetch can be given a file with IDs. See efetch -h:
$ efetch -h | grep -m1 file
  -input         Read identifier(s) from file instead of stdin

So you can pass your file of IDs (e.g. ids) directly:
$ cat ids
AB610939
AB610940

$ efetch -db nuccore -input ids -format fasta
>AB610939.1 Synthetic construct RNA, spike-in microarray control RNA CRM 6204-a, clone: 6204-a-500-1
CTCGACTAGTTAATACGGTACAGGATAACCGATCGGCTTGCAACATAACGGCGTTAAGAATGCGGGAGTG
CAGTTTCCGATTCTCACATCAATCGCCAATAAGGCCTTGTCGCAATATAGACTCAACGGTTCTAGTAGCT
GATCGGTATTACGTGACGCAACCGATTAGACATGCACAATTCCTTGGTCGCTATACTACGGAAATCGTCA
GGTACTATAACCCGTCGCAGGCCTAATACGTGTCGTCACATCGCCAACCTATCGTCAGTCGGAAAGACGT
TGCTGTCTACCATCGAAACTATTTACCGCTCCGAGATTCACGAGTACGAACTCACGAGGAAGTTGCCCTA
TGTAAGGTATCACTCCAGGTACTGCGCCGATAGTACCAGGTGATCAAACGGTTGCAAGAAGGCCACGACG
TATCGGGCTCTTTAGACGTACGCTCGAGATTAAACGCGCACTGATTCACTTTAGCCCGGAATGTCTCGGT
GCGATGTAGA
>AB610940.1 Synthetic construct RNA, spike-in microarray control RNA CRM 6204-a, clone: CRM6204-a-500-2
AGACTAAATCTCGGCGTCGGTTCATACGCGCGATCGTTTGCTGTCAGGGCATACTCGAATCCGGACTCCG
ACAATTATAGGCCATCCTGAATAGCCGATCATGCGAGTCACGATAAGGCAGGCTCTGCGATATCCCGATA
TACTGGAGAAGCTGAATCCCACCTAGAGCGAACTGTCAGAGGATCGACCTCAGGCTCGCTATCATCATAA
CGGCGGACGACCTGTGTCACATTCCGAACGCTACGTGACGATATTATCTGTCGAAAGGCATAGAACGCCG
GTCAATATCCTGCGGCATTCTCTTTATCACCGGCTATAACTACTAGGTTCCGCAGATATAGACTGCGCAC
GGAACATGTAGATAGATCGAGTAGGGTAGCGATTTAACGACTCGACTTACAGACAGAGACGTAGAACGTC
AGACGAGTGGTATGCCCACCAGAGGCGATACAGGCTGTACCTGCGTAGCACTAGAGTCGTGCGTCATGCG
GACCCTATCT

If that were not possible, another option would be to put all your IDs as a comma separated string. If your IDs are in a file named ids, one per line, you can convert that to a comma-separated string like this:
$ tr '\n' , < ids | sed 's/,$/\n/'
AB610939,AB610940

That assumes GNU sed. If your sed doesn't handle \n, you can use this instead:
printf "$(tr '\n' , < ids | sed 's/,$//')" "$(cat ids)"

Armed with this, you can use command substitution to pass the comma-separated list to efetch:
efetch -db nuccore -id $(tr '\n' , < ids | sed 's/,$/\n/') -format fasta

If there are limits to te number of IDs you can fetch at once, or as may well be the case, if this results in a command that is too long, then you can use a loop as suggested in another answer.

However, be aware that while looping seems to work for efetch, for reasons that are not clear to me, the esearch developers have chosen to make it consume stdin. This means that the obvious, straightforward approach fails with esearch. See https://unix.stackexchange.com/q/682748/22222 for details.

If you're happy looping through each id in turn and each id in the single file (idfile.txt) is separate by a newline:
while read -r x; do
  efetch -db nuccore -id  "$x" -format fasta
done < idfile.txt > fetchout.fa

Speculation If the ids are not separated by a newline but a comma on a single line that might give the simultaneous efetch you're wanting. I strongly suspect that each id would need to be surround in quotation marks in such a file. Again this is just speculation and the above solution is preferred IMO.

Sorry, my bad. That issue seems to be specific to esearch, I cannot reproduce it with efetch. I deleted my comment.

@terdon world-first from the BioInfo 'nix tools police ;-)


text:
How to identify [hyper]thermophilic proteins by BLAST search?
I need to find homologs of certain proteins by BLAST search or similar. But I want filter (select) homologs from organisms that are classified as thermophilic or hyperthermophilic. I am wondering if there are suggestions for me how approach to this. Thank you.

It seems like the main issue is finding a list of species to filter your blast results on. I will therefore address this list of species issue.
There are some resources such as Microbe Directory that contain organismal characteristics similar to what you describe. Their main site appears nonfunctional but their github still works and appears to have the full database. This database includes an extreme_environment field, which appears to be what you are going for.
I feel like PATRIC used to have something similar but I can't locate it now.
For more ideas you might look at the Microbe Directory paper to see if they cite other similar databases.

Warning: Hearsay

As the others have said, filtering by a list of taxonomy ids is the
only one. However, there is one extra consideration if this is for wet
lab applications:  expression yields even when codon optimised may be
affected by phylum choice.
Archaeal protein encoded by codon-optimised genes for E. coli generally perform less well than the equivalent bacterial thermophilic protein. I have not read a study on it, but I have heard its tale over and over and had issues myself. GC-content and the amino acid composition do correlate, but I have heard of good results in E. coli with Geobacillus (high-GC thermophilic firmicute) as with Bacillus (high GC mesophilic firmicute) after all. I don't know if its isoelectric points close to pH 7 or what.
But mid GC-content bacterial thermophiles, such as Aquifex or some clades of Betaproteobacteria, may be a better call than Archaea.
So as you need to filter your blast search by a series of phyla, you may want to be selective if for expression in a different host...

You could be onto something, I agree - but I suspect someone would have already noticed. I had a brief look at this problem and the level of convergent evolution between bacterial families is extraordinary looking at the results of metagenomic studies in thermal environments. @MaximilianPress had the easiest way forward. There's a lot more going on, particularly at the protein level.

I don't know, I do honestly think it's one of those many industry anecdotes: things that are of importance to white biotech but not to the academic publication machinery. I'll grant that they mostly involve the world of fermentors —e.g. E. coli clumps with added thiamine which messes up fermentor OD readings, which if as a reviewer I saw as a study I'd end up rejecting for not being scientific enough :D.

Okay, I've probably misread the aims ... as an evolutionary question its hard.

No no. I guessed as the aims: I speculated the OP was interested in searching for high T_M proteins from an engineering perspective as I've done so myself. An enzyme's K_M is close to physiological concentration of the substrate and K_M generally decreases with decreasing temperature, so it's a classic cheat to get enzymes with lower K_Ms.

Got it .. thanks

I recommend use the local BLAST search. You can form the request for desired taxa in Entrez or through NCBI website & download fasta-file onto your disk.

You know any taxa that are classified as hyperthermophilic except the obvious ones? Doing BLAST is trivial. For me the important thing to find the right databases.

I am scarcely sure that such database exists. There was ExtremeDB but it is inaccessible now. You could try to write a letter to the authors, maybe they changed the domain name.


text:
How to display novel genome assemblies or uncommon genome assemblies using the UCSC Genome Browser?
I want to display E.coli BW25113 (GenBank: CP009273.1) strain in UCSC browser. This strain is not listed in http://microbes.ucsc.edu/ browser. How can I display E.coli BW25113 assembly in the browser?

I found two methods:

You want to visit this page for instructions for novel assemblies: http://genomewiki.ucsc.edu/index.php/Assembly_Hubs .
However, before doing that I highly recommend that you check Genome Assembly Hub first:
https://genome-test.gi.ucsc.edu/gbdb/hubs/genbank/ . There is a good chance that your genome or strain is listed there. Here is an example of how you can do that. Example specifically shows how you can display E.coli BW25113 strain in the browser.

Instructions:

Go to NCBI Nucleotide: https://www.ncbi.nlm.nih.gov/nucleotide .
Enter your strain name: E.coli BW25113 (GenBank: CP009273.1)
Find Biosample ID on the page, located under DBLINK: BioSample: SAMN03013572
Go to GenBank assembly hubs: https://genome-test.gi.ucsc.edu/gbdb/hubs/genbank/ .
Select bacteria. This will take you to "Bacterial assembly hub"
at https://genome-test.gi.ucsc.edu/gbdb/hubs/genbank/bacteria/bacteria.html .
Look for E. coli hubs. At the time I looked #15-18 contain E. coli strains.
Click on each hub and search for SAMN03013572 using browser search function (Command+F on macOS). I found my genome in hub #16.
Click on "common name" in the column-3. This will take you to
https://genome-test.gi.ucsc.edu/cgi-bin/hgGateway?hubUrl=http://genome-test.gi.ucsc.edu/gbdb/hubs/genbank/bacteria/16/hub.txt&genome=GCA_000750555.1_ASM75055v1&position=lastDbPos , where FindPosition window will be populated with your genome Assembly.
Click on Go. This should take you to the UCSC Genome Browser with your genome assembly.

Update: The step-9 is no longer working since UCSC genome browser resources described above were experimental prototype proof of concepts
that have since been discontinued according to the Genome Browser Team. Once I find an alternative method, I will post it here.

I am not sure if this qualifies as an official answer. I have contacted UCSC bioinformatics team. They were were very kind and told me that they would add any track that was requested by users. After my request, they kindly added the track for E.coli BW25113. If you want your reference tracks to be added you could submit a request at this website:
Geneme assembly search and request page.
Then, at the upper right corner click on "Can't find your assembly?" button and submit your request. See the screenshot below from the webpage.


text:
Edit end column in GFF3 file
I am trying to edit the "end" of my sequences in a GFF3 that contains 10000 sequences (so I do not want to do it manually). I want to add 30 residues to every single fragment. How can I do this?
Currently, I have a file that looks like this:
##gff-version 3
sp_P45452_MMP13_HUMAN   SignalP-5.0 signal_peptide  1   19  0.994300    .   .   .
sp_Q8R2H5_PHLD_RAT      SignalP-5.0 signal_peptide  1   23  0.986089    .   .   .

I would like to edit it so that it looks like this:
##gff-version 3
sp_P45452_MMP13_HUMAN   SignalP-5.0 signal_peptide  1   49  0.994300    .   .   .
sp_Q8R2H5_PHLD_RAT      SignalP-5.0 signal_peptide  1   53  0.986089    .   .   .

Thank you!

That's what awk excels at! Simply set the input (-F '\t') and output (-v OFS='\t') record field separators to tab, and add 30 to the value of the fifth field on lines that don't start with a #:
$ awk -F'\t' -v OFS='\t' '/^[^#]/{$5+=30}1' file.gtf
##gff-version 3
sp_P45452_MMP13_HUMAN   SignalP-5.0 signal_peptide  1   19  30.9943 .   .   .
sp_Q8R2H5_PHLD_RAT  SignalP-5.0 signal_peptide  1   23  30.9861 .   .   .

yep this is the good stuff

Thank you so much! I am new to awk and this was super helpful!


text:
file order after nextflow groupTuple
Below is the code from nextflow site (https://www.nextflow.io/docs/edge/operator.html#grouptuple).
For the output, the first is A B C, and the second is C, A.
Is there a way to make the output like this
[1, [A, B, C]]
[2, [A, C]]
That is, the stuff inside [] is in some order.
Channel
    .of(
        [1, 'A'],
        [1, 'B'],
        [2, 'C'],
        [3, 'B'],
        [1, 'C'],
        [2, 'A'],
        [3, 'D']
    )
    .groupTuple()
    .view()

[1, [A, B, C]]
[2, [C, A]]
[3, [B, D]]

You can change the sorting criteria for the grouped items using the sort keyword. To perform a natural sort, just use sort: true. For example:
Channel.of(
    [1, 'A'],
    [1, 'B'],
    [2, 'C'],
    [3, 'B'],
    [1, 'C'],
    [2, 'A'],
    [3, 'D']
)
.groupTuple(sort: true)
.view()

Results:
$ nextflow run main.nf
N E X T F L O W  ~  version 22.04.0
Launching `main.nf` [special_darwin] DSL2 - revision: 2e641ccfd9
[1, [A, B, C]]
[2, [A, C]]
[3, [B, D]]

Otherwise, a custom sorting criteria can be specified using a closure.


text:
pacbio demultiplexing issue
I was wondering if there is a way to provide pacbio lima demultiplexing only a barcode, and let  lima demux sample even if only one end has the barcode, not both?

Can you please explain more about the context of your problem? What are you trying to do with the demultiplexing; what is the ultimate aim of the analysis?

Some examples of what your data look like might help

PacBio lima need both 5'end and 3' end barcodes. You can have a unique barcode in one end and multiple on the other end; but it need two barcodes.


text:
NMF - Subscript out of bounds
I've calculated the cophenetic coefficients for the NMF of gene expression data, but is giving error when performing the clustering information step as given below:
clustered_data <- NMF_estimate[[num.opt - 1]][[2]]

Error in NMF_estimate[[num.opt - 1]] : subscript out of bounds

Here num.opt is num.opt <- coph[which(coph[,2]==max(coph[-1,2])),1] which gives 2, 5 as outputs, and NMF Estimate is a list of 6 elements.
Please suggest a solution as soon as possible.


text:
Linking UK biobank patient IDs with their corresponding genetic data
I am analyzing ukbb data and I have patient IDs for both cases and controls. Now I want the genetic data of all patients based on their IDs.
I am trying to test with chromosome 5 as I have both bgen and sample file. I am using plink 2 for the conversion of bgen and sample file
plink2 --bgen ukb_chr5_v3.bgen ref-first --sample ukb22828_c5_b0_v3_s487239.sample --make-bed --out testfile
I am getting three files as output like .pgen .psam .pvar.
This is how .psam looks like,

FID
IID
sex

482828
482282
1

and .pvar file has five columns, chromosome, position, ID, ref allele, alt allele.
I am still confused which IDs have the corresponding alleles/genotypes.
Please guide me in further analysis. My goal is find genetic data of specfic IDs which I am looking for.
Thanks

The genotypes are in the .pgen file, but encoded as a binary, so you can't query it directly. What are you going to do with the genotypes?

My aim is to calculate polygenic risk scores using pre-calculated weights for each individual and for that I possibly need effect allele count (0,1, or 2)