text:
How to download multiple proteomes at once?
I'm looking for a way to download multiple proteomes at once from one clade, as fasta. Probably from the NCBI, because it looks more user-friendly than Ensembl, but Ensembl is okay too. In the best case scenario, I would like to constrain the download to a specific clade (eg. Actinobacteria, Chloroflexi etc.), maybe with a specific RefSeq Category (eg. reference/representative genome) or even depending on the genome assembly level. Command line would probably be the best option, but if there's a direct download option it would be nice to know it as well.
Does anything like this exist? Help would be much appreciated!

The FTP has it: ftp://ftp.ncbi.nlm.nih.gov/genomes/ therein look for the fa.gz file

You could try download-refseq-genomes. It will fetch all genomes from the NCBI FTP server that are in a specific subtree of the phylogeny.
For example, downloading the amino acid sequences from all proteins in each genome assembly for all species belonging to Chloroflexi (NCBI taxon ID = 200795):
download-refseq-genomes.pl -t faa 200795

This will download all assemblies with status=Complete. Use option -a to download all assemblies with any status.

You can use NCBI resources to download these data directly. For the entire proteome, you can download the sequences of all proteins annotated on the genome assembly in FASTA format. In the NCBI FTP directories, these files have the names with the suffix _protein.faa.gz. More details about the files in the FTP directories are described in the README.
Command-line method using EDirect
This method uses Entrez Direct to first query the Assembly database for all available genomes, parses the output to extract the ftp path and downloads data using wget command.
esearch -db assembly -query 'txid200795[organism] AND lastest_refseq[filter]' \
  | esummary \
  | xtract -pattern DocumentSummary -element FtpPath_RefSeq \
  | while read -r url ; do
    path=$(echo $url | perl -pe 's/(GC[FA]_\d+.*)/\1\/\1_protein.faa.gz/g') ;
    wget -q --show-progress "$path" -P protein_data ;
  done

Here, I am restricting the results to latest RefSeq data by using the latest_refseq[filter] but you can use any other Entrez filters to further narrow down your results. Along the same lines, you can modify the code minimally as shown below to download multiple file types:
esearch -db assembly -query 'txid200795[organism] AND latest_refseq[filter]' \
  | esummary \
  | xtract -pattern DocumentSummary -element FtpPath_RefSeq \
  | while read -r url ; do
    base_path=$(echo $url | perl -pe 's/(GCF_\d+.*)/\1\/\1_/g') ;
    for suf in genomic.fna.gz protein.faa.gz ; do
      wget -q --show-progress "$base_path$suf" -P ncbi_data ;
    done ;
  done

GUI method
You can use the same query as shown above in the NCBI Assembly portal and then use the 'Download Assemblies' button to pick and choose the type of files you would like to download. By default, if you do not select any assemblies in the results, this interface downloads files for all results.


text:
Using variables with fasterq-dump?
I am trying to download multiple fastq files from the SRA NCBI database which, conveniently enough, have their IDs in a range. Building a simple loop doesn't seem to do the trick, as bash appears to interpret my variable as an additional command? Help would be appreciated. This is what I have for now:
for i in {256..345};do
    sra_code = "DRR163""$i"
fasterq-dump "$sra_code" -p -v --outdir ./output/ -e 8

