text:
Seeking Guidance on Designing DNA Sequences with Specific Criteria
Dear Bioinformatics Community,
I hope you're all doing well. I am currently working on a project that requires designing some 10 nucleotide single-stranded DNAs, (preferably), with specific criteria.
The criteria for the DNA sequences are as follows:

An equal mix of nucleotides (A, C, G, and T)
Thermodynamically stable sequences
Sequences that do not self-hybridize or hybridize with each other

I have been exploring the NUPACK software suite to achieve these goals, but I am finding it challenging to use the tools effectively for this specific task. I understand that NUPACK is primarily designed for generating sequences based on target secondary structures, which may not directly align with my requirements.
I would greatly appreciate any advice, suggestions, or alternative methods you can provide to help me design these DNA sequences. Are there any other tools or resources that might be better suited for this purpose?
Thank you in advance for your insights and expertise. Your guidance will be invaluable as I work to complete this project.
Best regards,

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

Welcome @Hashem and thank you for your question. It would useful to understand the downstream application, i.e. why do you need the synthetic oligonucleotides and some understanding of the organism, possibly in general terms, e.g. just kingdom level and possible diversity. This appears to be primer design, for which there are a lot packages, within that you appear to be wanting nucleotide degeneracy.  The question is not hard to answer - but is not very clear.

I edited your question to make it clearer and to add relevant tags, please update if you do not consider my edits to be accurate.

It is not difficult to simulate every 10mer DNA oligo, here is some Python code that will do that:
(there is probably a more elegant way to do this)
from itertools import product

nts = ["A", "T", "C", "G"]
# create a tuple with 10 copies of nts for the cartesian product
ten_tup = (nts,)*10

decamers = []
# create and iterate over a stringified cartesian product,
# which consists of all 10mers.
for a,b,c,d,e,f,g,h,i,j in product(*ten_tup):
  decamers.append(
    a+b+c+d+e+f+g+h+i+j
  )

At this point, you have some choices to make:

how many 10mers do you want?
what exactly is your criterion for GC content?
what exactly is your criterion for thermodynamic interactions?

I'd suggest either primer3 or viennarna for evaluating thermodynamics. You can run each from the command line or from a web server, but they do have Python APIs, which I use below. ViennaRNA is definitely more trustworthy for sequences as short as 10, so I'll use it here. It will be a lot slower though as it's not cheating by using the nearest-neighbor approximations.
I have made up some numbers for the thresholds below, replace with your own.
from random import shuffle
GC_MIN = 4
GC_MAX = 6
DDG = -5
NUM_WANTED = 10

def test_gc(seq, gc_min=GC_MIN, gc_max=GC_MAX):
  gc_count = seq.count("G") + seq.count("C")
  if gc_count >= gc_min and gc_count <= gc_max:
    return True
  else:
    return False

def test_thermo(seq, other_seqs, ddg=DDG):
  intra = RNA.fold_compound(seq).mfe()[1]
  print(intra)
  if other_seqs:
    # find the most favorable interaction with any other sequence
    inter = min(
      [RNA.fold_compound(seq+"&"+other_seq).mfe_dimer()[1] for other_seq in other_seqs]
    )
  else:
    inter = 0
  homo = RNA.fold_compound(seq+"&"+seq).mfe_dimer()[1]
  if intra > ddg and inter > ddg and homo > ddg:
    return True
  else:
    print("thermo reject!!")
    return False

chosen_seqs = []
# shuffle the decamers to ensure that you get dissimilar sequences
shuffle(decamers)
i = 0
for seq in decamers:
  i+=1
  print(i)
  if test_gc(seq) and test_thermo(seq, chosen_seqs):
    chosen_seqs.append(seq)
  if len(chosen_seqs) >= NUM_WANTED:
    break

# chosen_seqs has your decamers that pass the threshold

In principle, you could even combine these 2 steps but I think it's cleaner to keep them apart unless storing 4^10 ~= 1 million decamers in memory is an issue, which it shouldn't be with modern computers.
Caveat emptor, this code may contain bugs and I haven't tested it other than making sure it runs.


text:
How can I calculate the coverage of a gene in exome sequencing / calculate coverage per gene?
I have a bam file of whole exome human data. I now want to check coverage of a list of genes in my data. Basically, I want to get output like, for example:
Gene    Percentage of coding region covered
A1BG               100 %
A1BG-AS1           100%
A1CF                99.94%
NEK8                100%
FANCI              100%
A2ML1             100%

I want to find out at which depth how much percentage of given genes are covered. for example AT 30X NEK8 gene covered 100%.
Interested in the canonical transcript and if canonical transcript is not present then we can consider other transcript in that case.

Please [edit] your question and add more detail. How do you define "coverage"? Is one read enough? Do you want a minimum coverage? What transcripts of the gene? All of them? Only the canonical? What if there is no canonical? How do you define a gene? Are you also interested in UTRs? How about non-exonic regulatory regions which will probably not be included in WES data?

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

I have edited the question and provide the information hope now its clear

Thanks. So, what data do you have? How can we know the locations of these genes? How do we know what transcript to use? Do you have a bed file with the regions of interest? Also, again, how do we define coverage? Do you want the average coverage over the length of the gene? The maximum? The minimum? Do you want a report by exon? By the entire gene length? What is your minimum coverage threshold to consider something as "covered"?

I have a bam file of human whole exome data and also a bed file of all regions of interest.  Yes i want the average coverage of gene at threshold 30x and report by gene.

source: Stack Exchange

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

sort: votes


text:
Automated method to figure out whether a recognition site is conserved or not using the results of simulated restriction digestion

Restrict from EMBOSS Suit version 6.3.1 with the following
parameters: snucleotide1, sitelen = 4, rformat = table, enzymes =
enzymes.txt.
From the 4379 enzymes present in REBASE, we selected the 650 restriction enzymes that were commercially available, since this assay is meant to be used in any laboratory. From the 650 enzymes, 152 digest all P. salmonis sequences and only 65 recognized conserved restriction sites in the complete set of sequences, generating the same/similar restriction pattern (same number of bands and similar sizes).

https://www.frontiersin.org/articles/10.3389/fmicb.2016.00643/full
If my reading is correct, these researchers figured out the conservation of the restriction sites from the resulting virtual patterns.
The coordinates of a cut can be different in two related sequences even if the restriction site is conserved. This is because base insertions/deletions could occur before/after the restriction site.
How do i figure out if a recognition sequence is conserved or not, if the coordinates of the cut are not the same in the sequences? How do i automate this?
We can consider this as an example:
I have two sequences of the 16s rRNA gene from two strains of a bacterial species. When i run EMBOSS's restrict by the following command:

restrict -sitelen 4 -enzymes HaeIII -snucleotide_sequence 1 -sequence
sequences.fasta -rformat excel -fragments true -outfile patterns.xls

Sequences link
The result is:

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


text:
How to use hashsolo for demultiplexing hashtags?
I am trying to use hashsolo and I want to make sure that I have done things correctly.
I did the below:
#import files
A= sc.read_10x_h5('filtered_feature_bc_matrix.h5', gex_only=False)
A.var_names_make_unique()
A.layers["counts"] = A.X.copy()
sc.pp.filter_genes(A, min_counts=1)

#hashsolo
hto=A[:,A.var_names.str.startswith('HTO')].copy()
hashsolo.hashsolo(hto)
hashsolo.plot_qc_checks_cell_hashing(cell_hashing_data) #nothing gets plotted
A.obs['most_likely_hypothesis'] = hto.obs['most_likely_hypothesis']
A.obs['Classification'] = hto.obs['Classification']

Then since I have two samples for which I did hashsolo (A and B) I concatenated the objects and continued with the scanpy analysis like shown below:
sc.pl.highest_expr_genes(A_B, n_top=20, )
sc.pp.filter_cells(A_B, min_genes=200)
sc.pp.filter_genes(A_B, min_cells=3)
print(A_B.n_obs, A_B.n_vars)

and I continued with the mitochondria QC and the rest of the analysis.
Is that correct?
I am wondering if I should have done it differently. Like shown below
A= sc.read_10x_h5('filtered_feature_bc_matrix.h5', gex_only=False)
A.var_names_make_unique() Then the
sc.pl.highest_expr_genes(A, n_top=20, )
sc.pp.filter_cells(A, min_genes=200)
sc.pp.filter_genes(A, min_cells=3)
print(A.n_obs, A.n_vars)

and following the QC (filtering and mitochondria) and make the hto list and use hashsolo and then concatenate the 2 samples (A and B).