done
```

You have an obvious error in your shell script that should have given you an error message:
sra_code = "DRR163""$i"

Should be:
sra_code="DRR163""$i"

You cannot have spaces around the = sign in a variable assignment. The line sra_code = "DRR163""$i" means "run the command sra_code with the arguments = "DRR163""$i". Since there is no such command, you see an error.
Here's a working version of your command:
for i in {256..345}; do
    sra_code="DRR163$i"
    fasterq-dump "$sra_code" -p -v --outdir ./output/ -e 8
done

which, conveniently enough, have their IDs in a range

I think it's a bit risky to rely on that - for this particular case they are not going to change but you may need to download more samples in the future.
I would suggest avoiding sra and fastq(er)-dump altogether and use ena instead.
In your case, it seems you want the runs in project PRJDB7736. From here https://www.ebi.ac.uk/ena/browser/view/PRJDB7736 you can get a tabular file (see link in Download report ... TSV) with the sample characteristics. Column fastq_ftp gives you the ftp link to the fastq files that you can pass to wget or curl.
You can also get the same report programmatically with:
curl -s -X GET "https://www.ebi.ac.uk/ena/portal/api/filereport?accession=PRJDB7736&fields=run_accession,sample_description,library_name,library_layout,fastq_ftp&result=read_run" > ena.tsv

You should edit the fields=... list to suite your needs. The fields available are those from:
curl -X GET "https://www.ebi.ac.uk/ena/portal/api/returnFields?dataPortal=ena&result=read_run"

I don't know what you think... I find this route much more useful and transparent than sra's.

source: Stack Exchange

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

sort: votes


text:
Software for taxonomic assignment?
I have a couple of hundred bacterial sequences of 2-30 genes of interest each, recovered from metagenomics. None of them encode rRNA. Normally I'd just BLAST the one gene I already know to be reliable phylogenetic marker, where present, but the scale of my current dataset makes it a bit daunting.
Is there any software for taxonomic assignment that you could recommend? It's fine if it requires command line use. There are so many options listed online that I didn't even know where to start... Please and thank you.

You could take a look at kraken2 It should fit your purpose.

yes, kraken2 performed well in an independant benchmark https://lemmi.ezlab.org/

It does work well, but it's less accurate at lower taxonomic ranks (as most tools are)

If you like blast, maybe you can try DIAMOND - the authors claim it's much faster than BLAST.

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


text:
How to measure the total size of a fastq file in base pairs?
Or Kbps/Gbps. It feels like it should be conceptually very simple, but I can't seem to figure out the right combination of keywords to find it via my search engine. Help would be appreciated!
I have BBMAP, SRAtoolkit and MEGAHIT already installed, and also use bash. I'd be very happy if this can be answered with software that I already have, but if not that's perfectly fine.

I've been using this:
cat file.fastq | paste - - - - | cut -f 2 | tr -d '\n' | wc -c
Explanation :

paste - - - - : print four consecutive lines in one row (tab delimited), to merge the info for each read

cut -f2 : print only the second column, to access the sequence after the paste

wc -c : count the characters

tr -d '\n': to remove from count the eventual newline characters

(a tip for your googling: try search for "counting number of bases in fastq file")

seqkit stats fastq.gz. Get SeqKit from https://bioinf.shenwei.me/seqkit/usage/

Using Perl:
Print the sequence lines (line number 2, 6, 10, etc). Remove the newlines with chomp. Count the bytes (here, bases) using wc -c:
perl -ne 'if ( $. % 4 == 2 ) { chomp; print; }' | wc -c

Using a fastq parser in Python (for instance the one provided by mappy), it is rather straightforward to compute the sum of the lengths of the sequences present in a fastq file:
$ python3 -m pip install mappy  # Only if mappy is not installed yet
$ python3 -c 'from mappy import fastx_read; print(sum(len(seq) for (_, seq, *_) in fastx_read("SRR077487_2.filt.fastq.gz")))'
2386161200

python3 -c '<insert some python code here>' uses the -c option of python3 that tells Python to execute the provided code (instead of starting an interactive interpreter, or running code present in a file).
Here the code has two steps (separated by a ;, in a standard script, we could use a new line instead): first making the fastx_read function available (because it is part of an optional module), second displaying the total length of the sequences.
fastx_read("path_to_a_file") generates (name, sequence, quality) triplets when a fastq file is provided, or (name, sequence) pairs when a fasta file is provided.
The (_, seq, *_) syntax is a form of "tuple unpacking" where we store the second element (the sequence) in a seq variable.
The _ in is just a way to ignore the first element in the triplet/pair. The *_ is a way to ignore whatever is after the second element in the triplet/pair.
len(seq) for (_, seq, *_) in fastx_read("SRR077487_2.filt.fastq.gz") is an expression representing the lengths (obtained with the predefined len function) of the sequences present in the file "SRR077487_2.filt.fastq.gz".
So the code uses the predefined sum function to compute the sum of values that are the lengths of the sequences extracted from the fastq file, and this is displayed using the print function.
This should also work with either fasta or fastq format, and with or without gz compression.
I tested more parsers on a slightly different problem here: https://bioinformatics.stackexchange.com/a/380/292 and you can find inspiration with the other approaches proposed in the other answers there.

The number of bases in a fastq file can be counted in bash with awk and wc
awk 'NR % 4 == 0' ORS="" fastqfile|wc -m

The awk code prints every fourth line (which is actually the quality scores, but that doesn't matter here). wc -m returns the total number of characters.
Edited to add ORS="" to prevent count of newlines (the double quotes are not strictly required)

With bbmap:
reformat.sh </path/to/fastq>


text:
Comparing groups of 16S rRNA sequences - how to?
I have several groups of 16S rRNA sequences associated with taxonomic groups (let's call them A, B and C). The sequences within each group have a common ancestor and are on average more closely related to sequences within the same group. The groups are also related, but less.
Ideally I'd like to have a table like this (where percentages represent the average similarity):
           Group A     Group B     Group C
Group A      99%
Group B      97%         98%
Group C      94%         95%         99%

Is there software that can calculate these percentages for me? I'm thinking that I could set up a many-against-many BLAST search and then try to use pandas in python in order to parse the resulting table...But if something like this already exists, I'd rather not reinvent the wheel. Thank you!

If you want to make an all vs all with blast this can be helpful. You can create a blast database with your sequences and then search against itself.
makeblastdb -in 16S_sequences.fasta -dbtype nucl -out my_16S_sequences_db

blastn -db my_16S_sequences_db -query 16S_sequences.fasta -outfmt 6 \
       -out allvsall_16S.tsv -num_threads "$(nproc)"

If you want a percent identity matrix for a pairwise alignment you can try with this approach: https://www.biostars.org/p/220154/


text:
Comparing gene abundances between metagenomes
My workflow until now:
Find fragments of a marker gene in unassembled metagenomes > download and assemble metagenomes > recover the gene neighborhood / gene set of interest
Right now I have a rough estimate of how abundant these genes are by the 'depth' noted on the assembled contigs (I use MEGAHIT for assembly). I was wondering if there's a more thorough/proper way to do this. I would like to compare the abundance of specific genes between a) samples in the same study, and b) different studies. I imagine that the size of the individual metagenomes should be considered in both cases, but point b) might add additional difficulties such as different sequencing techniques. I would appreciate your insights.

I would avoid using assemblies to answer this question, as there's no guarantee that you will be able to assemble your genes of interest; you can however estimate their abundance even if they are relatively rare.
How I understand your question as being one of estimating the abundance of either some specific genes (e.g. butyrate metabolism genes) or all genes in a microbial community across multiple samples for comparative purposes. In other words, not 16S or marker gene analysis for the purposes of estimating organismal abundance, which is a rather different problem (though in that case I would still not use an assembly).
A more standard workflow is:

align metagenomic reads against some existing database of genes annotated appropriately.
estimate the number of reads aligning against some gene or ortholog in some way (using for example KEGG Orthology or similar).
use counts from (2) as input to some statistical procedure, possibly summarizing across functional categories.

Some examples of how this has been done are here, here, here. I am sure that there are more recent/relevant references but I haven't been following the field closely in the last few years.


text:
Voom transformation of RNA seq raw counts data
I have RNA seq data of raw gene counts that I want to transform for linear modelling. I am trying to voom transform, to do a weighted analysis.
Data frame: prac_count_10
gene         S1     S2    S3   S4    S5
ENG0000456   0      1     10   145   24
ENG0000458   7      2     9    0     0
ENG0000657   76     12    56   10    2
ENG0000689   0      0     0    3     5

Code:
prac_prac_10 <- voom(prac_count_10[,2:5], design=NULL, lib.size=NULL,
  normalize.method = "none", span = 0.5

However this doesn't give an output with the gene sample names and logCPM values. Would like an output with gene id's and logCPM values for linear modelling.
Also tried:
cpm <- cpm(prac_count_10)
lcpm <- cpm(prac_count_10, log=TRUE)

````````````

Try to use DElist() function before you transform, and also make rownames first.
counts <- prac_count_10[,-1]

rownames(counts) <- prac_count_10[,1]

DGE1 <- DGEList(counts)

And then continue with your voom...

Thank you. what package is DElist in?

That would be edgeR.

If I do:   prac_count <- data.frame(prac_count, row.names=1), counts<- prac_count[,1], rownames(counts) <- prac_count[,1], DGE1 <-DGEList(counts)  If doesn't work... have I missed a step? Thank you so much

That's not the code I gave in my answer, is it?

Don't mention it! Good luck with it.


text:
WGCNA error codes in Network Construction & moduel detection
I am using WGCNA and code has worked with no problems with some of my datasets however with most recent attempt when completing gene network & identification of modules. I get following error code...
net__30= blockwiseModules(t5, power=6,
                              TOMType = "signed", minModuleSize=30,
                              reassignThreshold =0, mergeCutHeight = 0.25,
                              numericLabels = TRUE, pamRespectsDendro = FALSE,
                              saveTOMs = TRUE,
                              saveTOMFileBase = "PBCreponseTOMB",
                              verbose= 3)

I come up with error code saying:
            Number of items to replace is not a multiple of replacement length.

Does anyone please know what this indicts to? I have looked on all the WGCNA tutorial pages but cannot find this error message anywhere.
Thank you kindly

Could be a bug in the code. You will need to provide more detail. When you get the error, type traceback() and post what you get. Also, post the output of sessionInfo().

Error code:     Error in clustDist[ii, ] <- if (intNetworkType == 1) 1 - abs(cr1) else 1 -  :
number of items to replace is not a multiple of replacement length

sessionInfo()    R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.5 LTS
Matrix products: default

3: sizeRestrictedClusterMerge(datExpr, clusters = membership, clusterSizes = clusterSizes, centers = centers, maxSize = preferredSize, networkType = networkType, verbose = verbose, indent = indent) 2: projectiveKMeans(datExpr, preferredSize = maxBlockSize, checkData = FALSE, sizePenaltyPower = blockSizePenaltyPower, nCenters = nPreclusteringCenters,  verbose = verbose - 2, indent = indent + 1)

1: blockwiseModules(B_CC_4, power = 6, TOMType = "signed", minModuleSize = 10,       reassignThreshold = 0, mergeCutHeight = 0.25, numericLabels = TRUE, pamRespectsDendro = FALSE, saveTOMs = TRUE, saveTOMFileBase = "PBCreponseTOMB",
verbose = 3)

I can't add more of the sessionInfo() as too many characters.. can you let me know which parts are key to review. Thank you so much for your kind assistance. I am abit lost with the error and re-run over and over again and re-formatted my raw data, cant see any errors there, thank you

@holly You can [edit] your question to include these details and other that might be needed to give a good answer to your question (And help future readers to find it, as comments are not indexed by searcher)

Please try installing and using this development version of WGCNA: https://www.dropbox.com/s/fuu73li4ajjfr3z/WGCNA_1.67-90.tar.gz?dl=0 and see if it solves the issue. Thanks again for reporting the issue.

Thank you kindly for all your kind help, very much appreciate it

Looks like a bug in the code which I will try to hunt down and fix. In the meantime, I would suggest playing with maxBlockSize argument to blockwiseModules. Try increasing it as much as your available RAM allows (see the paragraph "A second word of caution concerning block size" on page 6 of WGCNA tutorial I, section 2c, for some guidance on how to set maxBlockSize). Hope this helps.

Thank you so much for looking over this for me, I very much appreciate your time and effort. I will try and play around with the maxBlockSize as suggested. Thank you kindly

I also got this error:
Did you manage to solve the issue @PeterLangfelder > > sessionInfo() R version 3.5.1 (2018-07-02) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 7 x64 (build > 7601) Service Pack 1 > > Matrix products: default > > locale: [1] LC_COLLATE=English_New Zealand.1252 LC_CTYPE=English_New > Zealand.1252 LC_MONETARY=English_New Zealand.1252 [4] LC_NUMERIC=C > LC_TIME=English_New Zealand.1252 > > attached base packages: [1] stats graphics grDevices utils > datasets methods base > > other attached packages: [1] WGCNA_1.67 fastcluster_1.1.25 > dynamicTreeCut_1.63-

Please try the newest version (1.68) of WGCNA from CRAN. This bug should be fixed.


text:
Is there a way to assemble contigs starting from a specific sequence?
My work involves searching for marker genes/fragments in metagenomic databases (like the Sequence Read Archive). Once I find these sequences, I would like to know more about the neighboring genomic region.
Is there a way I could assemble only sequences that create a contig which contains my region of interest? Contigs which don't contain this region are not useful to me. My organism of interest might represent a minority of the metagenome, and assembling everything in the dataset would use a lot of computing power.

Good question! I know that some of the assemblers for mtDNA are working on this principle - searching for reads mapping to a conserved mt gene and then extending the sequence, maybe you could tweak one of those tools.

This is one of the primary use cases for which spacegraphcats (preprint and code) was designed. The "neighborhood queries" discussed in the paper sound particularly relevant. I don't have any personal experience with spacegraphcats, but the run guide provides some examples of how to index the complete data set and how to query sequences of interest.

Yes, we'd be happy to chat (I'm one of the spacegraphcats authors). Just post questions over on our github repo!

You might also be interested in citations 21 and 22 in the spacegraphcats paper - metacherchant and a paper that does something similar.


text:
Using average of sequence similarities when delimiting genera - how to deal with outliers?
Sometimes I am trying to see where a (bacterial) genus ends and where another one begins based on 16S rRNA phylogenies. Of course, the length and support of the branch matter a lot, but this question in particular is about sequence similarity.
As far as I understand it, in general, sequences within the same genus are expected to have roughly >94-95% sequence similarity at the 16S rRNA level. But that leaves out a lot of variation. For example, in doodle A below, there may be <90% similarity between sequence 'a' and sequence 'e' ; however, there is always >95% similarity between, say, a and b, b and c, c and d, d and e...So all would belong to the same genus.

Previously I had considered just comparing every sequence with every other sequence in my group of interest, get an average and a standard deviation, and see if the average > 94-95%. However, it occured to me that if I have a tree shaped like that in doodle B (a couple of closely-related sequences from one genus and one very distantly-related sequence from another genus), I could get the very same average (and only slight changes in the standard deviation).
It feels like there should be a simple mathematical property to distinguish between A and B based on similarity values, and/or point out outliers like sequence 'a' in B, but I do not know what that would be. Help would be appreciated!
Mock dataset A:

sequence1
sequence2
sequence similarity (%)

a
b
99

a
c
93

a
d
92

a
e
95

b
c
99

b
d
98

b
e
97

c
d
95

c
e
94

d
e
88

Mock dataset B:

sequence1
sequence2
sequence similarity (%)

a
b
90

a
c
90

a
d
90

a
e
90

b
c
99

b
d
99

b
e
99

c
d
99

c
e
97

d
e
97

Wouldn't the average for all comparisons to "a" in dataset B be below 95? That would indicate that a is not the same genus as the others, following your definition of >94-95 similarity.

@Niklas Thank you. My issue is that this is simple enough for one sequence, but would get more difficult to see from the table if I add more (one hundred? five hundred? a thousand?). Of course, I can still spot the outlier easily enough in the tree...most times...but I feel that adds a degree of subjectivity that I am not quite content with. I would prefer something more rigorous.

@M__ I am not quite sure how to answer. I guess either? Basically, for me A and B correspond to a group of sequences (a-e), that could be seen as taxa, and can be easily represented as trees. So more or less interchangeable. I can visually spot outliers in the trees, but I hoped for something more objective, and working with the taxa at the level of sequence similarity felt like an option. Maybe it's not the only option?

@Laura I am not really sure if that is so complicated, basically you would just need to compute the average % identity across all comparisons that contain one of your taxa and see if that is below 95%? You could also plot it to make it more visual. It depends a lot on what you expect (is this a frequent or rare thing, will most species always be from one single genus?).. I think to "look at the tree" using some metric, total branch lengths (sum of all branches) could be informative for your problem, have you tried using that?

The basic question is where do genera start and stop within phylogenetics and the 16S 5% rule? This is not a trivial question within molecular taxonomy, especially for bacteria.

The key is to define the tree structure for 16S
Define the tree structure for other genes

Background and tree stats The difference between Tree A and Tree B is referred to as 'balance'. Tree A is unbalanced, whilst Tree B is balanced and that is a formal description. I would personally use the term Tree A was "heavily paraphyletic", or "extensive paraphyly" - at a technical level that is accepted, I have used frequently. Equally I'd describe Tree B as a "singular monophyly". The statistical basis for tree balance is called the Colless value.
"Colless statistic". It's fairly easy to calculate. The reason it is not often used is because the tree needs to be phylogenetically robust (below). Not much importance is given to it because if a branch isn't robust thats a headache for this statistic (it can't tell the difference between robust and non-robust).
Difference between branches The statistical difference between taxa A and taxa B is measured by bootstrapping (I'm sure you know this), where >= 80% is bifurcation is robust.
Colless statistic might be used, but its fairly obvious, bootstrapping is needed and iqtree is best package for this operation.

The second more difficult part of the question is the delineation of a genus. The 5% threshold of 16S has always been seen as a guide. I have never seen it rigidly enforced and it creates a conflict with phylogenetic theory as your example points out. Telling what is and ain't a genus ... you are right that a key is the tree structure.

Tree structure and bootstrapping ...

A declared monophyly must be supported by a bootstrap > 80%
A declared paraphyly ... each paraphyletic branch must be supported by a bootstrap > 80%

16S versus other genes and tree structure
The tree structure needs confirming particularly as the bootstrap values 16S will often be low. The starting point is to consider the traditional MLST loci (I forget what they are but stuff like atpD, etc ...), build trees from them and :

Don't combine 16S with protein genes (they evolve at different rates), although it could be done in a partitioned model
MLST protein genes and core-house keeping genes can be combined without partitioning the data (there's a test, but it'll basically work).

Take the MLST genes repeat the tree and assess the bootstraps. Basically it will produce the same tree, but: a) confirms the tree structure; b) can provide stronger bootstrap values (which confirm the tree).

5% stuff
All the above needs doing before any 5% stuff is considered. It is my personal view that using 5% 16S as a mean value is sort of bending the numbers, particularly in a paraphyly where the divergence is e.g. between 2 to 10%.
Monophyly vs paraphyly This where the tree structure is needed. Basically, a robust monophyly would be regarded as a taxonomic singular unit. In this case if its close to the 5% threshold genus rule:

even if the pairwise distance between taxa within the monophyly is say 8% via 16S;
and this is supported by other protein loci;
and there's no suspicion this could represent a family.

In the above criteria I'd say its evidence of a genus. 10% is a bit big distance however.
If the tree structure is paraphyletic and the maximum different is say 10% between two taxa - one at the top of the paraphyly and the other as a most recent child ... then thats tricky. It could easily be two genera within the paraphyletic tree structure. For example, the vertebrates as paraphyletic, viz. amphibian, reptile, bird, edentates, marsupial, mammals and that represents 6 monophyletic Classes. Thus if that was my result I would leave open the possible interpretations: maybe this is a single genus; alternatively it could be two (or more) genera. The key to this is defining the nodes (internal branches) via bootstraps to be confident of the tree structure.
Generally a single robust paraphyletic tree structure, i.e. heavily imbalanced, where each paraphyletic branch is robustly supported would not be regarded as a single genus especially when 16S genetic differences is >> 5%. Taxonomically, this scenario would need further phylogenetic analysis (separate question) and classical bacterial taxonomists would want phenotypes before splitting a genus.

Reviewing the question against the answer - what I would do for the monophyly (diversity 10%) is split off the basal groups (at 10%) view robust bootstrap support. Then forward the in-group as the genus and the full monophyly as a putative genus. So that 10% is not the divergence unequivocally defining a single genus. So thats quite easy on the picture drawn. Thus the in-group is <8% because 10% is a bit high.
It is harder to do this with a paraphyletic tree structure and thats a separate question on the advanced techniques needed to investigate this.

I would like to note out that the terms paraphyletic and monophyletic do not refer to trees, but to clades.

Hello, all terms of the "-phyletic" group inherently describe clades, a clade can either be mono-, para-, or polyphyletic. Since a tree has no inherent clades, these terms only make sense in combination when adding clades (e.g. taxonomy) to the tree. Clades are not necessarily monophyletic (which is why they can also be para- or polyphyletic).

Clades versus monophylys. I'm not getting into that discussion because of its acrimonious history. I have used "paraphyletic clade" and singularly used "clade" as monophyly as a published definition (a paper reviewed by Science journal). The semantics are not rigid and phylogenetics is more concerned with statistical robustness than "one person says this, another says that".

I am only raising this because you refer to the trees as  "heavily paraphyletic", or "extensive paraphyly"

Thanks, personal semantics as mentioned ("imbalance" better)

This is quite advanced and it will take a bit of time to wrap my head around, but it's very informative. I guess I'd hoped that there is something simpler and more objective than what I'd already been doing (looking at the tree shape, bootstrap support etc) but I guess that even when adding extra numbers it all comes back to that...One thing I'm not convinced about is that a robust paraphyletic structure should be split if you can 'step' from one end to the other using branches with >97% similarity to each other...but that is indeed another question. Overall, thank you!

Hi @Laura, glad it helps, you can simply ask the other questions, e.g. 97% as a separate question. because it depends what you mean by 97%, if that is base tip to the "most recent" tip, or distance between the ancestral node and the other tips.


text:
rRNA genes missing from metagenomic bins; is there a way to recruit them?
I have a couple of metagenomic bins which are okay in quality, but often missing rRNAs (16S, 23S...). I assume this is due to high population variability, combined with the high conservation of rRNAs, and reads of at most 300 bp. My binning software of choice is vamb.
Could this be fixed? Once a metagenomic bin is built, is there a way to see which 16S (or multiple 16S possible sequences) are likely to be associated with it? I considered using tools which explore the assembly graph, either manually (Bandage) or automatically (spacegraphcats), but I'm not sure if this is the right approach. After all, if the rRNA genes were linked in the assembly graph to genes that now clearly belong to a bin, they would have been added to that bin...It feels like this is a problem others may have encountered before, so I was wondering if anyone has an idea : )

I am not super familiar with vamb, but it doesn't appear to use the assembly graph: "Vamb is a metagenomic binner which feeds sequence composition information from a contig catalogue and co-abundance information from BAM files into a variational autoencoder and clusters the latent representation." Most binners I'm familiar with don't pay attention to the assembly graph. I'd expect that there is some assembly graph signal, possibly the issue is that there is a really nasty knot making it hard to navigate the graph properly.

@MaximilianPress Sorry, I was confusing. I meant using the assembly graph from my assembly program (MEGAHIT) and then trying to correlate the contigs with the bins. If you'd like to expand on the knot it would still be helpful! ('here are the reasons this would be difficult' is a valid answer : )

Ah, I mostly meant that "After all, if the rRNA genes were linked in the assembly graph to genes that now clearly belong to a bin, they would have been added to that bin..." isn't necessarily true if your binner doesn't pay attention to the graph. In conclusion: yes please do use spacegraphcats or a similar tool.

I considered mentioning this but considered it implausible. However, apparently people are actually able to use PE shotgun data to link rRNA genes to MAGs.
This is something that is quite trivial with linked-read or Hi-C data, but I am assuming that you don't have these data (and don't have the desire to go back and generate them). However you presumably have PE shotgun data that you used to generate your MAGs.

Really impressed that you linked something so recent, thank you!


text:
What are some ways to check metagenomic bin quality?
I am new to metagenomic binning. I've used CheckM in order to estimate completeness and contamination values, and most of my bins of interest appear to have good values. My workflow was pretty straightforward (BBDuk for quality control > Megahit for assembly > vamb for binning). What else could I do / where do I go from here? (especially regarding quality control)

Some papers check 5S, 16S, 23S and ≥18 tRNA genes for completeness. For example, you can see this paper.

There are many other tools. If you want to get nice statistics and graphs I recommend Amber, there is also Anvio for manual refinement of the bins.


text:
Is this a complete bacterial genome? (CheckM)
I've been working with a bacterial genome assembly. The initial stats look pretty good (3 main pieces, depth ~50x [20x - 100x] if I calculated it right). I'd expect this to be a virtually complete genome apart from some low-complexity regions.
However, when I use CheckM to estimate completeness / redundancy, I get ~99% completion, <2% redundancy. (out of 579 expected single-copy marker genes, 563 were present once, and 11 were present twice). I am especially concerned about the missing genes.
I can think of several possibilities:
a) Limits of CheckM. This software uses a sample of 82 genomes to select marker genes for this phylum: however, there are no available genomes in the same genus as the one I'm working with. It might be possible that this is natural variation.
b) The "single-copy" genes present twice might be due to misassembly. But while the sample was not axenic, the organism of interest was majoritary, and long reads were used as well, so this shouldn't be an issue?
c) The genes were somehow missed by sequencing. (but the depth is pretty good?)
I don't know how likely a) is, so it's difficult for me to assess whether what I have is a normal result, or something that requires additional sequencing in order to make sense. I would appreciate your thoughts!

See this table from the CheckM table to put this number in context. In general, your numbers sound excellent, even compared to isolate genomes.

I'd recommend for context taking some finished microbial genomes (E. coli K12 etc) and running CheckM. I bet you'd be surprised what you find.
I think that you can probably not worry about (c) as a possibility. If you are very worried about (b), you could map reads against your assembly, filter out the ones that don't map well, and assemble those to see what you find. Long reads suggest that your genome is in good shape. You can of course do manual curation of your genome using tools like Anvi'o, or by scaffolding/gap-filling your assembly with long reads or an orthogonal strategy like optical mapping or Hi-C.
I think that the influence of (a) is significant. Overall it sounds like you have a reasonable genome.

Don't read into the CheckM result too much. Genes having one copy in some species may have multiple copies in other species and genes present in some species may be absent in others. This may even happen between strains, let alone between species or between genera. It is rare to see 100% checkM completeness and 0% contamination. A common standard for "near complete" assembly is >90% completeness and <5% contamination. Your assembly is far above this standard.
Meanwhile, beware that checkM only evaluates single-copy genes but multi-copy genes are usually harder to assemble. There might be complex multi-copy genes in the three assembly gaps. Without getting a circular contig, it is hard to say for sure your assembly is really complete.


text:
Is is possible to predict ncRNAs from sequence and homology alone?
I'm working with a set of homologous genes (let's call it gene A) from several bacterial species.
I know (from previously published research) that in gene B (a close paralogue of gene A), there is a non-coding, antisense RNA. This antisense RNA regulates gene B by forming a duplex with the gene B mRNA.
I have noticed that, sometimes, an antisense RNA would be predicted for gene A in the same location in the IMG/MER database. This is not predicted for all of the available gene A sequences, only some of them. IMG/MER uses INFERNAL to predict RNAs, I think.
My question is - is there any way to assess how real this is without wet lab data? As it is, it could be that there's no actual antisense RNA in gene A, and this is a spurious prediction due to homology. Or it could be that the antisense RNA is present in all of the gene A sequences, and the RNA prediction software simply missed some of them. Is there a computational way of testing which of these options (if any) is more likely? Thank you for your time!

https://pubmed.ncbi.nlm.nih.gov/17631615/

@Chris_Rands I don't think that software would work for my current problem (even if it wasn't abandonware)... It seems to distinguish between 'coding' and 'non-coding', which it doesn't feel like it'd help when the antisense RNA is inside the gene (so, it is coding in one frame by default). Thank you for trying though!

Is there any conserved secondary structure in this RNA? One could align all known instances of this RNA and see if there is any significant covariation according to R-scape. That could serve as a computational signal that the RNA is real (although wet lab experiments would still be needed).
Another question is whether this RNA is already in Rfam? If so, then the job of making an alignment of known instances is easier.


text:
How to edit the headers of multiple fasta files from multiple folders? (recursively)
This is an expansion of my previous question, How to edit the headers of multiple fasta files from multiple folders?
My directories are organized as follows: one main directory, in which I have multiple project directories. Inside each project directory, there are multiple metagenome directories. Inside each metagenome directory, there is always a single .fna (fasta) file.
Names: project directories start with 'mgp', metagenome directories start with 'mgm' and end with '.3'. They end in numbers, but the numbers are not related. (eg. project directory 'mgp83581' contains metagenome directory 'mgm4729322.3').
What I want is a loop that would go through all these fasta files, and add to the headers codes corresponding to project and metagenome directory. For example, if the fasta file in project directory mgp83581, metagenome subdirectory mgm4729322.3 had the following headers...
>seq1
>seq2
...then I would want to change it to:
>83581_322_seq1
>83581_322_seq2
Basically, I want to add the number code associated with the project, and the last three digits of the 7-digit code associated with the metagenome (apart from the '.3' part of the file name).
Below is the answer that I got for my previous question, where I didn't have project directories, and where the 3-digit code for the metagenomes stayed in a manageable range. (thanks terdon!)
for i in {314..727}; do
    sed -i.bak "s/^>/>${i}_/" mgm4761${i}.3/*.fna
done
Using ubuntu. I would appreciate your help!

You can edit your for loop so as to include the paths to your files of interest.
The command below will give you the relative paths (relative to where you run the command, I think it is better to run this within your main directory) of your files of interest:
ls mpg*/mgm*.3/*.fna

You can use the output of this command (your files of interest) to iterate over via your for loop. Your code will look like:
for i in {...}; do
    sed -i.bak "s/^>/>${i}_/" i
done

Alot slicker than using Perl

Thank you! I've tried to do this. (though tried 'while' instead)
ls mgp*/mgm*.3/*.fna | while read line ;  do  header=$(sed -r "s/^mgp([0-9]*)\/mgm[0-9]{4}([0-9]{1,})\.3\/(mgm.*\.fna)/\1_\2/")   sed -i.bak "s/^>/>${header}_/" $line done
It seems to only edit the first of the files I'm interested in (but it does it correctly as far as I can tell). Any idea of what could've gone wrong?

I am not sure if you really need that complicated command within do, fasta headers always start with >, which seems to be parsed in your sed. My suggestion would be: ls mgp*/mgm*.3/*.fna | while read line; do sed -i.bak "s/^>/>${i}_/" $line; done.

I've tried. For some reason it doesn't work for me. All it does is add an underscore to the fasta headers in the edited files.

It is because I forgot to change the  ${i} part in sed while editing the code regarding your comment above with while read. You can do something like ${line:Ind:Len}, Ind being the zero based index and Len is the length of the substring you would like to extract from the full file path.


text:
Is there any value in scaffolding the output contigs of MEGAHIT assembler given a metagenomic dataset?
As far as I understood, for most assembly programs, the scaffolding step takes into consideration paired-end information in order to get from contigs (contiguous sequences) to scaffolds (longer sequences that might have some N-filled gaps).
My assembly software of choice, MEGAHIT, uses paired-end information to build the contigs, but it does not output a standard scaffold. So I am wondering the following, related things:
-is it meaningful to run a scaffold program on the output of MEGAHIT? I imagine there might be some instances in which paired-end information could span a gap.
-Which software would you recommend for it? (I've tried soapdenovo2 and SSPACE, but they appear not to be actively maintained so I have the issue of 'it doesn't work and I can't do anything about it')
-Could relevant information (eg. two contigs being connected by pair-end information) be recovered by alternative and perhaps more user-friendly means, such as exploring the assembly graph with Bandage?
Thank you for your time!

Update 2: It looks like your approach has actually been suggested here as one way to use the PE information. I guess megahit may not be really using the PE information anyways. I still believe that it's kind of weird but other people do suggest it, so maybe it's worth trying what Torsten suggests.
I think that scaffolding with PE information from the reads used as input to Megahit is somewhat sketchy. You could still try it, but I'd be worried about artifacts, just because de novo assemblers are so heuristic.
However, I think that it is perfectly ok to scaffold using orthogonal data. Here are some examples of orthogonal data:

A reference assembly (weird for metagenomics but still possibly ok).
long reads, e.g. ONT or PacBio (low coverage probably fine).
Proximity ligation data, e.g. Hi-C reads.
Mate-pair reads (also lead to artifacts at a somewhat high rate, but technically valid).
Linked-reads (e.g. 10X).
Probably other methods as well.

The tools used in each case would be somewhat different. For a little more information about how you might use this, here is a recent review of tech.
Full disclosure: I work for a company that sells Hi-C kits for such applications.
Update: Realized that I missed one part of the Q. I think that visually exploring the assembly using e.g. bandage is always a good idea. Quite possibly you can make some scaffolding decisions that way, but to me it sounds somewhat painful to do, especially in a metagenome where there are going to be a lot of multiple-branching collapsed regions.

You can use SOAPdenovo-Fusion to scaffold contigs produced by MEGAHIT as suggested by one of the developers: https://github.com/aquaskyline/SOAPdenovo2


text:
How good does the assembly of an NCBI prokaryotic genome have to be in order to argue gene loss?
NCBI has several labels for assembly completeness - Complete, Scaffold, Chromosome and Contig. Complete would be a circularized genome (or linear, rarely)
For a Complete genome it's fairly straightforward - if not present, a gene can safely assumed to be missing (or possibly present on plasmids, if those aren't available in the database).
But what about the others? With Contig-level genomes, I'd feel uncertain about making any statements (especially since some of these projects are split into thousands of contigs). But what about Scaffold, or Chromosome? Are they expected to be complete in terms of genes, but just missing repetitive sequences?

I would suggest that in the absence of "Complete" no very confident assertion can be made about gene presence/absence in a genome assembly. You need to actually  have the chromosome sequence to make strong statements.
I think that most people (me included) would reluctantly accept that genomes with very good BUSCO or CheckM numbers (depending on your organism's taxonomy) are probably more or less accurate in terms of presence/absence. For metagenomics, I would extend that to the MIMAG criteria, which are pretty stringent.
These criteria are of course independent of the NCBI categories, which are pretty coarse-grained.
That said, I think that you can make stronger arguments about gene presence/absence from the sequencing reads from which the genome assembly is derived. Find all the reads for your genome, map them to the genome assembly and also to the closest related ortholog of the gene in question. How well is the gene covered? How do the alignments look? Harder for metagenomics, but there are still probably tricks one could use.

Thank you! Regarding your last paragraph, are you saying that it is probable that genes would be lost during assembly / binning (therefore not fully trusting a non-complete assembly), but not so much during sequencing (therefore trusting the raw data)? I imagine that might not hold for low sequencing depths, though. I was wondering, how high do you thing the sequencing depth have to be in order for the raw data to be reliable in terms of gene absence?

@Laura yes, I believe that's accurate. Many things are lost/misassembled during assembly that are obviously present in inspecting raw data. You would be looking for a similar coverage to the rest of the genome, unless you expect some funny business in that regard. The issues that arise are in the case of repetitive or highly homologous sequences, which is more of an issue in metagenomics as you might have separately binned closely related or homologous genomes. For example if you have abundant genome A, that will tend to suck up a lot of reads for homologous sequences in genome B.


text:
How to use SRA Toolkit to blast SRA accessions?
I'm trying to look for a gene in a group of SRA files. Web BLAST can't cope with them, so I assume they're too large. I've tried tblastn_vdb from the SRA toolkit, as follows. It had previously worked well for me for WGS files.
module load anaconda3/personal
source activate sratoolkit

cd ~/sra-tools
PATH=$PATH:~/sra-tools/sratoolkit.2.10.5-ubuntu64/bin
tblastn_vdb -db "SRR6256958" -query ~/sra-tools/query/query.txt -out ~/sra-tools/out/results.txt

I get nothing in my output file except for 'TBLASTN 2.6.1+'. My error file goes as follows:
Loading anaconda3/personal
  Loading requirement: fix_unwritable_tmp fix_setxattr
VDB: 2021-07-16T12:03:22 . debug: KSysDirResolvePath_v1 = '/[path]/sra-tools/sratoolkit.2.10.5-ubuntu64/schema'
VDB: 2021-07-16T12:03:22 . debug: KSysDirResolvePath_v1 = '/tmp/pbs.4070476.pbs'
VDB: 2021-07-16T12:03:22 . debug: KSysDirResolvePath_v1 = '/tmp/pbs.4070476.pbs'
BLAST engine error: NCBI C++ Exception:
    T0 "/panfs/pan1/sra-test/vdb3/git/blast_vdb/sratoolkit_2_9_0/c++/src/algo/blast/api/prelim_stage.cpp", line 304: Error: BLAST::ncbi::blast::CBlastPrelimSearch::Run() - search cannot proceed due to errors retrieving sequences from databases

/var/spool/PBS/mom_priv/jobs/4070476.pbs.SC: line 7: 4032551 Aborted [then it lists the original code again]

Does anyone know what I am doing wrong? Help would be much appreciated!

When you say "SRA files", are these fastq? If these files are 'large', BLAST may not be the right tool for this

@Chris_Rands No particular format at the moment since I haven't downloaded them separately. It was supposed to happen automatically? There should be enough space at least for this first file (8.79 gigabases), but the download doesn't even start.

SRA contains raw sequencing reads, blast usually serves in aligning genes / parts of genomes / assembled sequences (not sequencing reads). But I don't know tblastn_vdb. You sure the command is supposed to do what you think it should?


text:
Is there an equivalent of Entrez queries for local BLAST?
I have been having some trouble with running command line BLASTn remotely. Therefore, upon recommendations, and considering I have >100 searches planned, I am considering just downloading the nt database (thankfully possible).
My issue is that one of the settings that proved valuable with my data in WebBlast is running an Entrez query. I get more useful results for each query when instead of one search, I run three, with each of the sub-searches delimited by subject sequence length (e.g. -entrez_query 700:1000[slen] to only get the results which are between 700 and 1000 nucleotides long).
WebBlast and remote Blast have that, local Blast does not.
Since I only really need 2 Entrez queries, I guess it would not be conceptually impossible to get my copy of nt database as FASTA, and then run a script to split it into three directories (700-1000nt, 1000-1400nt, the rest while discarding all<700 nt), then make three blastable databases out of it. However it does feel rather inelegant. I was wondering if anyone might have a better solution!

Comments are not for extended discussion; this conversation has been moved to chat.

I looked back at the question and have a more comprehensive solution.

There is no option for -entrez_query on a local blast search. The information is here

entrez_query   Restrict search with the given Entrez query. Remote searches only.

There is an option without -outfmt for length and that is worth investigating.

length means Alignment length

Perl one liner solution to NCBI's many parameter problem

a) Preprocess your fasta file into separate databases of the required length ... by turning a fasta flatfile db into one where all the subject sequences are 700 - 1000.
perl -0777 -nle 'while (/^(>.+\n^[AGCT-NRYMKWS]{700,1000})\n/gm){print $1,"\n"}' example.fa > example700-1000.fa

b) Now make the database on the restricted length fasta flatfile db,
makeblastdb -in example700-1000.fa -parse_seqids -blastdb_version 5  -out laura700-1000.db -title "Laura 700 - 1000" -dbtype nucl

This db will now contain an equivalent of the -entrez_query 700:1000[slen]. Obviously, you would need to make a new db for each change in the sequence length query. Note the option -taxid_map.
c) Finally, run the blast search on the restricted length blast db.
blastn -query laura.fas -db laura700-1000.db -evalue 1e-6 -num_threads 4 -out blasthits.txt

Code testing
Input is test.fa
>here
AGAGACGTAGAGAAAAAAAAAA
>here2
TAGAGA
>here3
AGAGACGTAGAGA

Modified perl oneliner because sequence length is different,
perl -0777 -nle 'while (/^(>.+\n^[AGC-T]{6,13})\n/gm){print $1,"\n"}' test.fa

Output
>here2
TAGAGA
>here3
AGAGACGTAGAGA

Thus the sequence that was outside the range 6 to 13 (sequence here2 and here3) is excluded (here). This seems alot easier to me that all those blastn options.
Note NRYMKWS refer to the degenerate nucleotide code

"Often, one needs to search ... a specific subset of sequences within an existing database. For these type of searches a convenient way to conduct them is by creating a virtual BLAST database"
From "Use blastdb_aliastool to manage the BLAST databases" https://www.ncbi.nlm.nih.gov/books/NBK569848/


text:
Prokka + table2asn (NCBI submission)
I have a couple of MAGs of interest that I am trying to submit to the NCBI (with annotations).
The process of submitting the annotations is rather challenging to me. I have the output of Prokka, which includes multiple files, but the one I think I need is the .gff (which is GFF3 format).
The NCBI website descibes how to obtain the required file format (.sqn) from .gff files (https://www.ncbi.nlm.nih.gov/genbank/genomes_gff/). I downloaded the required file (table2asn), and ran the following code:
"table2asn" -M n -Z -J -c w -t "template.sbt" -V b -j "[organism=Organism] [strain=strain] [gcode=11]" -i "fasta.fna" -f "annotations.gff" -o output.sqn

However, I got an error:
GFF annotation requires locus tags, which are missing from one or more genes, so the command line argument -locus-tag-prefix is needed. Which was a bit strange, since I could see that the original file had a locus_tag entry. The format liiked kind of like this:
S6C7    Prodigal:002006 CDS 244817  246085  .   -   0   ID=AJNLIECJ_06723;inference=ab initio prediction:Prodigal:002006,similar to AA sequence:ISfinder:ISNpu7;locus_tag=AJNLIECJ_06723;product=IS701 family transposase ISNpu7
S6C7    Prodigal:002006 CDS 246142  246702  .   -   0   ID=AJNLIECJ_06724;inference=ab initio prediction:Prodigal:002006;locus_tag=AJNLIECJ_06724;product=hypothetical protein

Anyway, I added -locus-tag-prefix AJNLIECJ to my code, which at least got me an output, but now the output is a touch weird:
 gene            complement(246142..246702)
                 /locus_tag="AJNLIECJ_006724"
 CDS             complement(246142..246702)
                 /locus_tag="AJNLIECJ_006724"
                 /inference="ab initio prediction:Prodigal:002006"
                 /codon_start=1
                 /transl_table=11
                 /product="hypothetical protein"
                 /protein_id="AJNLIECJ:AJNLIECJ_06724"

So basically, the /protein_id is repeated, and sometimes the numbers between locus and protein ID don't even match (some of the protein ids 'go ahead' for some reason).
My thoughts:

I have no idea how to fix this
Is it even necessary to fix it? Since all NCBI protein IDs tend to be in a consistent format anyway, so no matter what I add, it will get updated.
Should I remove the /protein_id field? How?

Thank you for your time!

Hi @Laura, contact the NCBI Helpdesk. They will respond with the exact answer within 1 day and can be within a few hours. You can copy paste, or just cite this question on the request form.


text:
In VEP annotation, how is the codon field interpreted?
After annotating with VEP a VCF file, we obtain different fields. One of them is called Codons which represents the affected codon in the transcript of the gene. Below is a screenshot of Insertions from a sample:

Variants and INDELs affect genes differently depending on the transcript we are looking at. My question is, for some notations for Insertions like in the screenshot, how can we interpret these results? Let's take the example of -/TC. Is TC inserted at the end of the previous codon? or the next codon of the transcript?
Below are the different representations I have in mind where | is the position of the indel and ^ is the possible insertion scenario of TC.
                     Before insertion               after insertion
                           v|
first scenario          xxx xxx xxx                   xxx TCx xxx

second scenario         xxx xxx xxx                   xxx xxx TCx
                            |  ^

My gut is telling me that this depends more on the strand I am looking at. If the strand is +, it is likely to be the second scenario. Otherwise, it would be the first scenario. Am I thinking correctly? or does this notation have another meaning? Or maybe VEP is not able its position?
Thanks for the help.


text:
How can I get the HGVS notation of a mutation from a FASTA?
In an old article, they gave a table of mutations found in a gene and its positions which points to another article containing the FASTA sequence.
This FASTA is the precursor sequence of the gene. How can I know which transcripts those mutations belong to? And how can I obtain their HGVS format?
For example, C-7111-G should have an HGVS notation like c.XXXC>G.

What's the gene? There are tools you can try, but we'd need to know the gene first. And no matter what you do, this will be a manual process since the old FASTA sequence you have is likely inaccurate anyway.

The gene described in the article is POMC. What I did is the following: first, I copied the same fasta sequence from the old article, then I pinpointed the exact location of the mutations from the table. Second, I copied a couple of unique nucleotides from the left and right of the mutations. Third, I matched those sequences from the actual gene fasta sequence to locate the real nucleotide, and I got its position in the current hg38 build. Then the rest is simple, I searched in the databases what mutations are found in those emplacements.

Yeah, that's probably the best thing you can do. I was going to shamelessly plug the product of the company I work for, https://varsome.com, which is free and can help you translate variant notations, but if all you have is the position relative to an old fasta sequence and don't know the actual transcript, we can't really help much more. You'll be able to get the proper HGVS notation easily though. Same for https://variantvalidator.org/ with which I am not affiliated and which is a great tool for validating HGVS.

Thank you @terdon. I will check those websites. Cheers!


text:
How to get the splicing regions (BED) of all genes for Human GRCh38?
How can I obtain a BED file that describes the splicing sites of all genes for GRCh38 human genome ?
Thank you

Interested to know if there are special file format for RNAsplicing.

Is it the alternative splicing events you are looking for? https://bioconductor.org/packages/devel/data/annotation/manuals/alternativeSplicingEvents.hg38/man/alternativeSplicingEvents.hg38.pdf

@PaulEndymion thanks for the link. Does the alternative splicing events mean all splicing regions, even for the same gene ? if that's the case, then yes, I am looking for something like that !


text:
Are there databases to annotate non-coding mutations?
I would like to ask if there are some known databases that allows to annotate non-coding mutations in a given region for the human genome?
Preferably, if the database is downloadable so I am able to automatize the processes of annotation.

What do you mean with non-coding mutations? SNPs? Check dbSNP.

What sort of annotations are you looking to include? Ensembl's VEP is great for annotation, and can include many different annotations applicable to both coding and non coding regions, e.g. minor allele frequency, conservation, pathogenicity predictions through CADD (via a plugin) and many other things

@benn yes SNPs, I know about dbSNP, I wanted to know if there are other known databases !

@JenG I have VEP installed on my cluster ! but CADD also comes with dbSNP plugin, but thanks I will have a look on it

We could give you better answers if you explained what sort of information you are looking for. What kind of "annotating" do you need?

Agreed. dbSNP can easily be downloaded and integrated into an automated annotation pipeline. Please provide more details about what precisely you're looking for.

Sorry for the late reply. I am looking to annotate non-coding mutations (introns, UTRs, etc... but also exons if possible) coming into active enhancers of the human genome of the pancreatic islets. For example, if a mutation comes into an active enhancer in an intron, what kind of predicted impact would that mutation do ?

VarSome allows you to annotate any variant, including non coding ones. While ACMG classification will not be available if the variant isn't in a gene's transcript (the ACMG criteria are only applicable to variants in transcribed regions), non-coding variants in transcripts will be fully annotated. For example, if you look up variant chr17-41197559-G-T which is a non-coding variant found in the UTR of one of the transcripts of BRCA1, you will get a lot of information (probably best to click on the link, the image has a lot of data):

For non-coding variants that don't fall in any transcript, you can still get frequency data, when available, as well as DANN pathogenicity predictions, GERP score, an overview of the genomic regions and other variants and/or srtuctural variants that can be found around the region, and any community annotations that may have been submitted by other users.
There is also an API you can use to automate the annotation, and VarSome Pro that lets you annotate an entire VCF file, but those are both paid services.

Disclaimer: I work for the company behind VarSome.

This looks interesting but yikes I don't think it would be much of use in my case since automating the annotation is in my needs. Thanks though.

@Law yeah, the API and/or the Pro version where you can send VCFs would probably be helpful for you, but those aren't free. I thought I'd mention them in case you do have a budget.


text:
What does the number mean in an HGVSp annotation?
Let's take the example p.Arg452Pro that I got from an annotated VCF file that contains an HGVSp column.
What does 452 mean in this case ?
Theoretically, is it possible to have the same HGVSp on the same gene but on a different transcript level ?
For example:
Chromosome      Position        Ensembl_Gene        Ensembl_transcript      HGVSp
4               4613173         ENSG00052342        ENST00141456            p.Arg452Pro
4               4613176         ENSG00052342        ENST00141458            p.Arg452Pro

Thank you in advance.

the 452 amino acid of the protein?

It's the position in the protein. p.Arg452Pro means that there is a variant that changes an Arginine residue at position 452 of the protein into a Proline.
And yes, it is absolutely possible to have the same HGVS p. notation on the same gene but on different transcripts. More than possible, it is quite common. Your particular example isn't actually a real variant since that position on chromosome 4 is intronic in both hg19 and hg38, but TP53:p.R175L1 is a real example of what you describe:

As you can see in the image or by visiting the link, this variant on TP53 has the same HGVS p. notation on 6 transcripts. It also has a different p. notation on two more transcripts that are not shown in the image as well as no p. notation for a few other transcripts in which it doesn't affect the protein product.

1Disclaimer: The link goes to https://varsome.com, a free human variant search engine that is developed by the company I work for.

I actually gave a random position lol... but thanks for the explanation. Then, I have another question, what does the protein position 452 mean? Is it the number of the protein for that gene on the exons for that specific transcript ? it's not the genomic position, right ? For example, if we have a different transcript, with the position 453, can I say that the previous amino acid is 452 Arginine ?

@Law it's the position in the protein that is produced when the transcript is translated. The number refers to the protein sequence, not the genomic sequence and not the transcript sequence. The 1st amino acid is at position 1, the second at position 2 etc. Your putative variant would be at position 452, it would affect the four hundredth and fifty second amino acid. Whether that is the same for another transcript will depend on what that transcript codes for.

Do you know if clinical diagnostics are transcript-dependent, or would the same protein change on different transcripts be considered equivalent for clinical purposes?

@BradyGilg that depends entirely on the transcript and protein in question. Do both proteins have the same function? Are they expressed in the same tissue? This isn't a simple question. You might want to post a new question instead.

What terdon has answered is right. For more details you can visit the official HGVS recommendations here


text:
Is it possible to obtain RefSeq transcripts from Ensembl database for GRCh37 to match Ensembl Genes?
I know that with GRCh38, it is possible to extract the matching RefSeq transcripts for Ensembl genes/transcripts/proteins using Biomart from Ensembl.
However, for GRCh37, this option does not exists. Using the following MySQL command returns Ensembl transcript IDs for matching RefSeq transcripts:
mysql -u anonymous -h ensembldb.ensembl.org -A -P 3337 -D homo_sapiens_core_99_37 -e "SELECT transcript.stable_id, xref.display_label FROM transcript, object_xref, xref,external_db WHERE transcript.transcript_id = object_xref.ensembl_id AND object_xref.ensembl_object_type = 'Transcript' AND object_xref.xref_id = xref.xref_id AND xref.external_db_id = external_db.external_db_id AND external_db.db_name = 'RefSeq_mRNA';"

Also, I was trying to play around with the databases schema but I am not able to find any solution to problem.
Ultimately, I wish to obtain a table like the following:
Gene stable ID      Transcript stable ID    Protein stable ID   RefSeq match transcript
ENSG00000120329     ENST00000239451         ENSP00000239451     NM_031947.4
ENSG00000164219     ENST00000419445         ENSP00000404676     NM_005023.4
ENSG00000150594     ENST00000280155         ENSP00000280155     NM_000681.4
ENSG00000113578     ENST00000337706         ENSP00000338548     NM_000800.5

Thank you in advance.

The RefSeq match option in BioMart is from the Matched Annotation from NCBI and EBI (MANE) collaboration between RefSeq and Ensembl. It has only been calculated for the up-to-date gene annotation on GRCh38 so cannot be obtained on GRCh37. You can get mapping from Ensembl to RefSeq transcripts through BioMart as RefSeq mRNA ID (refseq_mrna in R) but this is not a perfect match like the MANE, it is a mapping based on sequence similarity and similar genomic location, and there can be mismatches between them.

Thank you @Emily_Ensembl for your answer. I will continue discussing this topic here. Is it possible to obtain also the RefSeq mRNA ID with the versions using biomaRt ?

No, we don't have the versions on the old verion of the database.

Try the GRCh37 version of BioMart, and the External reference "RefSeq mRNA ID" / "RefSeq ncRNA ID":
https://grch37.ensembl.org/biomart
From the Perl code provided by BioMart, it looks like the field name should be "refseq_mrna" or "refseq_ncrna". I'm not sure how that translates to the SQL code.
I feel compelled to add a caution to this, though. GRCh37 is a very old build, having been replaced by GRCh38 in December 2013. If possible, please find a way to use GRCh38 for your analysis instead, and don't create any new software that depends on GRCh37. For example, gnomAD v3 supports GRCh38 (see bioconductor, or here).

Thank you for your comment. Unfortunately, I have to use GRCh37 since gnomAD were annotated with this old version using VEP from Ensembl. So, doing a lift-over or annotating again would lose any supports that the paper had at that time.


text:
What is a signal 'chunk' in the context of Nanopore sequencing?
Does one chunk correspond to the raw signal of a single read read? Why does messing with this parameter affect basecalling results?

Hi Nick, sounds like a valid question, maybe it would be helpful to link a page that talks about this parameter. Also, there are more basecallers for nanopore, no? Which one you use?

If you have access to the nanopore community website you can read more about this here: https://community.nanoporetech.com/posts/not-chunky-enough (I'm not sure to which extent I can copy paste from there so I won't).
A TLDR version: basecalling doesn't happen simultaneously for an entire read, but it is broken in segments which are processed separately and combined afterward. Those segments are chunks. Basecalling results are affected mostly by low-complexity regions, in which a heavily biased region (in terms of nucleotide content) leads to artifacts in the per-chunk basecalling. Those issues are solved by taking a chunking-size larger than most low-complexity regions.


text:
Get residue position from Pymol Selection
I can make a selection in Pymol using the resi parameter. Is there a way to get this resi data back out of a selection object? The only thing I found so far is how to get the position, as in the x/y/z coordinates from a selection.
For context, I have a table with a set of data points for each residue in the amino acid sequence of the displayed structure and I want to write a function where I can print the data based on the current selection. So my thought would be something like the following, where the content of the 'get_residue_positions_for_selection' function would be what I'm asking about.
def show_selection_position_data():
    residue_data = {
        "header": ['variation', 'rmsf']
        , "1": ['high', 0.1]
        , ...
    }
    selections = cmd.get_names('selection')

    if sele not in selections:
        print("No selection")
        return

    residues = get_residue_positions_for_selection("sele")

    for residue in residue_selections:
        print(residue_data[residue])

    return

This certainly a very interesting question and a cool answer. The advice I was given (to pretty much the same problem) was to use Pymol scripting directly rather than via Python.

Ok, after some more searching I found the answer.
You have to use the cmd.iterate function to iterate over the default selection ("sele", first parameter). With the second parameter you can push the variables you need into an array for later access. In this case I wanted the residue number so I used the variable 'resv'. A full list of the variables exposed by the iterate function can be found on the pymolwiki
Note that you can use a function call as the second parameter and if you need a more complex selection than 'sele' you can use selector.process in the iterate function.
Here is my final working function. Note that I used a set for storing the residue ids, since the iterate function technically iterates over the atoms and so the same residue will likely appear multiple times in your selection.
def show_selection_position_data():
    residue_data = {
            "header": ['variation', 'rmsf']
            , 1: ['high', 0.1]
            , 2: ['high', 0.2]
            , 3: ['low', 0.1]
            , 4: ['high', 0.05]
            , 5: ['middle', 0.2]
            , 6: ['high', 0.1]
            , 7: ['high', 0.1]
        }
    selections = cmd.get_names('selections')

    if "sele" not in selections:
        print("No selection")
        return

    userSelection = "sele"

    stored.residues = set()
    cmd.iterate(selector.process(userSelection), 'stored.residues.add(resv)')

    for residue in stored.residues:
        if residue in residue_data:
            print(residue_data[residue])
        else:
            print("{}: no data".format(residue))

    return

Nice. And thank you for taking the time to post an answer!


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

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

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

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

which resulted in the output
From    To
UniRef100_Q8RKD8        null

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

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

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

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

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

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

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

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

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

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


text:
Timeout when downloading the ncbi nr blast database
I am experiencing timeout problems when downloading the NCBI nr preformatted blast database using the update_blastdb script (version 504861).
I run the script with the following paramters
update_blastdb --decompress --passive --verbose nr

and I get the following error message (in verbose mode)
Downloading nr (45 volumes) ...
Downloading nr.00.tar.gz...Net::FTP=GLOB(0x5610fb59b8f8)>>> PASV
Net::FTP=GLOB(0x5610fb59b8f8)<<< 227 Entering Passive Mode (165,112,9,229,195,144).
Net::FTP=GLOB(0x5610fb59b8f8)>>> RETR nr.00.tar.gz
Net::FTP=GLOB(0x5610fb59b8f8)<<< 150 Opening BINARY mode data connection for nr.00.tar.gz (18745730730 bytes)
Net::FTP: Net::Cmd::getline(): timeout at /usr/share/perl/5.26/Net/FTP/dataconn.pm line 82.
Unable to close datastream at /usr/bin/update_blastdb line 202.
Net::FTP=GLOB(0x5610fb59b8f8)>>> PASV
Net::FTP: Net::Cmd::getline(): unexpected EOF on command channel: Connection reset by peer at /usr/bin/update_blastdb line 203.
Failed to download nr.00.tar.gz.md5!
Net::FTP: Net::Cmd::_is_closed(): unexpected EOF on command channel: Connection reset by peer at /usr/bin/update_blastdb line 101.
Net::FTP: Net::Cmd::_is_closed(): unexpected EOF on command channel: Connection reset by peer at /usr/bin/update_blastdb line 101.

The timeout happens after ~35 minutes and a file that is approximately 18GB big is being downloaded, which matches the expected filesize. The checksum file (nr.00.tar.gz.md5) is not downloaded. So I'm not sure which of the two files is actually the problem.
I tested downloading the nt database and everything seems to work fine, so I don't think the script is the problem. For comparison, this is the nt download output for the first file
Downloading nt.00.tar.gz...Net::FTP=GLOB(0x562575a73168)>>> PASV
Net::FTP=GLOB(0x562575a73168)<<< 227 Entering Passive Mode (165,112,9,229,196,51).
Net::FTP=GLOB(0x562575a73168)>>> RETR nt.00.tar.gz
Net::FTP=GLOB(0x562575a73168)<<< 150 Opening BINARY mode data connection for nt.00.tar.gz (4065912989 bytes)
Net::FTP=GLOB(0x562575a73168)<<< 226 Transfer complete
Net::FTP=GLOB(0x562575a73168)>>> PASV
Net::FTP=GLOB(0x562575a73168)<<< 227 Entering Passive Mode (165,112,9,229,195,107).
Net::FTP=GLOB(0x562575a73168)>>> RETR nt.00.tar.gz.md5
Net::FTP=GLOB(0x562575a73168)<<< 150 Opening BINARY mode data connection for nt.00.tar.gz.md5 (47 bytes)
Net::FTP=GLOB(0x562575a73168)<<< 226 Transfer complete

Any help would be appreciated.

The answer of the support stated that it is likely that the size of the first download is the problem and that I should try using another tool like rsync. Because I'm unfamiliar with rsync I decided to write a python script that did the job.
import urllib, tarfile, json

base_url = "https://ftp.ncbi.nlm.nih.gov/blast/db"

# download manifest file to get the filenames
manifest_file = 'blastdb-manifest.json'
urllib.request.urlretrieve(f"{base_url}/{manifest_file}", manifest_file)

with open(manifest_file) as f:
    manifest_data = json.load(f)

# download everything
for file in manifest_data['nr']['files']:
    # download checksum file
    checksum_file = f"{file}.md5"
    urllib.request.urlretrieve(f"{base_url}/{checksum_file}", checksum_file)

    # download the archive file
    urllib.request.urlretrieve(f"{base_url}/{file}", file)

    # check that the checksums match
    calculated_checksum = get_md5_for_file(file)

    with open(checksum_file) as f:
        for line in f:
            line = line.strip()
            if line.endswith(file):
                checksum = re.split("\s+", line)[0]
                if checksum == calculated_checksum:
                    checksum_matches = True
                else:
                    raise Exception(f"Checksum doesn't match expected for file {file}")

    # unpack the archive
    tar = tarfile.open(archive_file)
    tar.extractall(file)
    tar.close()

The function to calculate the md5 checksum is the following (borrowed from this question)
def get_md5_for_file(file):
    hash_md5 = hashlib.md5()
    with open(file, 'rb') as f:
        for chunk in iter(lambda: f.read(4096), b""):
            hash_md5.update(chunk)

    return hash_md5.hexdigest()

Note that this is not the exact script I use, so I can't guarantee it works exactly like this, but I tried to include all the important stuff.
In my version I include checks to whether my existing blast database needs to be updated (I keep the manifest files and compare the manifest['nr']['last_updated'] values) and I also try to skip already unpacked archives while iterating (check if the nr.xx.phd file exists) to save time in case the script fails due to problems in the download.

Same here. I used axel to download it a week ago, but fail to get finished. Also, the downloaded taxa zipped file (60Mb) doesn't match md5sum. Apparently, they update the FTP files every 30 minutes or 1 hour.

Ok, good to know. I wrote them an email and will update when I get an answer

They're usually pretty responsive. Have you heard back yet?

Sorry for the late reply, I posted an answer.


text:
Find similar proteins with similar tertiary structure?
I'm a mathematician studying knotted curves and am trying to find some applications of my work in proteins. Using https://www.rcsb.org/, how do I find proteins that have similar  (but perhaps not the same) tertiary structure? Preferably they would have similar primary structure and function...
Perhaps what I'm looking for is "similar proteins", if there is such a thing?
Also, is there a way in the database to select more "common proteins" (important or well studied, not just some rare exotic one (if there is such a thing)?

As a general question on the board, more detail is often good, because of the answers can be inaccurate due to misunderstanding of the biological context.
There are several ways to do this and at least 3 that I'm aware of.

Blast - sequence similarity searching the PDB database
Contact map identity, I use this approach, but its complicated - and its precision or utility is context specific, i.e. it could be loads of work and miss the vital detail you need.
Direct mapping one protein structure against another, this in particular is true for specific domains, or sites within the protein and the result is angstroms. Its not difficult and I assume this can be scaled to the whole protein.

Point 1 "Blackboxes"
By far the easiest is option one and I suspect I-Tasser (here) use this principle before instigating "threading". This package is considered "good" (others appear to consider it the best of its type), I believe it has done well in CASP tests (Critical Assessment of Techniques for Protein Structure Prediction). So the simple answer is dump your sequence in there an watch the output. Pyre is something similar. These are interesting tools but certainly not definitive, really there a starting point.
Point 1 "Formal Blasting PDB"
This is the easiest and most mature approach for your question, without knowing the details. I've run a test to ask the following question.
Which protein structure available on PBD is most similar to Japanese encephalitis virus (JEV) based on protein homology?

Step A
I've no specific interest in the virus BTW. The most robust Blast engine IMO is NCBI here . I dump a JEV sequence into the box, and click three important options (otherwise it doesn't work).

Database = protein databank (PDB), i.e. structural data
Organism "Japanese encephalitis virus (taxid:11072)": the taxid will come up automatically (if not just proceed as normal)
Next to the Organism box I click "Exclude", otherwise this whole approach doesn't work
Obviously I select "blastp", otherwise it definitely will not work

The output is/was here (but is now inactive). Summarising what it presented, E is so high its "0" (>2e-139) protein ID "77.59%" and just states "envelope protein" Chain A. I hit the link and the protein that comes up has a name which ends in "6A0P_A". This is the PDB reference and Percentage ID and PDB reference are the only two bits of info you need. The PDB reference needs the "_A" removing, thus its 6A0P. The _A just refers to chain A of a dimer.
Step B
You then go back to the Web site you mentioned and input the PDB ID you stated into the search bar. This results in this quite interesting page here. The structure is Usutu virus, which surprised me.
Point 2 . This is much harder and not only involves specific Python code, but a post unsupervised learning analysis and results in an average similarity of structure over the whole protein. There are a number of published algorithms for generating a contact map, or its easy to write your own. Its more advanced, its complicated and unless you've seen loads of them they not immediately understandable. Contact maps angle at a very specialist approach of predicting and understanding protein structure. It is not for the faint hearted and some biological details will be omitted, which may or may not be important for your purpose.
Point 3 . I don't have research grade expertise here, so I forget the exact button to click, but its quite easy via PyMol and other graphical viewers to perform spatial similarity for specific epitopes and this is something you need to learn within your sphere. Basically you will load two proteins into the image viewer and then they will be a "similarity aligning" function and you zoom in on your output to generate an answer in angstroms. You will need to learn this approach for the later parts of your project.
Summary, The orthodox way would be to perform point 1, select several targets and then examine them in detail via point 3. I definitely recommend direct Blast if you go down this route, because it is the most tractable and avoids potential unknown subsampling.  I personally would not recommend Point 2 in a general context because you have to consider an unsupervised learning approach towards processing a heap of contact maps. Ultimately, which ever approach you take you can never avoid Point 3 (which someone else does, in my case).
Note This may or may not be obvious to your but just because two sequence have the closest homology doesn't mean they will have the closest structural similarity.

An explicit model of a "knotted curve" would be very interesting by the way because it would enable an analysis of hinge regions, which are quite important in viruses. I can immediately think of a entire application thereafter via an algebraic solution based around modelling complex non-linear lines.


text:
How many entries in PDB have SS bonds?
Is there a way to figure out what percentage of entries in PDB have SS bonds?
In addition, what other types of bonds in PDB may be worth looking into (not H bonds)?

The easiest way is to check annotations in mmCIF or PDB files.
Looking at my page with statistics of mmCIF tags, last updated in March:
_struct_conn_type.id        59137×metalc 42706×covale 36195×disulf 10725×hydrog
_struct_conn.conn_type_id   5350402×hydrog 1733277×metalc 559092×covale 226472×disulf

Which means 226472 disulfide bonds in 36195 PDB entries at the time when the PDB had 161470 entries.

Great, thanks! Exactly what I needed.
What about ionic interaction (salt bridges), do you have the numbers?

No, I'm afraid I can't help with salt bridges.


text:
Building cancer gene database
I am trying to build a database of different cancer genes i.e lung cancer, throat cancer etc. Searched through NCBI website by searching for nucleotide sequences, but I am unable to find appropriate results. Any help in this regard is much appreciated

COSMIC has a set of mutations and the Cancer Gene Census already curates cancer-related genes. What is the purpose behind your venture?

@RamRS To perform clustering

Why not just cluster based on existing annotations? You could get information from the CGC, use the HGNC symbols or other identifiers to download additional data using biomaRt and then perform all sorts of analysis on that. Why build a new database?

@RamRS i have taken up this as part of my PhD , wherein I will be clustering cancer gene data treating them as pure categorical data and ultimately try to identify similar patterns across clusters if they exist

What will you be clustering (the term data is too generic and vague) and what sort of patterns will you be looking for? If this is part of your PhD, you definitely would have looked into literature on such protocols and the sorts of data available at the gene level on various cancer types. Can you give us any more detail on it?

@RamRS gene expression to be specific

And will you be clustering gene expression data from a public resource such as TCGA? I ask because if it's based on public datasets, this would not be a novel endeavor at all. In fact, there exist explorer tools that let you run correlations between various omics dimensions on the fly.


text:
Repeatscout installation problem
I've installed Repeatscout from https://bix.ucsd.edu/repeatscout/.
I've also installed TRF and nseg required by the program.
I followed the steps:
build_lmer_table -l 16 -sequence ~/file_fasta/od_1000.fa -freq od_1000.freq

 RepeatScout -sequence od_1000.fa -output od_1000_repeats.fa -freq od_1000.freq -l 16

and the I put this command:
cat od_1000_repeats.fa| filter-stage-1.prl >od_1000_repeats_filt1.fa

but this message appears:
bash: filter-stage-1.prl: command not found...

Why? It is all correctly installed, I can't understand the problem

If you're in the RepeatScout directory you need to execute the Perl script with a ./ prefix, like so.
cat od_1000_repeats.fa | ./filter-stage-1.prl > od_1000_repeats_filt1.fa

If you're not in the RepeatScout directory, you need to provide the full relative or absolute path when executing the script. Alternatively, you can put the RepeatScout directory in your $PATH environmental variable.
These are all standard UNIX/Linux concepts that many bioinformaticians take for granted, for better and for worse. If you're lost in all the terminology, a crash course in UNIX basics is probably where you would want to start.

I installed repeatscout on my PC, I created the Repeatscout directory and It worked without errors. The problem is that it can't execute a file larger more than 800 mb in a 16 gb ram of my pc. So I had to use the cluster but I don't know where is the Repeatscout directory.(If I put "locate Repeatscout" I don't find it). I would to find the path to check if all the programs are correctly installed..

It sounds like you need to get some help from your system administrator.


text:
blast command line
I would to use blastx from command line to analize a library produced with RepeatScout in format .fasta to find new transposable elements.
I use this command mylaptop$ tblastx OB150DEC_2_repeats.fa
but it reports an error:
USAGE
  tblastx [-h] [-help] [-import_search_strategy filename]
    [-export_search_strategy filename] [-db database_name]
    [-dbsize num_letters] [-gilist filename] [-seqidlist filename]
    [-negative_gilist filename] [-entrez_query entrez_query]
    [-db_soft_mask filtering_algorithm] [-db_hard_mask filtering_algorithm]
    [-subject subject_input_file] [-subject_loc range] [-query input_file]
    [-out output_file] [-evalue evalue] [-word_size int_value]
    [-qcov_hsp_perc float_value] [-max_hsps int_value]
    [-xdrop_ungap float_value] [-searchsp int_value] [-sum_stats bool_value]
    [-max_intron_length length] [-seg SEG_options]
    [-soft_masking soft_masking] [-matrix matrix_name]
    [-threshold float_value] [-culling_limit int_value]
    [-best_hit_overhang float_value] [-best_hit_score_edge float_value]
    [-window_size int_value] [-lcase_masking] [-query_loc range]
    [-strand strand] [-parse_deflines] [-query_gencode int_value]
    [-db_gencode int_value] [-outfmt format] [-show_gis]
    [-num_descriptions int_value] [-num_alignments int_value]
    [-line_length line_length] [-html] [-max_target_seqs num_sequences]
    [-num_threads int_value] [-remote] [-version]

DESCRIPTION
   Translated Query-Translated Subject BLAST 2.5.0+

Use '-help' to print detailed descriptions of command line arguments
========================================================================

Error: Too many positional arguments (1), the offending value: OB150DEC_2_repeats.fa
Error:  (CArgException::eSynopsis) Too many positional arguments (1), the offending value: OB150DEC_2_repeats.fa

I don't understand what does it mean too many positional arguments

A positional argument is one that doesn't follow something like -foo or --bar. That is, the variable that it should internally map to depends on it's position in the command. You're looking for something like tblastx -db something -query OB150DEC_2_repeats.fa

ok. I have to download a database to use?

Or point it to a remote one.

I've never used blast from command line before. Could you explain to me how to point it from remote one? I'm working with a cephalopoda genome

Just use nt and I think it will figure it out. We only use local databases, so if that doesn't work you'll need to read the documentation.

ok. I've tried but the message is BLAST Database error: No alias or index file found for nucleotide database [nt] in search path [/mydirectory/OB150DEC::]

Use the -remote option https://www.ncbi.nlm.nih.gov/books/NBK279680/

Thank you for your help, but it doesn't work. It reminds a lot of errors. I'll try downloading an nt database


text:
HVG and variation coefficient
I'm looking for highly variable genes from a data set. I have calculated the coefficient of variation of the expressed genes of my dataset and inserted them in a histogram graph. Now I should select only the highly variable genes. Should I select genes with a coefficient of variation >1, or how?
thank you for your help,
regards

Heuristically, I agree with @Jonathan Moore's answer. However, you could try to take a more disciplined approach and model the distribution as a mixture of two distributions, and infer the parameters of the two distributions from that, which will allow you to make somewhat more informed decisions. A basic intro to mixture modeling in R is here.
However, a worry of mine when I see data such as this is that you have some confounder that leads to that spike at CV ~ 1. I would actually plot the mean against the SD of the genes in question, and make sure that it is not due to observations of e.g. one transcript across samples when in most samples there are zero transcripts (see last section here). This would almost certainly be noise, and then what you would be enriching for is genes with very low expression, leading to a very small denominator for the CV calculation. If you are pseudocounting the expression (for example adding 1 to all values to avoid dividing by zero), I would try a few different values for the pseudocount to ensure that the value of the pseudocount is not significantly affecting your CV distribution that you show here. That would be a major warning sign of artifacts.
Therefore I would also suggest setting some sort of cutoff for minimal expression of genes in that distribution. Maybe an average of 1 transcript per observation (sample [bulk RNA-seq] or cell [scRNA-seq]) might be a reasonable cutoff for genes to consider.
I would do this data cleaning and filtering before making any decisions about this CV distribution.

thank you very much for the suggestion. i have already made a cutoff, but probably it was not enough

There seems to be a dividing line around CoV=0.9.  Perhaps this is a good place to start for a cutoff.