text:
mapping nucleotide sequences
I am trying to do mapping with multiple sequences against a reference genome. I tried the following code: The outputfile is .fastq.sam. I need only .sam. How do I get it?
for i in test/*.fastq;
do
 ./minimap2 -ax map-ont reference.fna $i > $i.sam;
 done

The command is fine, it is just a suffix. If you want to trim the suffix then use ${i%.fastq}.sam for the output, or even better convert the output from minimap2 to bam directly as this saves space:
./minimap2 -ax map-ont reference.fna "${i}" | samtools view -o "${i%.fastq}".bam

piping through samtools sort will save time as well as space

@gringer Thank you for the suggestion.

Honestly, I wouldn't do it with a one-liner like this.  Make a proper bash script, so you can spend a couple of lines to get the naming just right.  Then you can use the bash script over and over again.  Easier to document too if the command line used is in a file, not just in terminal history.
Also, yes, piping to convert to bam is almost certainly what you want to do.  But piping into a sort commmand is better.  And I think samtools sort will output a bam by itself if you ask it to.

I would recommend
ls test/*.fastq|sed s,.fastq$,,|xargs -i echo ./minimap2 -ax map-ont reference.fna {}.fastq \> {}.sam|sh

Replace trailing "sh" with GNU parallel for concurrent jobs on the same node; replace with asub for submitting array jobs to Slurm/LSF/SGE.
Once you get used to this, you can do a variety of tasks quickly, like renaming files in batch, gzipping all files in parallel, format conversion in batch, aligning 100 files on a cluster, etc.


text:
How to access specific SRA sequences within an SRP file?
How do you get a specific 100 SRA sequences from a Bioproject which contains 400 runs, when you don't need all of them?

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

Welcome @kirtikulanthaivel, the answer is very likely to involve the SRA toolkit (freeware), if it is doable. How to apply this is more challenging because the question isn't clear. Could you supply a reproducible example? It does not need to be the data set(s) you are interested in.

I actually wanted to know how a specific number of SRR runs, for example 100 of them from a total of 400 in an accession list can be downloaded at a time. Problem solved, thanks to a friend. Thank you all anyway.


text:
Get BioProject code of a GEO dataset
In publications the GEO code is usually shared and not the BioProject/study code.
I want to get the BioProject/Study Accession code of a GEO dataset (to later download the data from ENA).
I know that I can go to the website and see find it but I would like to automate the process.
GEOquery didn't provide it, and the data doesn't seem to be available via the FTP.
Is there a programmatic way to obtain it?

This must be doable because it has an API here https://rdrr.io/bioc/GEOquery/api/ - which is exactly what is needed here

Can you please explain more about why going to the website is not good enough? For a single GEO dataset, it would be quicker and easier to use the website.

Because I want to download the files from multiple GEO datasets and match the data with other information. The background history is that my group currently do not have a local copy of some data which I know was published in a paper, that paper has several GEOs and I want to download them and find which samples are from the missing data.
I think I identified the relevant one, but I might have missed it and I want to have a solution for the future (I think this will happen more :( )


text:
SRA toolkit and SRA- prefix
I am trying to download raw sequence data from the SRA database using the SRA toolkit (v.3.0.2). Unfortunately, the authors of the article I am reading gave SRA numbers (i.e. with the SRA prefix, not the SRR). This does not seem to work with the prefetch command:
prefetch --max-size 100000000 SRAxxxxxx -O downloads_SRAxxxxxx/
I ended up going on NCBI and manually finding the relevant SRR number for each sample; however, I was wondering if I could have done it programmatically. I appreciate your help.
Edit: Thank you for your comments, here are a couple of examples:
SRA030738 ; SRA030736 ;
Searching for them on SRA I found the SRR numbers, respectively:
SRR135604 ; SRR135602.
The paper is Hasan et al, 2012:  DOI: 10.1073/pnas.1207359109

Please [edit] your question and add a couple of example SRAs so we can test our solutions out without needing to go and look one up.

SRA is an unusual prefix for an accession within SRA; I'm used to SRR (for runs), SRX (for experiements), SRS (for samples), and SRP (for projects). Can you please provide a link to the paper, or an alternative example?

Can confirm, Supplementary Document #1 has codes starting with SRA. Searched first one in NCBI and I get two runs (starting with SRR...) from the same BioProject that seem to be correct.
Maybe the SRA-starting code is something mainly meant for the submitter? I wonder if there's a way to programmatically turn these codes into something more useful (BioProject code, run code...)

Using Entrez Direct E-utils, run the following command for your two examples:
esearch -db SRA -query 'SRA030738 OR SRA030736' | efetch -format docsum | xtract -pattern DocumentSummary -element Runs | xtract -pattern Run -element Run@acc

(I know, it seems strange to run xtract twice, but the formatting on the plain xml output is dodgy. The first pass parsed out the unwanted &gt or &lt or &quot characters. The second pass worked on the re-formatted stdout)

Perfect, the solution works. Thank you for taking the time to explain the double xtract

@Max_IT and @user3479780 Thanks for bringing this issue to our attention. It has been fixed in EDirect 19.5, which is now available on our public ftp site (See installation instructions here). Please update your installation and try again.
The problem was that using capitalized -db "SRA" bypassed normalization of the record. This was easily solved within EDirect by lower-casing the -db argument value to make it case-insensitive.
Normalization unescapes the XML-encoded ExpXml and Runs field contents to integrate them as regular XML within the docsum in situ, before passing them to the user. Thus:
efetch -db SRA -format docsum -id SRA030738,SRA030736 |
xtract -pattern DocumentSummary -element Run@acc

will now work properly, without the need for chaining two xtract commands together.
When creating structured records, NCBI internal software enforces conformation to the ASN.1 specification or XML document type definition. This is normally beneficial, blocking syntactically invalid records (and preventing new software bugs from silently corrupting production pipelines), but it limits the ability to add experimental structured data that is not (yet) in the spec.
Since xtract does not use a spec or DTD, it can find any element within the newly-unescaped contents by object name. (It still checks for properly-balanced opening and closing tags.)
Note also that you can skip the esearch step if you know the identifiers, and can instead pass them directly to efetch in an -id argument. For larger sets of identifiers, you can now pipe them from a file, one ID per line, without the need for an epost step.
EDirect has logic to split large lists into smaller chunks, sized appropriately for each database, and will internally loop through the fetch requests until all of the IDs have been processed, automatically retrying in the event of transient network or server failure.

Great, thank you, it sounds like this corrects the issue. The explanation on the usage of the new version is much appreciated.


text:
Correlate ssGSEA score with Seurat's scaled expression data
I want to check the correlation between the expression of a given "test_gene"* in each cluster within a scRNA-seq database with the expression of a signature (i.e., gene sets whose expression is known and validated to be a “marker” of a specific phenotype) in a specific cluster, that is:

Corr.: ssGSEA of signature in cluster X vs test_gene in cluster A
Corr.: ssGSEA of signature in cluster X vs test_gene in cluster B
Corr.: ssGSEA of signature in cluster X vs test_gene in cluster C
and so on...

Since I don't have much experience in transcriptomic data analysis, I was wondering what would be the best method to achieve this. The database I'm using (public) contains data from 19 different samples (patients), so I was thinking of doing a "sample-wise" Spearman's test. See this dummy example to exemplify my point:

test_gene in cluster_A
ssGSEA score of signature in cluster_X

Sample_1
2.1
2.3

Sample_2
1.63
1.76

Sample_3
0.39
0.45

. . .
. . .
. . .

Sample_19
0.95
0.81

test_gene in cluster_B
ssGSEA score of signature in cluster_X

Sample_1
1.91
1.21

Sample_2
0.03
1.85

Sample_3
3.13
0.9

. . .
. . .
. . .

Sample_19
0.09
2.12

And so on with clusters C, D, E, etc...
In this example (not real data) you can tell that gene_a's expression in cluster_A looks like could be linearly correlated with the signature's expression in cluster_X (they go "hand in hand"), but the same is not true for cluster_B. Here, a Spearman's test should output a p-value lower than 0.05 (signif.) for "test_gene in cluster_A vs signature_in_X", and higher than 0.05 (not signif.) for "test_gene in cluster_B vs signature_in_X".
All the samples in the actual database are from treatment-naïve renal cell carcinoma (RCC) patients, all of them with the same sub-type. Concerning the test_gene, it is actually a "placeholder" for each of the 20+ genes I want to test for a correlation. The cluster_X here would actually be the malignant cell cluster, being all the other clusters (A, B, etc.) the many cell types in the tumor micro-environment (TME). As for the data points, there would be 19, one per sample, but in the example I only put 3 just to make the point. There is one sample per patient, so sample = patient.
Are there any special considerations (like the type of normalization or data treatment) I should have to achieve this?
*The average of the normalized expression for that gene in all cells from the same cluster and sample.


text:
Is it normal to get truncated identifiers after applying SQL scema of RegaDB?
This question was also asked on Biostars
I am trying to set up RegaDB on my CentOS system, and I have installed PostgreSQL via conda. I followed the RegaDB installation manual to configure the PostgreSQL database. However, when I have applied the postgresSchema.sql schema from the RegaDB software package, I noticed several NOTICE messages about truncated identifiers.
Here is a sample of the output I received:
psql:postgresSchema.sql:25: NOTICE:  identifier "query_definition_parameter_type_query_definition_parameter_type_ii_seq" will be truncated to "query_definition_parameter_type_query_definition_parameter_type"
psql:postgresSchema.sql:27: NOTICE:  identifier "query_definition_run_parameter_query_definition_run_parameter_ii_seq" will be truncated to "query_definition_run_parameter_query_definition_run_parameter_i"
psql:postgresSchema.sql:29: NOTICE:  identifier "resistance_interpretation_template_resistance_interpretation_template_ii_seq" will be truncated to "resistance_interpretation_template_resistance_interpretation_te"

I understand that PostgreSQL has a maximum identifier length of 63 bytes, and identifiers exceeding this limit will be truncated. My question is whether these truncated identifiers will cause any issues with the proper functioning of RegaDB. If so, what would be the best approach to fix this issue?
I would appreciate any guidance on this matter, and if there are any additional precautions or testing measures I should take to ensure RegaDB functions as expected.
Best regards,

Hi @Vovin I would make sure RegaDB is not abandonware. Its github hasn't been updated for 9 years. PostgreSQL will have been revised many times since then and will not be the version used to develop RegaDB. Upvoted anyway.

@M__ I suppose that newer versions of PostgreSQL (as well as Tomcat) actually caused a bunch of my fail attempts to run RegaDB. But the maximum length of identifiers is a substantially base thing in PostgreSQL, as it is set by the constant NAMEDATALEN at the compilation step.


text:
Gene rank scores from LINCS database: Evangelista et al. 2022
The website https://maayanlab.cloud/sigcom-lincs/#/SignatureSearch/UpDown can be used to "Identify reversers and mimickers from over 1 million signatures by entering up and down gene sets .. " I am looking to get these "reversers" and "mimickers", aka chemical perturbagens.
So on this website, I paste my list of a selection of Up-regulated genes, in their respective box, and the same for my list of down-regulated genes. These genes are from an RNA sequencing experiment.
After the website has done its thing, I expand the search for chemical perturbagens (LINCS L1000 Chemical Perturbations (2021)). One of the data tables I can download contain something called a Gene Rank Matrix. This table consist of one row for each gene and columns of the chemical perturbagens found. The values are integers.
A small subsection of the data table would look like this:

My questions are what these integers (12, 144, 411, 641) represent. Are they how much the perturbagens (In this case REP.B025_YAPC_24H_K01_resiquimod_2.22uM) will affect the genes?
The paper describing this website searches (https://academic.oup.com/nar/article/50/W1/W697/6582159) mentions "The expression tables provide the computed CD differential gene expression coefficients for each signature, while the rank matrices provide the integer gene ranks for each signature as determined by the coefficients"
and
"Matrices that represent the gene rankings of the L1000 signatures are stored as HashMaps by the data API in random access memory (RAM). HashMaps are data structures designed to efficiently search for key-value pairs, and this implementation drastically cuts the fetch times of the gene rankings by UUIDs, thus improving the speed of the search"
Does this mean that the scores found wont have any meaning in regards to how well the perturbagens found will affect the genes? Or can some meaning be derived from this?
Thanks in advance!

Summary answer
Although I'm unhappy about how the paper is written and possibly the methodology, the following sentence is likely (almost certainly) to be the answer.
The coefficient (aka weight) is how important that gene is for the classification of that signature. Thats it! :-)
Essentially, there are loads of genes involved in identifying which genes represents which signature. Not all those genes are equally informative. The positive or negative skews in that calculation, called a "classification", are identified through their weights (aka coefficients). Thats simply the answer. The bigger the coefficient the more important that gene is in the classification, i.e. its associated with the signature.
There is no need to be more complicated than that. If you are interested in a critique of that paper, which I had to do to get to the answer and thereby understand the analytics of that paper, this is presented below.

Critique of Evangelista et al. 2022: Nucleic Acids 50:W697–W709
Note 1 I don't have time to read the paper in absolute detail, particularly the appendices, where some of the information needed may have be placed.
Note 2 In the answer below, a "expression coefficient" is what I've called an "expression weight" its the same thing.
Weightings and analytical concerns I do understand weightings in deep learning, particularly machine learning. I read the paper and the ANN (artificial neural network) was very basic. I don't believe it had multiple layers, because it would have said so and it was a very odd description for a deep learning algorithm. It would also be odd to have a single layer in deep learning, very odd. The loss function is also not described at all which is extremely important in deep learning (it will not work otherwise). Essentially, what the authors call "deep learning" appears to a machine learning neural network using an Adams gradient descent.  Adams gradient descent is a standard approach in deep-learning admittedly. Thereafter in the paper the details were very sketchy, maybe they are in the appendix (should be checked). Any deep-learning, machine-learning requires a demonstration of the robustness of the training and this appears absent. The two basic stats are accuracy and AUC-ROC, which were definitely absent. Universal stats are used to quickly understand how good the algorithm is in comparison to other deep-,machine- learning algorithms for other data sets.
Analytical overview Generally I would assume the authors have performed a supervised learning algorithm, but this may not be correct (below). If it is supervised learning they have performed they would have identified their 'signatures' against a very large data set, then trained the different types of signature against that data and tested it. The learning could be terrible - but it's just not stated at all. I'm surprised because Nucleic Acids is at the pinnacle of algorithm publications.
Summary. The authors performed deep-learning (which looked like a neural network machine learning algorithm). The description is at best poor, but further information could be present in the appendix. Therefore its not absolutely clear what a "weight" would represent because there are two possibilities.
Weights possibility: unsupervised learning The one thing that might have happened is they've performed an unsupervised learning algorithm and the weights are a sort of "maximum variance". This very trendy. That would explain the lack of accuracy and AUC-ROC type stuff. IF this is true and I don't believe it's explained - I could be wrong - then the weights represent maximum information. What does that mean? Good question, its context specific, its the gene which is most important in sectioning up the data - e.g. housekeeping gene which is required 'all over the place', not just part of the cell cycle is an example. Basically, in unsupervised weighted output high ranking weights are a mix of really boring and really brilliant stuff. So you just go through it with biological understanding, is this gene boring yeah/neigh, is this one - wait wow! unexpected ... That is actually fashionable analytics :-)
Weights possibility: supervised learning Is already answered at the beginning. This is almost certainly the correct answer - it would be the orthodox approach - and we just accept the robustness checks for the algorithm and possibly the algorithm itself were simply insufficient.
Why isn't the paper clear (at best)? What I assume has happened is the team have brilliant frontend and genomics team backend. The analyst in the middle is unlikely to be a specialist and thats why it's poorly described, or possibly performed. I assume if the reviewers were of similar composition then that's how its got through Nucleic Acids. Thus, I would personally treat the output with caution because deep learning (machine-learning) will be the essential bottleneck between their backend development and your data. Maybe "bottleneck" is not quite an accurate description, because its speeding up the classification, so maybe a "gate-keeper" is a better description.


text:
Detect monomer from from multimer structure in PDB files?
I have a bunch of .pdb files that consist of protein-peptide structures. However, before start working with them I need to clean them. Now, I can easily remove water, ion, or other stuff like BMET, BASN, BARG info, and anisotropic temperature information. Now, if the PDB flies consist of more than one instance of protein-peptide then it will be hard for me to work with them. The chains of protein and peptide numbers are arbitrary. For example, in a PDB file, a protein may have 4 chains(A, B, C, D) of proteins and 4 peptide chains(E, F, G, H). In that case, it may have 4 instances and we can select A and E to get a single instance from it. But, in pdb file, protein may have 4 chains(A,B,C,D) of proteins and 2 peptide chains(E,F). In that case, it may have 2 instances.
Now the question is, how to select one instance of protein-peptide from the PDB files without doing manual curation?

set a cutoff lenght for peptide, everything smaller thn lenght X  would be a peptide.  Then to get the protein-peptide pairs , I guess you'll have to calculate wich protein chain is the peptide closest one  and assign all of the peptide to a single protein chain. Does it sound reasonable ?? No idea how to accomplish all of that though


text:
Filter VCF file for most differentiated SNPs (Fst) in VCFTools
I have a filtered VCF file that I want to subset for the 500 most differentiated SNPs based on estimated Fst values (calculated in vcftools). I am not completely sure how to approach this. I have calculated the fst values across all populations (10 in total), but am not sure how to specify that I only want 500 SNPs with the highest fst value.
I used VCFTools, but would BCFTools be more appropriate?

Please [edit] your question and give us an example we can use, and the output you would want from that example. We cannot help you filter data we cannot see. Show us a valid VCF file (that means all headers plus a few variant lines) and then tell us which variants you want to keep from it.


text:
How to export the list of soft-clipped regions from a region of interest in IGV
I want to export the list of soft-clipped regions into a csv file from a region of interest in IGV. Additional information regarding the soft-clipped sequence, coordinate and nucleotide length will be helpful.

PFB the output file format

Gene
Variant

FLT3
chr13:28034082C>CAAACTCTAAATTTTCTCTTGGAAACTCCCATTTGAGATCATATAG

FLT3
chr13:28034113T>TTTGAGATCATATTCATATTCTCTGAAATCAACGTAGAAGTACTCATTATCTGAGGAGCCGGTCACCTGTACCACC

Does this have to be done from IGV?

Also, please [edit] your question and add an example of the output format you want to see. What exact fields should be in the csv?


text:
How to identify all descendents of a branch on a phylogeny
I'm trying to simulate the evolution of binary traits on phylogenies in R to end up with a dataset of the traits held by each tip. Starting at the root, at every timestep, new traits are added with a certain probability, and existing traits are lost with another probability.  This happens independently for every branch that exists at that timestep before moving to the next timestep. Those traits are added/removed from the final trait vectors of each species that is descended from the current branch/timestep.
The issue is I don't know how to identify which tips (final species) originate from which current branches (at a specific timestep).
Is there a function for identifying, on a phylogeny, at a specific time, tips that share a current ancestor/tips that are descended from the currently selected internal branch?
Or, do you have a better way of organizing the trait data through time?
I've tried using other software (TraitLab in Matlab) but not sure that it does exactly what I want, which is why I'd rather code it in R. If you have other suggestions of software to use, that would also be much appreciated! Unfortunately most trait evolution simulation software require you to know the total number of traits from the start, whereas I want new traits to appear and be lost.
I want to do something similar but different to getMRCA. Starting with one selected species, at a specific time in the phylogeny, I want to identify all other species from which it hasn't yet split.
Traits are gained under this model. The way I planned it is that every species has a vector of traits, each position in the vector being for 1 unique trait (the positions indicate the same traits across all species). So when a trait is gained, a "1" is added to all vectors of species descended from that branch, and a "0" is added to all other species. When a trait is lost, the "1" for that trait in all species descended from that branch is replaced with a "0". I'm using the ape library, I just can't figure out how to identify which species are descended from a specific branch.
Do you know of a function that can do that?
Thanks!

Okay there are other ways in using complete packages outside R, but if it's strictly R I just dunno. getMRCA is a good idea, but it's the opposite of what is wanted here I think. I believe the OP wants MRCA -> tips (OTUs), rather than tips -> MRCA.  Possibly from the bio.infer library get.otu? Just guessing.

Maybe you want extract.clade in ape that gets the subtree at a specified node? Its description is hidden in the drop.tip documentation.

It sounds like you do know the total number of traits from the start; it's just the length of the vector of traits. As I understand it, what you're calling gaining or losing a trait, is really just changing the state of the trait between 0 and 1.


text:
Where do I get a large reference VCF?
I would like to download a large .vcf file containing many (hundreds or thousands) of samples. Ideally, I would download different population-specific .vcf files, but the ability to sort/filter by ancestry group is fine. Where do I get such a file? I prefer GRCh37 for consistency with other files I'm using.

Update: I tried the HGDP file--I'm trying the 100 genomes one now, in addition to trying http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000_genomes_project/release/20190312_biallelic_SNV_and_INDEL/ now.
I downloaded the file, but when I opened it it seemed to have no samples (instead of having many samples). Here is a full row:
 1  10145   .   AAC A   151.43  PASS    AC=4;AF=0.0127389;AN=314;BaseQRankSum=0.72;ClippingRankSum=0;DP=24144;ExcessHet=160.275;FS=0;InbreedingCoeff=-0.1558;MLEAC=5;MLEAF=0.002604;MQ=1.63;MQRankSum=0.764;NEGATIVE_TRAIN_SITE;QD=3.52;ReadPosRankSum=0.085;SOR=0.458;VQSLOD=-0.9391;culprit=QD;VQSRMODE=INDEL;NS=157;ExcHet=0

This question has also been asked on Biostars

@BigMistake please [edit] your question and give us some context on what you need. You have now mentioned you also need "ancestry information" (which isn't very clear, what kind of information?), and of course you need samples, but please explain all this in your question so we can know what you are after. Ideally, tell us something about how you will use the file so we will be able to give you better answers.

Hi @BigMistake thank you for posting your answer. Could you kindly mark this as "accepted"?

Solution: Use https://bochet.gcc.biostat.washington.edu/beagle/1000_Genomes_phase3_v5a/b37.vcf/, then to get the ancestry information (i.e., the ancestry group of each sample, e.g., "HG00566" or "NA19758") use https://www.internationalgenome.org/data-portal/sample.

For anyone who arrives here looking for a large reference VCF but not wanting GRCh37, you might consider the HGDP+1KG subset of the Genome Aggregation Database's dataset. You can find a VCF file per chromosome aligned, I believe, to GRCh38, here.


text:
Easier / lazy way to convert an existing bash pipeline to nextflow pipeline?
Not going to lie, my head hurts when I read the NextFlow documentation (No dig on the documentation, I think my head is losing its processing capacity lol).
I had written a pipeline in Bash, essentially iterates over a text file containing full paths of R1 and R2 reads, followed by a series of commands and liberal use of gnu parallel to improve efficiency.
Is there any template I can follow, where I can slot in my commands and get things running quickly?
EDIT - The skeleton of pipeline goes as follows (As I'm not sure if I can put our inhouse pipeline in its entirety here) -
Bash code which loops through R1 and R2, aligns with BWA-mem2 and sorts the bam, then using gnu parallel, I essentially follow the steps as outlined in GATK Best Practices (looking to drop BQSR, but the steps remain) to obtain a gvcf file. The steps in Best Practices have been turned to functions, rather than each step being called sequentially.

Please [edit] your question and show us your code. Nobody can help if they don't know what your code actually does. But, that said, every pipeline is different so I really doubt there will ever be a "template".

Hi KarthikNair could you kindly supply the code or else a reproducible example? Otherwise there are a lot of NextFlow questions and answers on Bioinformatics SE, particularly answered by user @Steve, and those would likely provide pointers. To be honest it's not possible to answer the question as it stands in my personal opinion. I've a reasonable understanding of NextFlow, but I'm not an expert.

I have mentioned the brief gist about the pipeline, essentially we are following GATK best practices for variant calling (humans).

If your objective is to convert Bash scripts into Nextflow code, you might find Viash (https://viash.io) a particularly useful open source tool. It automatically generates independent Nextflow modules from each of your scripts, along with the corresponding Docker containers for each module.
Docker containers are generated on the fly based on your config description in Viash.
These modules can be coupled together to build a (modular) Nextflow pipeline.

Thanks for answering the question as written! Super helpful!

Generating Docker containers on the fly is part of the behaviour of NextFlow per se, which is why it's so cool. However, I know zero about Viash. It would be useful if someone could give feedback on Viash.io because it may not possibly work as advertised here .... The question was flagged for hidden issues, i.e. it may not work as advertised, so I would recommend proceeding with care.

The gist would be something like:
fastqs = Channel.of(...some annoying groovy code to make this a list of tuples...)
workflow{
  align(fastqs)
  makeGVCF(align.out.bams)
  ... collect() the output from above joint call them all...
}

process align {
  input:
    tuple val(sampleName), path(reads)
  output:
    tuple val(sampleName), path("${sampleName}.bam"), emit: bams
  """
  bwa ... | samtools sort ...
  """
}

process makeGVCF{
  input:
    tuple val(sampleName), path(bam)
  output:
    path "${sampleName}.gvcf", emit: gvcfs
  """
  ...haplotype caller code...
  """
}

The usual way would be to use either the fromFilePairs factory method, or to read from a samplesheet using the fromPath factory and the splitCsv operator, e.g. Channel.fromPath( params.samplesheet ).splitCsv( header:true ).map { row -> tuple( row.sample, [ row.fastq1, row.fastq2 ] ) }.set { readgroups }. Then align( readgroups ).

I really thank you both for your inputs. But my head started hurting when I reached the end of the comments. I guess I just have to grit my teeth and push through it.

Since you asked about templates, if your task is something rather generic, there is a good chance that a module, or several, that you can use to build your workflow exists in https://nf-co.re/modules
Generally I would recommend taking a look at the source of nf-core pipelines (e.g. https://nf-co.re/sarek for variant calling) to get an idea how to get your samplesheets in.

Nextflow (https://nextflow.io/) is a workflow manager / language used to implement analysis pipelines, nf-core is an open source collection of modules and pipelines written for / in nextflow

Hi @Niklas thank you - I see your point. Its leveraging NextFlow and providing a link to a NextFlow template on GitHub. Got it. FYI I can use NextFlow. Its Docker container on the fly, is neat, but I agree with DevonRyan its Groovy (an obsolete language) I personally think is a bit weird.


text:
Tidying MEGAN Taxonomy for Python/R Analysis
I am analysing some WGS data in MEGAN and would like to do some additional analysis in Python/R
I am having trouble Tidying the Taxonomic data in a format which would be conducive to this. Originally I tried exporting Taxon Path to count and using the Split function in Pandas to break up taxonomy into new columns as this worked well for KEGG Pathways. Due to varying specificity of taxonomic assignments this didn't work.
I've had more luck Exporting Data as a BIOM file format and converting it to TSV resulting in a file with the following column containing all taxonomic info (See below)
Ideally , I'd like to generate columns for; Domain, Phylum , Class , Order, Family and Species by accessing the  features contained in Taxonomy. I'm unaware of any functions which would let me do that.
an idea of what I'm trying to achived is

Domain
Phylum
Class
Order
Family
Species

Bacteria

uncultured bacterium

Bacteria
Bacteroidetes
Bacteroidia
Bacteroidaceae
Bacteroides
Bacteroides acidifaciens

From

Taxonomy

d__Bacteria; s__uncultured bacterium

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides acidifaciens

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides caccae

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides cellulosilyticus

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides faecis

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides finegoldii

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides fragilis

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides intestinalis

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae;

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae;

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides salyersiae

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides stercoris

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides thetaiotaomicron

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides uniformis

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides xylanisolvens

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides cellulosilyticus CAG:158

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__uncultured Bacteroides sp.

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides sp. 2_2_4

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides sp. 3_1_23

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides sp. D2

Thanks @dunc4n, pandas is a good solution. It's fairly easy. I will not have time in the next two weeks to answer this, particularly coding from an inline table. If you supplied some code that would be helpful, particularly the data you want as a Python list rather then a inline table.e

Save the following as a python script (e.g., named splitty.py). Make the script executable by doing chmod +x splitty.py. Then, if your MEGAN data are in a file called taxonfile.csv, run splitty.py by doing ./splitty.py taxonfile.csv. Get help on the script by doing splitty.py -h.
#!/usr/bin/env python3

"""
Splitty

Use this to split up a taxonomy table
"""

def argsparser():
    """
    argsparser
    """
    import argparse
    parser = argparse.ArgumentParser(prog="splitty",
                                    formatter_class=argparse.ArgumentDefaultsHelpFormatter,
                                    description="Uset splitty to get taxon table.")
    parser.add_argument("table", help="Taxon table in semicolon separated format.")
    return parser

def parseline(line_in):
    """
    Parse each line_in as return a dict
    """
    transl = {"d": "Domain",
             "p": "Phylum",
             "c": "Class",
             "o": "Order",
             "f": "Family",
             "g": "Genus",
             "s": "Species"}
    return {transl[val[0]]: val[1] for val in [i.split("__") for i in line_in]}

def main(args):
    """
    The main routine.
    """
    import pandas as pd
    df_list = []
    with open(args.table, "r") as input_handle:
        for line in input_handle.readlines():
            linelist = parseline(list(filter(None, [i.strip() for i in line.rstrip().split(";")]))[1:])
            df_list.append(pd.DataFrame(linelist, index=["Taxon"]))
    print(pd.concat(df_list).to_csv(sep="\t", index=False))

if __name__ == "__main__":
    main(argsparser().parse_args())

Thank you this got things sorted great !


text:
How can I detect the Apo B gene mutations in my raw dna data?
Familial defective apolipoprotein B-100 (FDB) is an autosomal dominant disorder associated with increased plasma cholesterol concentration and may thus increase the risk of premature coronary heart disease. To date, three apo B mutations, apo B R3500Q, apo B R3531C and apo B R3500W that cause defective apo B, have been reported.
I have my raw DNA data from 23andMe (where each line corresponds to a single SNP). How can I search it for those related Apo B gene mutations?
A single line in my DNA data looks like this:
rs12567639 1 756268 AA

I have a programming background, but none in the bioinformatics area.

There are a fair number of members that could answer this. They might think this is too difficult to easily explain because this is an easy answer, but it's alot of explaining given your background. It would be helpful to try and provide the key information, a couple of lines of the output placed into the answer - you've probably got a simplified VCF. The reference genome used is really important (the company should supply and the original VCF). If there's no answer then VarSome is one company that could do this (no affiliation).  I'm no longer involved in eukaryotes apart from immune stuff.

Anyway its a good question (in my opinion) and good luck with the analysis. Upvoted.

For some reason 23andMe uses genome positions from the previous genome assembly. Anyway, what counts are the identifiers, like: rs12567639

I think you mean these.  The amino acid numbering has changed, but that happens as we discover more transcript variants.
http://useast.ensembl.org/Homo_sapiens/Variation/Mappings?db=core;r=2:21005788-21006788;v=rs5742904;vdb=variation;vf=184237596
http://useast.ensembl.org/Homo_sapiens/Variation/Mappings?db=core;r=2:21005696-21006696;v=rs12713559;vdb=variation;vf=186851457
http://useast.ensembl.org/Homo_sapiens/Variation/Explore?db=core;g=ENSG00000084674;r=2:21001429-21044073;v=rs144467873;vdb=variation;vf=251617275

If you are interested in ApoB gene variants you can:

go to https://www.ensembl.org/Homo_sapiens/Gene/Variation_Gene/Table?db=core;g=ENSG00000084674;r=2:21001429-21044073

sort the table by Clin.Sig. column

select the name of the variant, i.e. rs5742904

search for that variant in your 23andMe file

compare the variants from the ENSEMBL with what you got in your result

https://www.ensembl.org/Homo_sapiens/Variation/Explore?db=core;g=ENSG00000084674;r=2:21001429-21044073;v=rs5742904;vdb=variation;vf=184237596
But I assume all of this is moot since 23andMe should give you some summary of the results and not just the raw file, right?

23andMe doesn't give any health recommendations to their clients outside of the US. All I got was the raw file, some ancestry data and human traits probabilities.

OK, just to make sure: did they send you any non high cholesterol related potentially pathological trait probability?  Anyway, if you want to extract relevant markers/sequence variation from your results or rather intersect these the above ApoB pathogenic markers  not one by one you probably should start with getting either command line on Mac or WSL (Windows Subsystem for Linux).

To give you some idea about the task: one can obtain all variants for the human APOB gene and their putative clinical significance using BioMART from ENSEMBL and ENSG00000084674 gene id. The result is a TSV (tab separated file). Selecting the column with variant IDs and getting the relevant lines from your 23andMe text file is trivial on Linux/Mac command line (awk and grep commands). I am not familiar with Power Shell but looks like Select-String does not have the same functionality as grep.

No pathological traits from 23andme whatsoever. I'm fairly proficient with the command line. Just need to devote a couple of hours analysing your answer :)

of course 23andMe won't give you health recommendations.  That would be irresponsible.

@PawełGościcki here are the annotated APOB gene variants from ENSEMBL 109: https://drive.google.com/file/d/1CECqvG_oWDi_MoV_h7GsjL9aqnKPskeC/view?usp=sharing  You just need to check how many  of the Variant_name identifiers are present in your 23andMe file, look at  pathogenic in column 3. I will need to dig into the 23andMe result format to get back to you.


text:
Low fraction reads in cell metric in snRNA-seq data
For some of our snRNA-seq samples we are finding low fraction of reads detected in single-nuclei rna-seq samples from cellranger, while the other metrics are perfect.
While I understand this could be caused by library preparation and abundance of ambient RNA. Just wondering if this behavior is expected for single-nuclei datasets. My confusion is specifically coming from the fact that 10X website also has such a sample that has low number of fraction reads in cells.
https://cf.10xgenomics.com/samples/cell-exp/7.0.0/5k_mouse_brain_CNIK_3pv3/5k_mouse_brain_CNIK_3pv3_web_summary.html
Would love to know if this behavior is expected or seen often who are more experienced in snRNA-seq. The metric documentation page is here
Here is snapshot from 10X report

What do you mean by low fraction? Fraction of what (i.e. what's the denominator)? What were you expecting? What did you see instead?

Great questions, sorry for being vague, uploaded the screen-shot for cell-ranger. low fraction - 39%, expected according to 10x 70% (if you follow the link). I am not sure how they are computing the metric though, I added the man page from 10X

Yes, "low fraction reads in cells" is more common with snRNA-seq from frozen samples, I wouldn't worry about it by itself (you're just throwing away some sequencing data). You should still check though that the other QC metrics are good. In particular, sometimes such samples have an unclear seperation of cells from empty droplets on the UMI rank plot and the consequent estimated cell count is off, requiring alternative approaches to distinguish cellular droplets from ambient RNA etc.

Thanks a lot Chris. I am thinking about more sophisticated pipeline like CellBender etc. I just wanted opinions for professional more experienced than me because I have not dealt with non-standard snRNA-seq data a lot.


text:
Salmon Pseudo count when dealing with male and female RNA-seq data
I've generated a quant seq data that I intend to use to compare male and female gene expresion, with a focus on sexual chromosome.
For my species (three-spined stickleback), it is a classic XY sex determination system, and i have a male and a female reference transcriptome, including or not the Y chromosome.
I am not sure how to handle that with SALMON, and I have 2 questions in particular:

Given that I use two different transcriptomes depending of the sex, can I still compare the obtained normalized pseudocount of would it affect the estimation too much?
I want to compare expression levels of genes present in both the X and Y reference (i.e genes that are not recombining but have not been lost or degenerated on the Y). For that, I'd need to add the gene counts of the Y and X copy. Is it Okay to add  pseudo count from different gene copies, or will that be issue? I do not understand the pseudo-count statistical model enough to have an idea about that

My other option is to go with a STAR / HTseq pipeline, but my mapping rates seems about 20% lower with this tools.


text:
regex: samtools command to "refine" PacBio IsoSeq data?
I am working with preprocessed data (IsoSeq PacBio) and I cannot understand one of the steps:
samtools view -b -e [rq] >= 0.9 flnc.bam

What exactly was filtered from the bam file?
To obtain flnc they used the standard IsoSeq pipeline (lima -> isoseq3 refine).

The samtools view command you've provided just uses a filter expression to discard records that do not meet the specified expression criteria. In this case, the expression is looking for alignments where the value of the 'rq' (read quality) tag is greater or equal to 0.9. The isoseq tags documentation says that the read quality of the alignment is really the predicted accuracy of the polished isoform. Note that the command produces BAM output, so be sure to redirect stdout or write to a file for example using:
samtools view -b -e '[rq] >= 0.9' -o out.bam flnc.bam

thank you very much!


text:
How do you go about grabbing proper genotyping infomation from dbSNP using entrez API?
I made an R function to grab rsID genotypes from dbSNP using the entrez API. For example, here is the xml output that I am trying to parse:
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=snp&id=rs2656176&retmode=xml
the 'SPDI' tag contains the genotype info. so I can easily use xpathSApply to get this tag for many IDs. But for this dbSNP there are two different values, T>C and T>A.
The T>C change is the one that has more support in all pop freq databases, and this is easily grabbed from the dbSNP webpage:
https://www.ncbi.nlm.nih.gov/snp/rs2656176
The "Alt Allele" frequencies show A being nearly 0, while C is the major alt allele.
However, I don't understand how I can parse the data from the XML from above using this... the 'global maf' tag in the XML shows some MAF values, but the 1000G value is only for the reference allele, not the alt.. There are other dbs that show the alt "C" allele but that might not always be the case.
Does anyone have an idea of what tags I should look at? Or is there a better way to hit a different API to get this data instead?
Thanks!


text:
maxquant error: Unhandled Exception: System.IndexOutOfRangeException: Index was outside the bounds of the array
An error occurred while I was running the program maxquant (2.0.3):
Reading search engine results (SP)
Unhandled Exception:
System.Exception: Exception during execution of external process: 1228404
Unhandled Exception:
System.IndexOutOfRangeException: Index was outside the bounds of the array.
at MsLib.Search.AndromedaPeptide.FromString (System.String line, BaseLibS.Ms.Data.Protein.ProteinSet proteinSet) [0x000a9] in :0
at MsLib.Search.AndromedaQuerySp.Add (System.Int32 count, System.String line, BaseLibS.Ms.Data.Protein.ProteinSet proteinSet, System.Int32 minLength) [0x00001] in :0
[ERROR] FATAL UNHANDLED EXCEPTION: System.IndexOutOfRangeException: Index was outside the bounds of the array.
at MsLib.Search.AndromedaPeptide.FromString (System.String line, BaseLibS.Ms.Data.Protein.ProteinSet proteinSet) [0x000a9] in :0
how can i fix it, thank your very much!

Thank you for your reply. I only changed the modification and added the fasta  raw file  in the configuration file, Other parameters were default.

Dear @jingzhao thanks for posting back. Could you kindly "accept" your answer (green tick button) - it does help with our site stats

****this was the xml file that i configured to maxquant. i changed modification and added the fastq and raw file. Other parameters were default. however, 'Unhandled Exception: System.IndexOutOfRangeException: Index was outside the bounds of the array' was reported when I run the command 'maquant mq.xml'. but, When I only kept the N modification, the program ran successfully, but a small amount of peptides were identified.


text:
Where is the galaxy results file path?
I want to wrap galaxy into a docker image, and I want to delete files that are older than a certain date on a regular basis, so I want to know how to generate the files and where to store them.


text:
Mapping statistics from the bam files
I would like to find the mapping statistics from the sorted bam files.
Samtools flagstat gives the output only for a single file. What is the easiest way to find the mapping statistics for all sequences? How do I find the mapping percentage for all sequences? Is it possible to use the samtools flagstat if I merge all the files in to a single bamfile?

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

Hi @Kam and welcome. What sort of stats are you thinking of? E.g. coverage, average read length? Also what sort of sequence data is this?

I would like to find the mapping percentage for all the reads after mapping. I have around 300 bam files after mapping. Is it possible to find the mapping percentage for all the reads?

run it in a loop. something like for f in bamfiles; do ...; done

Which mapper are you using? Have you checked the logs of your mapper to see if these metrics are already provided there? Are you sure you are retaining unmapped reads in the final BAMs in order to calculate such rate?


text:
What statistical test to apply for DE after CibersortX deconvolution?
This question was also asked on Biostars
I am running CibersortX in high-resolution mode (which yields estimates of gene values per sample).
After that, I want to perform DE between two conditions on the resulting gene estimates.

What statistical test would one need to apply to perform DE after high-resolution mode de-convolution?

If I use the output of DE for ranking purposes only (in pre-ranked GSEA), does it make sense to use a simple t-test?

The data is in TPM (CibersortX received TPMs as its input). I transform to log2 scale before t-test. There are only 7 samples available (4 control, 3 treatment)
Here are some genes which get the smallest p-value after t-test (but do not pass FDR)
ETFBKMT 0.08727244  0.2549699   0.1277664   0.000000    0.597810    0.4984627   0.4672444
AFF3    3.28551359  2.9969676   2.8793547   3.013827    2.483018    2.5070328   2.3119873
PDLIM1  4.26118568  4.0109785   3.9223976   3.856899    3.388990    3.5155023   3.2931998
MYOF    6.17754486  5.9286025   5.7452668   5.845756    5.276884    5.3981318   5.1251924
LINC02210 3.22922278    3.1785117   3.5162067   3.307364    2.736807    2.9091521   2.8052582
CD274   5.20308488  5.0063801   4.8715090   4.922224    4.440247    4.6192961   4.5082217

And these are the values before log
ETFBKMT 1.062360    1.193311    1.092601    1.000000    1.513417    1.412707    1.382466
AFF3    9.750753    7.983202    7.358209    8.077041    5.590659    5.684498    4.965666
PDLIM1  19.175412   16.122220   15.162099   14.489126   10.475814   11.435934   9.802840
MYOF    72.381287   60.909803   53.641094   57.510614   38.770421   42.169610   34.900902
LINC02210 9.377626  9.053727    11.441519   9.899558    6.665934    7.511766    6.989834
CD274   36.837031   32.141830   29.273209   30.320557   21.709387   24.578008   22.756735

And that is the a distribution of the expression values (after deconvolution) before log2, in one of the samples
 Min.   :    1.00
 1st Qu.:    5.01
 Median :   18.03
 Mean   :   89.08
 3rd Qu.:   47.82
 Max.   :34406.87

To specifically answer the question - a T-test is a parametric test which requires a normal distribution. So it's not a good idea in context to this data. However, the non-parametric equivalent of the T-test is the Wilcoxon and thats absolutely perfect and described below together with the entire rationale for parametrics versus non-parametrics. Basically, I agree with your approach, except I would simply shift to non-parametric statistics.

I noted the response on Biostars from the link recommends DSeq2 and EdgeR. Both tests along with the additional test described below, along with the T-test, are used to identify whether each gene in the analysis is significantly up-, or down-regulated and can then be ranked. The underlying issue is the validity of the tests associated P-values.
Both DSeq2 and EdgeR have been notably criticised by Li et al 2022, in a paper entitled

Li et al 2022 Exaggerated false positives by popular differential expression methods when analyzing human population samples" Genome Biology: 79

Li et al (2022) needs to be read carefully before making a choice. Li et al 2022 is recommending using a non-parametric test, via Wilcoxon, which is implemented in the method NOIseq from the R-Cran Bioc library (described below). Genome Biology has a current impact factor of 18 (exceptional). I was told that one member of the DSeq2 team produced a later paper advocating non-parametrics methods.
The underlying issue is the dependency on a parametric distribution for both DSeq2 and EdgeR. Basically, the central issue is the validity of parametric statistics (assumed distribution) versus non-parametrics statics (no assumptions about the distribution).
In parametric statistics outside DE, i.e. in general, when these are used numerous tests are performed to assess the fit of the data set to the parametric distribution assumed. For example, for ANOVA it's the Fmax test. For a T-test there's a number of tests for the normal distribution. If the data set fails the test, a multitude of transformations are possible usually aimed at normalising the variance.
From this perspective the concern is no comparative tests were devised for DSeq2 or EdgeR given the age of the papers DSeq2 is from 2014 Love et al (2014). If there's no test how could the data use alternative transformations to fit the assumptions? Again this is unusual given the age of papers because such as tests would be expected.
Non-parameteric DEG
The non-parametric method mentioned, NOISeq, is from the R/ Bioc package by Tarazona et al 2015, Nucleic Acids:43. This was the approach use by Li et al 2015 and compared against both DSeq2 and EdgeR. There are more recent examples, again I believe one of the DSeq2 team produced a non-parametric DEG as well.
The GitHub page for Li et al 2022, is https://github.com/xihuimeijing/DEGs_Analysis_FDR

I understand the criticism recently stated about Wilcoxon (non-parametrics in general) is their reduced power (response on Biostars). What happens is you get fewer significant 'hits' that are under or over-expressed genes. This is the basic trade-off between parametrics vs. non-parametrics. The key issue is if you've trapped your key results within a non-parametric test, that is a very solid result. If the equivalent results are trapped using EdgeR or DSeq2 there is always going to be question mark whether it's a false positive.
When to use EdgeR/DSeq2: if you are doing alternative statistical tests on the expression data which are directly comparable with under- and over-expression (e.g. ML and DL analysis) thats okay because we know Edge/DSeq2 gives false positives. When I've done this there's a 'centre ground' which can give conflicting result, its the central quartile type area, i.e. 25 - 75% region (to be honest it can be a bit outside that).
If you are singly dependent on EdgeR, DSeq2 or NOISeq it is better to be conservative in my opinion, i.e. fewer hits which you are more certain about.
Final point The final issue is simply how close to the DSeq2 and EdgeR data sets is the data you are working with, i.e. same species? What I would caution is these distributions have been applied broadly across species. The extent of the false positive rate must vary in my opinion.


text:
Protein fold topology search web-server
I'm finding a web server (or standalone program) that allows the upload of a custom PDB for protein fold topology search.
Do you know any of them?
More details:
Actually I started off with Foldseek but have a question in mind:
how do I know if I have a novel fold? (because similarity search will often return partially similar hits).
Then, I tried to think of topology, as an illustration:

which simplifies structures to spatial arrangements of secondary structures (or topology).
I want to know if there's any database of protein fold topologies and tools for topology similarity search.
ADD COMMENT • link • edit •

Have you tried foldseek?

I started off with foldseek but I need lower resolution search.

There are a few, but the main candidates for this would be :

CATH - http://www.cathdb.info/
SCOP - https://scop.mrc-lmb.cam.ac.uk/
Ptgl - http://ptgl.uni-frankfurt.de/

CATH and Ptgl have search services on their website, but not sure about SCOP.


text:
WES variant calling with DNBSEQ-T7: technical quality assessment
I recently received whole exome sequencing samples that were sequenced on an MGI sequencing instrument, DNBSEQ-T7. I am interested in somatic variant-calling on the paired tumor-normal samples. As this isn't the canonical Illumina sequencer, an issue that I am concerned with is that the technical artifacts arising from the sequencer are not yet well characterised and this may lead to spurious artifacts in the variants.
I asked the company that did the sequencing for a panel of normals (PON) for the instrument to run Mutect2, but they simply replied that they use the standard GATK PON for their somatic variant calls.
How should I go about assessing whether the variants are "biased" in any way, as a result of the instrument's technical properties? Is there a way to quantify the magnitude of these technical errors from the VCF files?

What exactly do you have? Fastq files? Bams? You mention VCFs, how were those generated? The PONs are mostly used to try and guess if a variant is germline or somatic, as far as I know, why would they help with artifacts?

The fastq files contain the raw seq information. They also performed some secondary analysis as "part of the package", and this includes the bwa aligned bam files and a germline haplotypecaller vcf output. The purpose of the PON is to infer possible technical artefacts - see https://gatk.broadinstitute.org/hc/en-us/articles/360035890631-Panel-of-Normals-PON-

Okay it seems the company has given you the VCF (maybe), you want to re-run it via Mutect2.  Mutect2 is part of the GATK suite, it's an upgrade, and they are refusing to run the appropriate PON. I think this is more personnel management than a major technical hurdle.
I would reply to the tech and CC the technical director  of the company explaining Mutect2 vs. regular GATK PON. You simply use a standard process of escalation within technical justification. Maybe send them the commend to demonstrate how simple it is:
     gatk Mutect2 \
     -R reference.fa \
     -I tumor.bam \
     -I normal.bam \
     -normal normal_sample_name \
     --germline-resource af-only-gnomad.vcf.gz \
     --panel-of-normals pon.vcf.gz \
     -O somatic.vcf.gz

The company will have to justify their position regarding PONs or accommodate your requirements either now or the future. There may simply be a technical misunderstanding between you and the tech.
You've paid for the service - thats the standard you expect. It's easier for them to assist you than all the agro of a complaints procedure and I'm pretty sure they'll help.

Thanks! Actually fyi, the company provided me with germline VCF files as part of the secondary analysis, which isn't useful given that I am interested in the paired timor-normal comparisons.

@kane9530 good luck and glad it helped. You are working in a really interesting area of research and Broad Institutes' tools/work are second-to-none (really good). I understand your concern given how meticulous the Broad Institute is and the technical issues in separating cancer cell DNA from healthy cells within the tumour (tumour) sample.


text:
Blast output file only shows 500 lines -outfmt 6
I had created databases of different sets of metagenome datasets - one with 6 runs, other with 48 runs, another with 100 runs, etc using the accession list for each of these datasets and makeblastdb command. Since I wanted to check the abundance of one specific 16S rRNA gene sequence, I ran a blastn command using this 16S sequence (F7) as the query against each of the created databases with a txt file (F7_droughtS.txt, F7_Dy.txt, so on) as the output. Ideally I should have got differently sized files with different number of OTU matches. I did get different files since the SRR ids were different for each database, but the results did not go beyond 500 lines for any database search. Can anyone explain why this is so?
Following is the command I used for blastn of:
sequence.fasta against database firstdb :
blastn -query sequence.fasta -db firstdb -outfmt 6 -out output1.txt

Also ran the following to know the default number of lines :
blastn -help -max_target_seqs
Number of database sequences to show one-line descriptions for

Not applicable for outfmt > 4 Default = 500'     * Incompatible with:  max_target_seqs  -num_alignments <Integer, >=0>    Number of database sequences to show alignments for    Default = 250' * Incompatible with: max_target_seqs

I suggest you have a look at the settings for the blast version you are using. These can usually be accessed with 'blastn -help'. In particular, '-max_target_seqs' may be of interest.

Please [edit] your question and i) remove the images of the files, that isn't telling us anything. You can show us the output of wc -l  file to show the number of lines. ii) show us the exact commands you used. You are probably using a version of blast that shows the top 500 hits by default, but we need to know what blast you use and what exact commands you ran in order to help.

@Laura, Thank you for your suggestions. When I did run the code provided by you, this is what I got :
Number of database sequences to show one-line descriptions for
Not applicable for outfmt > 4
Default = 500'     * Incompatible with:  max_target_seqs  -num_alignments <Integer, >=0>    Number of database sequences to show alignments for    Default = 250'

Incompatible with:  max_target_seqs

How can I change these settings?

Not an expert on 16S, but is BLAST really the correct tool for this(?) Maybe one of these tool suites like QIME

QIIME is an interesting suggestion. Here I'd say BLAST - although there are performance issues - offset parallelising. BLAST is a known approach to scoop all the genetic diversity for a given locus and 16S is a taxonomic goldstandard. This isn't the analysis - the analytics are in the next step, but how you've prepped the data for the next step needs to be understood and be solid. There are a lot of other approaches but blast has flexibility e.g. to shift tblastn .. etc ... and its why I like it.

Personally, IMO "modulisation" custom data pipelines in general, i.e. loads of steps - gives more control and each developer specialises on their component rather than try and do everything at once. The latter is the QIIME approach, but its an interesting debate.

I've reproduced the issue. The solution is to simply set -max_target_seqs to whatever number of hits you want. E.g. -max_target_seqs 1000 gives 1000 and so on ... The problem  is specific to the -outfmt option and not default searches.
blastn version  I'm using v2.12.0, Feb  8 2022

Proof
My database contains around 260000 fasta sequences.
makeblastdb -in mystuff.aln.0.fa   -parse_seqids -blastdb_version 5  -out test.db -title "test" -dbtype nucl
blastn -query query.fas -db test.db -evalue 1e-6 -num_threads 4 -out blasthits.txt

Then wc blasthits.txt

>5000 hits

If we now use -outfmt 6, but remove the stringency of -evalue 1e-6 then ...
blastn -query query.fa -db test.db -outfmt 6 -num_threads 4 -out blasthits.txt

Now perform wc blasthits.txt

500

So the OP is correct, I didn't know -outfmt 6 restricts the output to 500.
If we now introduce the max_target_seqs 1000 parameter, suggested by @Laura then,
blastn -query query.fa -db test.db -outfmt 6 -max_target_seqs 1000 -num_threads 4 -out blasthits.txt

Then wc blasthits.txt

1000

If its 550 wc ... is 550 etc...
Thats it. Note in my blasts I usually parallelise across 4 processes i.e. -num_threads 4.

Thank you very much, it worked!

Thanks @KirtiKulanthaivel glad it worked and much appreciated.


text:
Generating 3D cartesian coordinates of a peptide and computing van der Waal interactions (1-4 atoms) for a given set of dihedral angles?
I wanted to generate 3D Cartesian coordinates of a peptide sequence for fixed bond lengths and bond angles for a given set of backbone dihedral angles and side dihedral angles. I also wanted to compute non-bonded van der Waal interactions (1-4 atoms) for each of the peptide generated from the dihedral angles. I initially tried frag builder module in python although it works fine. There are few issues,

I could not fix bond length and bond angle of my own choice.

It generates PDB file and save it inside the hard disk. So I need to read it again and parse it to compute non-bonded van der Waal interactions (1-4 atom). It consumes lot of time

Also I find fragbuilder is not that fast in generating peptides.

Reading , parsing and computing the interactions takes around 30 mins to complete the task for 200,000 3D ensemble of structures of given peptide sequence of length 18.

Kindly suggest me a better module in python which does the job faster and also give reference to it.

Don't ask multiple questions in a single post. Focus on solving one problem at a time.

@M__    Yes Fragbuilder has many shortcomings as pointed out in the questing itself. So I am looking for better alternative. Yes, I am limited to python. Is there any other modules which could handle this inside python ?  Yet, I want to know whether there is any other method outside python as well. Thanks in advance !

Thanks @vigneshwarankannan updating the question would be cool.

The PDB module in Biopython is worth checking out for obtaining angles, distances, torsion angles. It is described from p208 at http://biopython.org/DIST/docs/tutorial/Tutorial.pdf. It doesn't do van der Waal interactions.
For example,
from Bio.PDB.PDBParser import PDBParser
from Bio.PDB.Chain import Chain
from Bio.PDB.internal_coords import *
from Bio.PDB.PICIO import write_PIC, read_PIC, read_PIC_seq
from Bio.PDB.ic_rebuild import write_PDB, IC_duplicate, structure_rebuild_test
from Bio.PDB.SCADIO import write_SCAD
from Bio.PDB.PDBIO import PDBIO

parser = PDBParser(PERMISSIVE=1)
myProtein = parser.get_structure("1a8o", "1A8O.pdb")
myChain = myProtein[0]["A"]
IC_Chain.MaxPeptideBond = 4.0
myChain.internal_coord = None
myChain.atom_to_internal_coordinates(verbose=True)
r2 = myChain.child_list[1]
r2ic = r2.internal_coord
r2chi1 = r2ic.get_angle("chi1")

I wanted to compute the 3D structures of the peptides (i.e) x,y,z coordinates of the atoms for the given set of input Ramachandran angles for a fixed bond lengths and bond angles. I don't require to compute the Ramachandran angles of a given PDB file.

@vigneshwarankannan thanks for the feedback. What you want is simply to extend the Biopython module with appropriate code for the trigonometric calculation, i.e. angles and distances. It's not my area - but it would be doable.


text:
editing dataframe in r
One of the columns in my R data frame looks like this:
[1] "SVLEN=1007"   "SVLEN=."      "SVLEN=."      "SVLEN=."      "SVLEN=204"    "SVLEN=-379"
 [7] "SVLEN=-492"   "SVLEN=-825"   "SVLEN=-713"   "SVLEN=."      "SVLEN=3553"   "SVLEN=432391"
[13] "SVLEN=."      "SVLEN=-52"    "SVLEN=>1615"  "SVLEN=."      "SVLEN=6899"   "SVLEN=."
[19] "SVLEN=."      "SVLEN=4468"   "SVLEN=."      "SVLEN=."      "SVLEN=."      "SVLEN=627"
[25] "SVLEN=."      "SVLEN=."      "SVLEN=."      "SVLEN=1171"   "SVLEN=."      "SVLEN=625"
[31] "SVLEN=."      "SVLEN=."      "SVLEN=-237"   "SVLEN=."      "SVLEN=."      "SVLEN=-30"
[37] "SVLEN=39"     "SVLEN=."      "SVLEN=-30"    "SVLEN=-377"   "SVLEN=120497" "SVLEN=."

and I would like it to look like this:
[1] "1007"   "."      "."      "."      "204"    "-379"
 [7] "-492"   "-825"   "-713"   "."      "3553"   "432391"
[13] "."      "-52"    ">1615"  "."      "6899"   "."
[19] "."      "4468"   "."      "."      "."      "627"
[25] "."      "."      "."      "1171"   "."      "625"
[31] "."      "."      "-237"   "."      "."      "-30"
[37] "39"     "."      "-30"    "-377"   "120497" "."

I could go in to the excel file I got it from and edit it but is there an easy way to do this in R?
ANSWER:
I found a solution to this probelm going back to the R documentation (https://www.statology.org/gsub-r/) use gsub to section off characters in a string. My command for this was: gsub(".*=","\\1",df$SVLEN_INFO)
And that worked!

Instead of adding that extra section to the question, can you please create an answer and put it there? It'll help the StackExchange algorithms to do the right things.

Can you check the \\1? I'm not an R expert but you need brackets for \\1 to work, otherwise just use '' as the replace term not with brackets. Thus it should be gsub(".*=(.*)","\\1",df$SVLEN_INFO). Alternatively gsub(".*=(-?\d+)","\\1",df$SVLEN_INFO) might work better, not sure. I don't know how type casting works in R. @RamRS has a simple solution.

Adding an answer for now:
You can use sub to remove the SVLEN= part. You don't need gsub and you don't need regular expressions, let alone capture groups.
df$SVLEN_INFO <- sub(x = df$SVLEN_INFO, pattern = "SVLEN=", replacement = "")

This is a good answer @rimo

this worked great! thank you

Following gringer's advice I am also adding this solution:
Used gsub to section off characters in a string.
My command for this was: gsub(".*=","\\1",df$SVLEN_INFO)
Thanks for all the help!

Does this work? You're using a backreference (\1) without a capture group (like (.*)), which usually does not work. Also, please use code formatting to format your post properly, I've done it for you this time.

I checked, this only works because "\\1" is empty since you don't have a capture group. Just use "" instead of "\\1". Otherwise, this is NOT a good solution as you're essentially using a bug.

worked for me at the time but again the answer above worked as well...so I'll use that to format moving forward...

Worked for you because it's a bug. If they fix the bug, it will stop working. Many reg-ex based tools error out or behave differently if you use backreferences without capture groups. IIRC sed retains the whole matched expression as the default backreference. Use predictable tools predictably instead of depending on bugs to hack your way through problems.


text:
How to get the nucleotide sequence for protein NCBI accession code, programatically?
It is a lot easier in my work to do blast searches against protein databases when I can, since I work with distantly related organisms. But sometimes, I need to work with the actual nucleotide sequences. Until now, I've often done it manually, but it is getting unfeasible due to increasing data amounts.
The NCBI website is structured such that I suspect automating the process (e.g. input a list of protein accession numbers, get out a fasta file with the corresponding nucleotide seqs) should be possible. At the same time, I am not quite sure what tools to use (biopython?), and if a tool already exists, I would rather not reinvent the wheel. Would you happen to know of anything like this? Thank you for your time!
Update and WIP:
I have tried Entrez Utils. In part it works, but I am having some issues scaling it up. For example,
elink -db protein -target nuccore -id WP_244350345.1 | efetch -format fasta_cds_na | awk '/^>/ {printf("\n%s\n",$0);next; } { printf("%s",$0);}  END {printf("\n");}' | grep WP_244350345.1 -A 1
works. (elinks to connect the protein accession WP_244350345.1 with the nucleotide database | efetch to grab the nucleotide sequence, which is a bit of a pain because it grabs the entire genome rather than the cds I want | so I need a bit of cleaning, including using awk to turn the resulting multi-line fasta to single-line fasta | and grep to actually get the cds sequence associated with WP_244350345.1. As I said, it works, and it results in:
>lcl|NZ_JAFIRA010000020.1_cds_WP_244350345.1_35 [gene=bchM] [locus_tag=JX360_RS09120] [protein=magnesium protoporphyrin IX methyltransferase] [protein_id=WP_244350345.1] [location=35695..36360] [gbkey=CDS]
ATGCAAGAGAAAAGCATCGTCCGCGACTATTTCAACTCCATTGGCTTCGACCGTTGGCGGCGCATCTACGGAGAGGATGAGGTCAATTTTGTCCAAAAAGACATTCGCTCCGGCCATGCCCGTACCGTAGCCACTGTGTTGGAATGGCTGGGGGATCCCTCAGGTTTGCAGATTTGTGATGCGGGGTGTGGGGTCGGTAGTCTTAGTTTGCCGTTGGCTGCCAAAGGGGCCCGCGTCTTCGCCAGCGACATTTCCGAACAAATGGTGAACGAAGCCCGCCGTCGCCAGCAAATCCAGTTGGGATCCACCGAGAACCCCCAATTTCGGGTCTCGGACTTAGAGGAACTAAGCGGAGAATATGATGCCGTCATTTGCCTGGATGTGATGATCCACTACCCAGAAGCCGATGCCCTGCGCATGCTGGAACATCTAACCCGCCTGGCCCGCTCTCGCTTGATCTTCAGTTTTGCCCCCAAAACCCCTCTGTTGACGTTGCTAAAAAAAGTGGGAGAGTTCTTCCCCGGCCCCAGCAAAGCCACCCGTGCCTACCAACACCGAGAAGCGATTTTGGTGAGCAAATTAACCGAATTGGGCTGGAACGTTCAGCAACGACAGACCATTCGTAGCCGCTTTTACTTTGCCTGGATTCTGGATTTGGCCCGGTGA

However, trying to read multiple accession numbers from a file, in a loop, is not working. My current code (essentially the same as before, just replacing WP_244350345.1 with ${LINE}:
while IFS='' read -r LINE || [ -n "${LINE}" ]; do
echo "${LINE}"
elink -db protein -target nuccore -id "${LINE}" | efetch -format fasta_cds_na | awk '/^>/ {printf("\n%s\n",$0);next; } { printf("%s",$0);}  END {printf("\n");}' | grep "${LINE}" -A 1
done < "results.txt" | more

The echo "${LINE}" part works fine, the rest does not. It has also worked perfectly fine with other parts of Entrez Utils (see below), so I am confused and I would appreciate some help.
while IFS='' read -r LINE || [ -n "${LINE}" ]; do
organism=$(elink -db protein -target taxonomy -id "${LINE}" | efetch | tr --delete '\n' | sed 's/    species, /_/g'| sed 's/^.*\. //g' | sed 's/ /_/g')
echo "${LINE}"
echo "$organism"
done < results.txt | more

results in:
WP_244350345.1
Thermostichus_vulcanus_cyanobacteria
WP_073547903.1
Chroogloeocystis_siderophila_cyanobacteria
WP_088429095.1
Halomicronema_hongdechloris_cyanobacteria
...

Can you please add any solutions you find as an answer instead of as additions to questions? This keeps the Stack Exchange automated algorithms happy, so they don't try to get more answers for an already-answered question.

Yes its easy using Biopython.
What I'm requesting is the sequence must be >9500 otherwise it's discarded. In your case this value will be MUCH bigger, because these are bacterial genomes.
from Bio import Entrez
from Bio import SeqIO, AlignIO

def getgenbank (smallGenDict):
    bacteriadb = []
    for bacteria in smallGenDict.keys():
        Entrez.email = "name@email.com"
        with Entrez.efetch(
                db="nucleotide", rettype="gb", retmode="text", id=bacteria
                ) as handle:
            seq_record = SeqIO.read(handle, "gb")
            # record = get_CDS(seq_record, bacteria, smallGenDict.keys)
            record = seq_record # better to remove the above function because its outside the question
            if len(record.seq) > 9500:
                bacteriadb.append(record)
    return bacteriadb

def get_CDS(record, bacteria, smallGenDict):
    for feature in record.features:
        # custom manipulation of input sequences
        # better to forget this bit of code
        # this whole part of the code has been skipped
        # its capturing the coding sequence (not relevant)
return record

def main(smallGenDict):
    db = getgenbank(smallGenDict)
    SeqIO.write(db, '/Users/username/Desktop/myoutput.fasta', "fasta")

if __name__ == '__main__':
    smallGenDict = {'ABxxxxxx':'AGACTGATA', 'ABxxxxx2:'GAGATGA'} # better to automate this bit (different question)
    main(smallGenDict)

Thats it.
Notes In the code below, the file smallGenDict is a dictionary of Genbank accession No.s and their sequence.
{'ABxxxxxx':'AGACTGATA', 'ABxxxxx2:'GAGATGA'}
This better written here as a list
['ABxxxxxx','ABxxxxx2' ...] say call it smallGenList.
Where it says smallGenDict.keys(): just replace it with smallGenList and remove the method .keys(), otherwise the code will through a bug.
The function get_CDS is processing the sequence file to remove the non-coding regions. I've not included it because its specific stuff.

You can use R Rentrez, in the link you can find a tutorial for the package usage.

Thank you! I've tried Entrez (R is not something I am familiar with), but I got stuck.


text:
Create bed file where each column is the number of intersections with another bed file
For example I have 3 bed files:
window.bed:
chr1    0   1000
chr1    1000    2000

lp1.bed:
chr1    20  50
chr1    30  40
chr1    1500    1900

lp2.bed:
chr1    68  786
chr1    899 987

How can I create this output:
chr1    0   1000    2   2
chr1    1000    2000    1   0

In this output columns 1, 2, 3 are the standard bed format cols, of the chromosome, start, and end of the features taken from window.bed above. Column 4 is the number of intersections between window.bed and lp1.bed per feature in window.bed. For example, the first element in column 4 is '2' because there are two intersections between window.bed and lp1.bed. Column 5 is the same but for the number of intersections between window.bed and lp2.bed for each feature in window.bed.
In the case that a feature is a subset of another feature in either lp1.bed, or lp2.bed, they should be counted as two intersections. For example, the second feature of lp1.bed is a subset of the first feature. This is reflected in the first element of column 4 of the output.
I suspect that it is doable with bedtools intersect, but I cannot work out the syntax.

The relevant bedtools option seems to be -c:
    -c  For each entry in A, report the number of overlaps with B.
        - Reports 0 for A entries that have no overlap with B.
        - Overlaps restricted by -f, -F, -r, and -s.

We can get the number of overlaps per region per file with:
$ bedtools intersect -c -a window.bed -b lp1.bed
chr1    0   1000    2
chr1    1000    2000    1

$ bedtools intersect -c -a window.bed -b lp2.bed
chr1    0   1000    2
chr1    1000    2000    0

So we're almost there, we just need to combine them, and that's a job for paste:

$ paste  <(bedtools intersect -c -a window.bed -b lp2.bed) \
         <(bedtools intersect -c -a window.bed -b lp1.bed)
chr1    0   1000    2   chr1    0   1000    2
chr1    1000    2000    0   chr1    1000    2000    1

But that gives us extra fields, so we can use cut to remove them:
$ paste  <(bedtools intersect -c -a window.bed -b lp1.bed) \
         <(bedtools intersect -c -a window.bed -b lp2.bed | cut -f 4)
chr1    0   1000    2   2
chr1    1000    2000    1   0

The paste command just joins lines from its two input files in the order it receives them, but the bedtools commands will always output the same order as the lines in window.bed so this should be fine.
The <(command) construct is called "process substitution" and is a way of passing the output of a command as a file to another command that expects files as input. Note that it isn't supported by all shells, but if you're using bash or zsh which are the most common these days, it should work.

Via bedmap --count, which gives the option to specify a custom value when regions do not map:
paste <(bedmap --echo --count --unmapped-val 0 --delim '\t' window.bed lp1.bed) <(bedmap --count --unmapped-val 0 window.bed lp2.bed) > answer.bed

Ref. https://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html#element-and-overlap-statistics


text:
Unable to install an R Bioconductor package on Apple M1/M2
I have an Apple latest arm64 machine. Running the latest R 4.3.0 version. However, I'm unable to install an R package - S4Arrays.

Doesn't look like there is an arm64 build for this package. Anybody know how to resolve it?

Can you please paste your error message as text into your question? Images of text are difficult to read, and harder for screen readers to understand.

The easiest solution is to install the x86 version of R in your M1/M2 Mac here.
It will be then emulated by the Rosetta. After that you will have no problems with Bioconductor packages

This is a huge pain once one goes down the ARM path for R and libs. I had to revert back when BioC did not support ARM but now, there are other hacks.

Yes, I agree it is a pain! However R 4.3 is out for a couple of week, forcing you anyway to start from scratch with the packages' installation, in this situation I think it'll be easier this solution than the various hacks around


text:
Creating a scatterplot by a json data in js
I have a json file of protein data of mouse and human like
{
    "name": "P04202",
    "full_name": "Transforming growth factor beta-1 proprotein",
    "symbol": "TGFB1_MOUSE",
    "length": 390,
    "species": "mouse",
    "function": "Transforming growth factor beta-1 proprotein: Precursor of the Latency-associated peptide (LAP) and Transforming growth factor beta-1 (TGF-beta-1) chains, which constitute the regulatory and active subunit of TGF-beta-1, respectively.",
    "nr_modifications": 2821,
    "residues": [
        {
            "pos": 1,
            "residue": "M",
            "modifications": []
        },
        {
            "pos": 2,
            "residue": "P",
            "modifications": []
        },
        ...
        {
            "pos": 60,
            "residue": "A",
            "modifications": [
                {
                    "mod": "[127] Fluoro",
                    "type": "Chemical derivative",
                    "prob": 1
                },
                {
                    "mod": "[53] HNE",
                    "type": "Post-translational",
                    "prob": 0.25
                }
            ]
        },
        ...
        ]
}

How do I create a scatter plot [![enter image description here][1]][1]
Thanks for any idea
So far I gathered a code to fetch the data from an online json
export const load = ({ fetch }) => {
  const fetchProteins = async () => {
    const res = await fetch('https:address/assets/ptm.json')
    const data = await res.json()
    return data
  }

  return {
    Proteins: fetchProteins()
  }
}

This code generates a scatter plot but I will need to hover around the dots
script>
    import { extent } from 'd3-array';
    import { scaleLinear } from 'd3-scale';
    export let data = [];

    const margin = 30;

    $: xScale = scaleLinear()
        .domain(extent(data.proteins.map((d) => d.length)))
        .range([0, 800 - (margin * 2)]);

    $: yScale = scaleLinear()
        .domain(extent(data.proteins.map((d) => d.nr_modifications)))
        .range([0, 400 - (margin * 2)]);
</script>

<svg width="800" height="400">
    <g transform="translate({margin}, {margin})">
        {#each data.proteins as datapoint}
            <circle
                cx={xScale(datapoint.length)}
                cy={yScale(datapoint.nr_modifications)}
                class:human={datapoint.species === 'human'}
                r="20"
            />
        {/each}
        <circle
                cx={10}
                cy={10}
                class="human"
                r="10"
            />
            <circle
            cx={10}
            cy={40}
            class="mouse"
            r="10"
        />

    </g>
</svg>

What did you try so far?

I gathered a code to fetch the data from an online json

"fetch data from online JSON" is basically downloading a file and maybe subsetting it. What else have you tried in actually generating the plot?

Yesterday I accomplished to plot a scatter chart but I will need to deploy hover around the dots

It seems like you are already on the right track. The code you provided fetches the protein data from the JSON file and generates a scatter plot with D3.js based on two properties: length and nr_modifications.
To add an interactive hover event, you can add a tooltip to your scatterplot. The tooltip will display additional information when you hover over each dot.
Below is an example of how you can add a tooltip:
First, add a
 element in your HTML that will be used as the tooltip:

Protein:

Modifications:

Length:

Then, in your script, add the code to handle the mouseover and mouseout events for each circle:
{#each data.proteins as datapoint}
    <circle
        cx={xScale(datapoint.length)}
        cy={yScale(datapoint.nr_modifications)}
        class:human={datapoint.species === 'human'}
        r="20"
        on:mouseover="{(e) => handleMouseover(e, datapoint)}"
        on:mouseout="{(e) => handleMouseout(e)}"
    />
{/each}

Define the handleMouseover and handleMouseout functions in your script. These functions will show and hide the tooltip respectively. For the handleMouseover function, it also updates the tooltip content with the protein data:

In standard HTML or JavaScript, the {#each} syntax and usage of {} for expressions and variable binding won't work as they are specific to Svelte. Same applies to the usage of on:mouseover and on:mouseout. In standard HTML or JavaScript, we use onmouseover and onmouseout.
Here are some things which I haven't checked, but you might want to consider looking at:

It looks like you are trying to load JSON data from an external
source, but the function fetchProteins() is never called and its
result (the promise returned) is not being awaited. The way it's
structured now, Proteins in the returned object will hold a Promise,
not the actual data.
The code assumes that the returned JSON data is an array stored in
data.proteins, but the provided JSON example suggests that the data
structure might be different. It doesn't show any 'proteins' key in
the JSON. The code needs to align with the actual structure of the
returned JSON data. data.proteins implies that data is an object and proteins is a property of that object containing an array.
In your SVG transform attribute, you are using braces {} though the correct syntax for a
translate transformation in SVG is translate(x,y).

You might want to try looking at the following changes:
export const load = ({ fetch }) => {
  const fetchProteins = async () => {
    const res = await fetch('https:address/assets/ptm.json')
    const data = await res.json()
    return data
  }

  return {
    Proteins: fetchProteins
  }
}

  $: xScale = scaleLinear()
        .domain(extent(data.map((d) => d.length)))
        .range([0, 800 - (margin * 2)]);

    $: yScale = scaleLinear()
        .domain(extent(data.map((d) => d.nr_modifications)))
        .range([0, 400 - (margin * 2)]);

<g transform="translate({margin}px, {margin}px)">
        {#each data as datapoint}
            <circle
                cx={xScale(datapoint.length)}
                cy={yScale(datapoint.nr_modifications)}
                class="{datapoint.species === 'human' ? 'human' : ''}"
                r="20"
            />
        {/each}

</g>