text:
Can Bioinformatics be used in developing a marketable product?
As we all know computer science, electrical engineering, etc. knowhows can be used in developing software and electronics/electrical products that can be sold to end-users (B2C).
Can bioinformatics be used to develop such products other than drug discovery?
Or, is it just about scientific/academic research meant for B2B clients?

This question appears to be asking for opinions, rather than an answer to a specific problem, and as such doesn't fit well with the StackExchange format. I recommend you try elsewhere for an answer, e.g. the Bioinformatics subreddit.

It is pivotal to vaccine research at multiple levels:

Deep theoretical vaccine design, which takes years to get to market (my stuff's gone into that and still stuck in advanced clinical trials).
Looking for variants that might have evolved immunoevasion against current vaccines involving SARS-CoV-2 (why do we need so many boosters?) and influenza.

Bioinformatic patents are tricky.
The reality is both points are not systematically automated in 'point and click' user friendliness. Its definitely of commercial value at an industrial/institutional level. However most work is submitted to the public domain and code freely available or held under a commercial firewall.
Generally, I think you underestimate how much bioinformatic work takes place in modern clinical practice and when fighting a virus how it just keeps evolving around all the obstacles thrown at it.
More generally, a lot of bioinformatics is used to direct wet-lab research because bioinformatics gets to the answers quicker.


text:
Understanding protein-peptide docking analysis
I want to explore some potential reactive peptides and seems a common step to conduct some form of docking analysis. I have a question about the theory and a practical question.
Practical question:
An interesting secondary step seems to be to calculate the binding energy with something like PRODIGY, although I'm struggling to understand the .pdb files to run the webserver.
Once I have created a .pdb on either HPEPDOCK or CABS-dock, how do I find the naming convention for the peptides ( Ligand ID ) within the .pdb file? So that I can run in Prodigy-Ligand.
Theory question:
I'm struggling to understand what to take-away from the docking analysis.
I've been using HPEPDOCK ( & CABS-dock, although the former is significantly quicker and seems a good starting place ).
I've been trying to read some papers and see how they use docking analysis. I get the feeling, particularly for peptide-protein ( known receptor ), it is simply a visual inspection to see if it makes sense. Adding the abstract docking score to double check it is low enough ( < -100 ).

I'm struggling to understand further use - most papers seem to simply state they did such analysis and ~ it looked okay, or ~ chose the good ones.
If anyone can recommend some literature that might help me understand, that would be great!

'Reactive peptides' (covalent peptidomimetics) are generally far from a medchemist's first choice.
A ligand is not a peptide, it is a small molecule. A residue in compbiochem is a unit and can be AAs residues, NCAAs residues, ligands or ions. Caffeine has the residue name/id CFF. Therefore, you need to use the first tab: protein-protein.

source: Stack Exchange

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

sort: votes


text:
Are the GIAB reference genomes (HG00N) fully phased?
Specifically, the NISTv4.2.1 HG002 dataset here (more info on Genome in a Bottle). According to their README, the HG002 v4.2.1 VCF should be fully phased.

For HG001, HG002, and HG005, we have phased the v4.2.1 benchmark vcfs by transferring phasing information from a phased assembly-based trio-hifiasm dipcall vcf to the v4.2.1 vcf

I ask because there appear to be FAR more heterozygous variants on the second haplotype than the first (0|1 >> 1|0).
For example:
tail -n 500000 HG002_GRCh38_1_22_v4.2.1_benchmark.vcf | \
        cut -f 10 | cut -f1 -d: | sort | uniq -c
 302960 0/1
   1104 1/0
 190042 1/1
   2037 1/2
   3857 2/1

Is this simply an artifact of how these truth sets were produced? For example, if the GRCh38 (haploid) FASTA was meant to be haplotype 1, then it makes sense that it's nearly identical to HG002's haplotype 1, and haplotype 2 has more differences. Is that the case, or are only certain variants phased (such as those that pass a minimum confidence threshold, or are in a high-quality BED region)?

As the two alleles in a genotype are separated by / but not |, the main v4.2.1 VCF is not phased. The README says the phased VCFs are available in the SupplementaryFiles directory. It contains a few versions of VCFs with v4.2.1 genotype calls but with phasing from additional tools such as whatshap and hifiasm. I don't know which version is the most accurate, though.

Thanks, I was aware of the convention that / means unphased and | means phased, but the fact that it labeled some variants as 1/0 really threw me off. Not sure why they'd do that.

I had interpreted the README to mean that they included the trio hifiasm dipcall VCFs in the SupplementaryFiles directory, not the resulting phased reference benchmark. Thanks!

@TimD1 I asked a question about that recently and got an interesting possible reason, but I haven't had the time to look more closely. Here it is, in case you find it interesting: Why does GATK produce both 0/1 and 1/0 genotypes in the same file? Are the two not equivalent?

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


text:
deseq2 full and reduced model interpretation
I would like to know how to interpret the output of the formula although i been thorough loads of literature but I'm not confident yet.
So this is my full model I'm running
dds <- DESeqDataSetFromMatrix(countData=rpkm_ordered, colData=coldata, design= ~ FAB + Sex + Age + TMB + WBC + BM_percentage)

So if I have to interpret the output of this how it would be ?
What I have read is this
Here in the above formula we are trying or controlling the FAB variable while we look for variation that is coming from other variables?
Now my confusion is when I made this dds object the FAB was the variable which I fixed while I look for variation coming from other variable.
Next thing is I would like to run a reduced model.
I tried this
test.Age <- as_tibble(results(DESeq(dds, test="LRT", reduced= ~Sex + Age )))

Now my question is since I'm using the same dds object here where the FAB is the variable i fixed and  with new formula which is
reduced= ~Sex + Age this so how would i interpret this model?
Any suggestion or help would be really appreciated

You know that RPKM is inappropriate for DESeq2?

yes..are you rpkm_ordered telling about this i suppose its just dataframe i named in my code which i didn't change..

@swbarnes2 sorry for this  rpkm_ordered  confusion since once i had to reorder my dataframe...which I did quite long time back so its just in my code..

By fitting the full model ~FAB + Sex + Age + TMB + WBC + BM_percentage you can estimate the effect of FAB (or any of the other variables) on RNA expression after adjusting for differences caused by Sex, Age, TMB, WBC, and BM percent. See this post for exactly what it means to adjust for other variables
When you fit a reduced model (in your case ~Sex + Age) and then conduct a likelihood ratio test (LRT), you are comparing the fit of the full model (~FAB + Sex + Age + TMB + WBC + BM_percentage) to the fit of the reduced model (~Sex + Age). If the fit is substantially better, the likelihood of the full model will be larger than the reduced model and the p-value will be significant. This means that adding the extra variables that are in the full model (FAB, TMB, WBC, BM_percentage) improved the fit of the model to the data because together these variables have some type of relationship with RNA expression even after adjusting for the variables in the reduced model.
The challenge with using a LRT is that it doesn't tell you which of the variables added in the full model are responsible for the improved fit - it just indicates that the full model is better.  This type of test is useful when you are trying to make broad statements about several variables having a relationship with RNA expression. From the DESeq2 vignette:

The LRT is therefore useful for testing multiple terms at once, for example testing 3 or more levels of a factor at once, or all interactions between two variables. The LRT for count data is conceptually similar to an analysis of variance (ANOVA) calculation in linear regression, except that in the case of the Negative Binomial GLM, we use an analysis of deviance (ANODEV), where the deviance captures the difference in likelihood between a full and a reduced model.

Definitely read the DESeq2 vignette for more details. In your case if the LR test is significant, it means that FAB, TMB, WBC, and BM percent combined are associated with RNA expression.
Most people are interested in specific comparisons when performing differential expression analyses with DESeq2. These comparisons are usually like do males have higher gene expression than females (comparing sex differences), do patients with higher TMB have greater expression than patients with lower TMB (comparing effect of TMB on RNA expression). If you want to do a specific comparison, you will likely want to use the Wald test, which calculates p-values for individual coefficients (male vs female coefficient, TMB coefficient etc). The Wald test determines whether a coefficient is equal to 0 or some other null hypothesis (see lfcThreshold parameter).
The Wald statistic and accompanying p-values are reported by default when using results() and not specifying test="LRT"
There are some additional differences between Wald tests vs LRT, including accuracy differences when operating on small sample sizes and computational efficiency issues when operating with large sample sizes. Discussion on Wald vs LRT
TL;DR Use LRT to test if several variables or factor levels are associated with RNA expression. Use the Wald test to determine if specific comparisons are significant.

thank you for the detailed write up...this is more than enough for my data reporting which would be now be more statistical rather than using loose terms thanks a ton

Little doubt "The design formula should have all of the factors in your metadata that account for major sources of variation in your data. The last factor entered in the formula should be the condition of interest." are these saying the same thing you described?  https://hbctraining.github.io/DGE_workshop_salmon_online/lessons/04a_design_formulas.html

The last factor in the formula is the variable whose comparison (if sex is the last variable then male vs female) will show up by default when you call results(). You can change which comparison you get with the contrast argument.

if I have a single factor in my design for example age and i ran deseq2 what would be my interpretation  would I be saying here Im testing the effect of age on the differences in gene expression?is that statistically correct statement?

Yes that is correct. If the age variable is continuous, the age coefficient for a gene $g$ would represent the log2FoldChange of $g$ when age increases by 1 unit. If age has no effect on gene expression for $g$, the coefficient will be very small and not significant (assuming adequate power)

in case of numerical such as age what becomes my intercept? for categorical Im now clear ,in case of age is it the lowest age becomes the intercept ?

The intercept would be the value when age is 0. Based on your questions it would probably be helpful to learn more about linear models as they are very common in bioinformatics. DESeq2 uses a negative binomial linear model. I would start with linear regression and then checkout this book

'it would probably be helpful to learn more about linear models as they are very common in bioinformatics' yes i been vigorously digging and reading trying to wrap my head on the concept and the biology going on


text:
How to get strain names/ids contained in a multi FASTA file using seqkit?
FASTA files can be very big and unwieldy, especially if lines are at most 80 characters, one can't speed up browsing them by using less with -S to have one sequence every two lines.
How can I extract just the strain names (or sequence names, i.e. the string on the line starting with >) into a list? I'd like to use seqkit since I'm already familiar with it and know some of the flags and idioms it uses.
I googled the question but couldn't find an answer using seqkit - the examples all use awk, grep, sed or similar, which are alright for copy/pasting but not very memorable.

The nice thing about seqkit is that the answer is really simple and memorable:
seqkit seq -n sequences.fasta

This will output a text file containing one sequence name per line:
MPX/PT0001/2022
MPX/PT0002/2022
MPX/PT0003/2022

Thanks for posting back @CorneliusRoemer thats useful.

You can trivially do this with standard *nix tools. Here are a few options:
sed -n 's/^>//p' file.fa
awk 'sub(/^>/,"")' file.fa
grep '^>' file.fa | tr -d '>'
perl -ne 'print if s/^>//' file.fa

The advantage of these is that they will work on any POSIX-compliant system and most other *nix flavors, even if they're not POSIX such as Linux, out of the box without needing to install any external tools.
Another thing to take into account is how fast a command is, especially when dealing with large files. I ran a few tests using the solutions presented here on the human hg38 genome assembly which can be downloaded here and seqkit version 2.2.0 which was downloaded from here. The file was stored locally on my laptop's SSD drive, and the laptop itself is:
$ inxi
CPU: quad core Intel Core i7-6820HQ (-MT MCP-)
speed/min/max: 2896/800/3600 MHz Kernel: 5.17.5-arch1-1 x86_64
Up: 25d 10h 27m Mem: 15035.9/31958.1 MiB (47.0%)
Storage: 476.94 GiB (103.1% used) Procs: 452 Shell: Bash inxi: 3.3.15

To get more accurate results, I used a little function I have which will run the command you give it a specified number of times and then report the average real time across all runs. Here's the code if you want to try it:
timethis(){
    max=$1;
    shift;

    for com
    do
        printf '\nCOMMAND: %s\n' "$com"
    c=0;
    while [[ $c -lt $max ]]
    do
      let c++;
      ( time -p eval "$com" ) 2>&1 | grep -oP 'real.*?\K[\d\.]+'
        done |
        awk -vm=$max '{k+=$1}END{print (k/m)}';
    done
}

And the results. These are the average "real" (i.e. wall clock) times taken across 10 runs:
$ timethis 10 \
  "rg '^>' GCA_000001405.29_GRCh38.p14_genomic.fna | tr -d '>' >/dev/null" \
  "grep '^>' GCA_000001405.29_GRCh38.p14_genomic.fna | tr -d '>' >/dev/null" \
  "sed -n 's/^>//p' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null" \
  "perl -ne 'print if s/^>//' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null" \
  "./seqkit seq -n GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null" \
  "awk 'sub(/^>/,\"\")' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null" \
  "perl -nle 'push @res, $_  while s/^>//; END {print  join \"\n\",@res}' GCA_000001405.29_GRCh38.p14_genomic.fna"

COMMAND: rg '^>' GCA_000001405.29_GRCh38.p14_genomic.fna | tr -d '>' >/dev/null
0.474

COMMAND: grep '^>' GCA_000001405.29_GRCh38.p14_genomic.fna | tr -d '>' >/dev/null
0.583

COMMAND: sed -n 's/^>//p' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null
5.62

COMMAND: perl -ne 'print if s/^>//' GCA_000001405.29_GRCh38.p14_genomic.fna  >/dev/null
6.408

COMMAND: ./seqkit seq -n GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null
4.326

COMMAND: awk 'sub(/^>/,"")' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null
13.518

COMMAND: perl -nle 'push @res, 3  while s/^>//; END {print  join "\n",@res}' GCA_000001405.29_GRCh38.p14_genomic.fna
8.455

As you can see above, the grep or rg and tr approaches only took around half a second on average, so that is clearly much faster than any of the other choices. However, seqkit is pretty good since it only took around 4 seconds, faster than any of the other ad hoc solutions.

Seqkit is arguably a bioinfo standard tool ;)

@CorneliusRoemer I guess. I have never used it myself and have been doing bioinformatics for more than 20 years. But by "standard" I mean "standard *nix tools", not bioinformatics tools. I edited to clarify.

@CorneliusRoemer there are weakness on reliance on a developers 'kit'. Flexibility, if you wanted to count the number of pangotype A.1 in SARS-CoV-2, my code is easy to modify to />.A.1./ and introduce $count++ , whereas you would need to search seqkit to see if the developer had foresight. Bugs: are 'kits' completely bug free, how can you be sure? Upgrades: the promised upgrade ain't ever released (very, very common). Distributing your algorithm: this is a headache because you've now got a dependency which you don't have permission to redistribute.

@M__ of course there are use cases for *nix tools. But there are also use cases for *kits. In this case I wanted to use a memorable solution that will not require lookup every time. The two other answers are arguably off topic because they don't address the question - which explicitly asked for seqkit usage. Another *kit would have been more on-topic than *nix tools. I don't mind them, it's just a bit typical for SO, you ask a question and you get a solution that's for a different question ;)

@CorneliusRoemer answers providing different approaches are welcome across the Stack Exchange network (note that this site is not [SO]). Especially when there is already a fine answer with using the requested tool. As for what is memorable, that is entirely personal. You are familiar with seqt so that seems easy and memorable for you. I have been using *nix tools for 20 years, so nothing is more memorable than a simple sed s/foo/bar/. But that's fine! That's the beauty of it: we can all provide answers and people can choose whichever approach works well for them. Everybody wins!

Fair :) though it's seqkit not seqt

Heh, you see how memorable I find it, @CorneliusRoemer! :P

@CorneliusRoemer you've got to agree with terdon: seqkit is far too slow. It is worrying. This would cause problems on the data sets in question where a million strong fasta file is normal. Take your pick but realistically its the end of this discussion, there's no point using an inefficient algorithm against huge data sets, sooner or later they'll be a RAM issue.

That's a good reason to use an alternative tool, I agree, if speed matters. If you use the command line interactively though, terseness can be more important. It depends on the use case I suppose. I'd be curious how fast rg would be. The slowdown in seqkit seems to be related to file access? Large sys time.
Can you add a link on how to download the benchmark data so I can compare on my machine and try rg?

@CorneliusRoemer good catch! The test results I show were using a genome file stored on a network drive, which would of course affect the times. I had assumed that wouldn't be relevant since all tools were accessing the same file over the same network. However, I downloaded a different genome file locally and on my local SSD, seqkit is far more competitive. I'll update in a few minutes with the new test results, but do you have any idea why one tool would be affected by the network I/O times more than another? That's very surprising to me.

Ah interesting! I thought you were the bioinformatician with 20y of experience. I barely have 1 :D
Seqkit could read in chunks? No idea, one should maybe even open this as a bug on Github. By the way there are other similar tools out there.
Please try rg too, it's a Rust port of grep and likely faster than grep. Also best to link to the data file so it's reproducible. Performance could differ on macOS and *nix. Ah yes, that's an advantage of seqkit it could possibly work even on windows...

@CorneliusRoemer Ha! Serves me right for throwing my age around :) I do indeed have 20 years of experience, but I still come from the world of biology, not programming so there will be holes in my knowledge. I have updated the post with the new results, including all data needed to reproduce them, and seqkit now only loses out to grep + tr but is faster than all the others. So much for my years of experience, I have to eat my words now!

@M__ my initial numbers don't hold up when doing the same thing with a fasta file stored on a local SSD. So it looks like seqkit is slower than expected when reading data over a network but is faster than everything except the simple grep + tr when running on a local file.

@CorneliusRoemer if it is chunking it at source, its not a bug. It would make sense. Does it parallelise? Chunking and parallelising would be a technique to speed the operation and minimise RAM issues. Obviously it will fall over via a network.

Thanks @terdon its surely grep/ripgrep. If seqkit is chunking at source, that is really interesting (its clearly not pulling the whole file in and then chunking).

Or a bit more verbose ...
perl -nle 'push @res, $_  while s/^>//; END {print  join "\n",@res}' myfile.fa


text:
Publicly available population specific reference panels
The publicly available reference panel for GWAS analysis is the 1000 genomes reference panel. The reference panel consists of several populations. However, I am searching for a population-specific publicly available reference panel. I want to perform some analysis like Linkage disequilibrium, inheritance, etc.
Now there are some imputation reference panels for the specific populations such as UK10k (European), CAAPA (African), etc.
However, to access this panel we need to go through a lot of agreements and processes. Can anyone give me a list of publicly available population-specific reference panel, preferably, South Asian, European, Eastern Asian, African and Hispanic?

Can you describe why one of the populations in the 1000 genomes reference panel won't work for you?

are you telling me to subset the 1000 genome panel?

If that's your question then that's a great question. Answer to the question is that, I'm working on low coverage sequencing(<1x). To achieve the expected accuracy on the imputation of a specific population's sample, it's important for me to have samples of more than 1000 in the reference panel. Now, I can subset the 1000 geome panel. However, it consists roughly 400-800 samples per population which is pretty low for my cause. That's why I need reference panel with larger samples.

Very much all population or nation-specific large reference panels are park of national biobanks and you will not get the genotype-level data without having to do the process of requesting the data , since these are protected in almost all cases,


text:
Calculating the effect of differences between two groups of nucleotide sequences at every position
I'm investigating two different groups of short nucleotide sequences. The behaviour of a certain DNA-binding protein is likely influenced by common differences between these two groups leading high or low activity. Some positions are hypothesized to be more important than others. All sequences have the same length and are very similar.
I'm trying to find out how well changes at certain sequence positions explain the between-group difference and which positions are of the highest importance for this relationship. By looking at the sequence logo for both groups I can easily identify crucial positions. However, I'm at a loss how to actually score the effect each position has and how well it correlates to membership in either group. I thought about calculating Cramers V at each position but that doesn't seem to be very useful.
Can someone here point me in the right direction?

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.

are you talking about a sequence alignment of homologous sequences? Can you include a sample of what your data look like? Recall that the logo is based on sequence entropy, why not use that measurement?


text:
somatic SNV tool for ONT samples
I'm looking for a tool to call somatic single nucleotide variations on reads from an Oxford Nanopore Technologies (ONT) run.
I have both tumor and normal samples and have already aligned them to a reference genome.
Popular tools like mutect2 have been developed for short reads, and according to this paper it is not recommended to use tools specifically designed for short reads with long reads.
I know of tools that work with ONT reads like clair3 which is used in the epi2melabs workflow for SNPs wf-human-snp,  but clair3 takes only one sample at a time, so seems for germline calling only.
Another approach would be to simply call the SNVs in both tumor and normal samples independently with clair3, and then subtract normal SNVs from tumor SNVs, but according to what I have read in a few papers and online resources like this from gatk, this is a "naive approach", and a somatic caller should use both normal and tumor samples together.

Do you know of any tool that takes both tumor and normal long read samples to detect somatic variants (similar to mutect2 but for long reads)

Have you used tools for short reads like mutect2 with long reads?

Have you attempted the "naive approach" of simply subtracting SNVs?

Thanks,
Walter.

Note that tumor/normal is only one way of doing somatic variant calling. It is unquestionably the best way, but in the clinical space, it isn't as common as you'd expect and in many cases you only have the somatic sample without the normal control. The main reason you don't want to use standard variant callers is that you need to be able to call variants with very low allelic balance (proportion of reads supporting the variant) because of the heterogeneity of somatic samples which will be a mix of tumor and normal cells, as well as the heterogeneity of the tumor cells themselves.

ClairS is now released for long-read somatic variant calling. https://github.com/HKU-BAL/Clairs

Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.


text:
How to make pan-core genome curve through command line on linux?
I´m working with a dataset of 566 genomes to analyze a pangenome. So I was working with PANWEB to create this pan core genome curve, however, there is too much sequence to work with this webserver. Well, specifically I´m looking for this kind of plot:

So, I just want to ask for some ideas of programs to do that with a big dataset and with a command line.
Here is the reference: https://www.researchgate.net/figure/Pan-genome-and-core-genome-Graph-representing-the-pan-genome-blue-and-core-genome_fig1_317390138

The command create_pan_genome_plots from Roary can create a plot like this:


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

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

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

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

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

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


text:
How do I find unique file names without their extensions in shell script?
I need to generate a list of unique file names without their extensions or using the shell script command "find"
Say I have a list such as:
file2022_col.R1.twx.gz
file2022_col.R2.twx.gz
file2022_col.bwa.bam
file2332_col.R1.twx.gz
file2332_col.bwa.bam
file3442_col.bwa.bam.bai

I want a list of unique file names such as:
file2022_col
file2332_col
file3442_col

I'm trying to use find -type f -name but I'm not sure how to construct the option to restrict the find command.

It is simpler to use cut command to remove filename extensions:
ls -1 | cut -d '.' -f1 | sort | uniq

This seems nice and clean, but it will also pick up any directory names. I only want to extract from file names.

ls -1 -p | grep -v '/' could filter directories. The -p flag appends / indicator to directories. Then the -v flag only matches files.

ls | perl -p -n -e  's/(.*)(\..*){1,3}$/$1/' | uniq

or
find . -name "file*" | perl -p -n -e  's/(.*)(\..*){1,3}$/$1/' | uniq

find per se is a really powerful command, but its a while since I've used it full functionality.

Thanks, M__. I'll test it.  I notice you call perl.  I'm hoping to find a way to run it entirely in shell, but this could work too.

Bash has a regex, so it can be rewritten, but dunno if it has the {1,3} syntax

Use find, remove all extensions with perl, and make unique with sort:
find . -type f -printf '%P\n' | perl -pe 's{[.].*}{}' | sort -u

The Perl one-liner uses these command line flags:
-e : Tells Perl to look for code in-line, instead of in a file.
-p : Loop over the input one line at a time, assigning it to $_ by default. Add print $_ after each loop iteration.
s{[.].*}{} : replace [.] (a literal dot), followed by .* (any character repeated 0 or more times) with en empty string.
SEE ALSO:
perldoc perlrun: how to execute the Perl interpreter: command line switches
perldoc perlre: Perl regular expressions (regexes)
perldoc perlre: Perl regular expressions (regexes): Quantifiers; Character Classes and other Special Escapes; Assertions; Capture groups
perldoc perlrequick: Perl regular expressions quick start


text:
Single-cell RNA sequencing proper abbreviation
I wonder which form of a single-cell RNA sequencing abbreviation is correct:

scRNA-seq

scRNA-Seq

In the same single-cell sequencing Wikipedia article, both forms appear.
There is also inconsistency among papers I've seen - some authors prefer the first form, others the second one.
If there's any reason behind the different form, I'd like to hear about it.


text:
How to interpret a negative beta within a binary GWAS?
I have run a GWAS on a binary trait - "Are you following a carnivore diet?" - so the cases are individuals who follow this diet, and the controls are individuals who do not follow this diet. I have a genome-wide significant SNP, with a negative beta in reference to Allele 2. Does this imply that individuals with Allele 2 are less likely to follow this diet?

Can you provide some outputs and provide some details about exactly how you performed your GWAS and estimation? What is your statistical test, is there e.g. a logistic/logit transform involved?

Negative beta means that the trait under test is negatively correlated with the allele. What this actually means in relation to your specific question will depend on how the trait is coded and how the allele is coded.
It might, for example, be the case that "yes" is coded as 0 and "no" as 1, in which case a "no" answer would be associated with a high allele count for allele 2. You will need to check the algorithm to find that out, or (better) do your own manual / independent checking of alleles to discover the direction of association.
Take care when interpreting the results of these association tests, especially when deviating from statements beyond what the results actually say. You should also pay attention to is meaningless, whereas a beta value of -3.43 means there is a strong correlation between the trait and allele.


text:
Installation of Biovia DS visualiser on ubuntu
I am trying to install Biovia DS visualiser2021 on ubuntu 22.04. I used the procedure from the following website
https://askubuntu.com/questions/1175772/cant-install-discovery-studio-visualizer
When I run the installer ./install_DSClient.sh I got the following error
Calling LicensePack Installer

Verifying archive integrity... All good.

Uncompressing License Pack and extracting files, please wait...
./install_lp.sh: 73: Syntax error: redirection unexpected
Copying 1 files to ./share/license/data/tempdata.dat Please wait...
Copying 1 files to /home/user/BIOVIA/DiscoveryStudio2021/Logs/dsclientinstall.log Please wait...
Completed installing  DiscoveryStudio2021 at /home/user/BIOVIA/DiscoveryStudio2021

Check log file [ /home/user/BIOVIA/DiscoveryStudio2021/Logs/dsclientinstall.log ] for more information. Return code: 0

How to solve this error?


text:
snakemake tutorial stuck in basics-step-5 (variant calling)
I am learning to use snakemake. I am following the tutorial (https://snakemake.readthedocs.io/en/stable/tutorial/tutorial.html). In the section called "Basics", steps 1  to 4 worked well, processing one sample at a time, from mapping reads until indexing alignments.
snakemake --cores 2 sorted_reads/{A,B,C}.bam.bai

However, when I add the rule for joint variant calling (step 5):
snakemake --cores 2

I get the following error:
"WorkflowError:
Target rules may not contain wildcards. Please specify concrete files or a rule without wildcards at the command line, or have a rule without wildcards at the very top of your workflow (e.g. the typical "rule all" which just collects all results you want to generate in the end)."

It looks like I need to do this "rule all" step, any idea how to set this up?
This is how the snakemake file looks like:
SAMPLES = ["A", "B", "C"]

rule bwa_map:
    input:
        "data/genome.fa",
        "data/samples/{sample}.fastq"
    output:
        "mapped_reads/{sample}.bam"
    shell:
        "bwa mem {input} | samtools view -Sb - > {output}"

rule samtools_sort:
    input:
        "mapped_reads/{sample}.bam"
    output:
        "sorted_reads/{sample}.bam"
    shell:
        "samtools sort -T sorted_reads/{wildcards.sample} "
        "-O bam {input} > {output}"

rule samtools_index:
    input:
        "sorted_reads/{sample}.bam"
    output:
        "sorted_reads/{sample}.bam.bai"
    shell:
        "samtools index {input}"

rule bcftools_call:
    input:
        fa="data/genome.fa",
        bam=expand("sorted_reads/{sample}.bam", sample=SAMPLES),
        bai=expand("sorted_reads/{sample}.bam.bai", sample=SAMPLES)
    output:
        "calls/all.vcf"
    shell:
        "bcftools mpileup -f {input.fa} {input.bam} | "
        "bcftools call -mv - > {output}"

Read the snakemake documentation. all is a pseudo-rule that accepts a set of one or more final output files as input and does nothing with them.

One clarification: Snakemake defaults to trying the first rule in the file (whether called the idiomatic all or anything else) and since in your case that has a sample wildcard, it doesn't know what to fill in for that wildcard.  (A typical all rule might use something like input: expand("mapped_reads/{sample}.bam", sample=SAMPLES) for example.)

Even more to the point: move bcftools_call to the top of the file and it will behave as your all rule.


text:
How to get coding sequences (CDS) of UniProt reference proteomes?
I downloaded the following reference proteomes:
https://ftp.ebi.ac.uk/pub/databases/reference_proteomes/QfO/
'Drosophila':'UP000000803_7227.fasta',
'Arabidopsis':'UP000006548_3702.fasta',
'C.elegans':'UP000001940_6239.fasta',
'Zebrafish':'UP000000437_7955.fasta',
'human':'UP000005640_9606.fasta',
'mouse':'UP000000589_10090.fasta',
'yeast':'UP000002311_559292.fasta'

Now, I would like to have the corresponding coding sequences (CDS). My temporary solution is to download CDS files from ENSEMBL FTP server, translate those and match with the ref proteome. However, this only covers 50-90% of protein sequences.
Thanks!

I went through diverse README files and found the answer:
Coding DNA Sequence FASTA files (*_DNA.fasta)
=============================================

These files contain the coding DNA sequences (CDS) for the protein sequences
where it was possible.
The format is as in the following example (UP000005640_9606_DNA.fasta):

>sp|A0A183|ENSP00000411070
ATGTCACAGCAGAAGCAGCAATCTTGGAAGCCTCCAAATGTTCCCAAATGCTCCCCTCCC
CAAAGATCAAACCCCTGCCTAGCTCCCTACTCGACTCCTTGTGGTGCTCCCCATTCAGAA
GGTTGTCATTCCAGTTCCCAAAGGCCTGAGGTTCAGAAGCCTAGGAGGGCTCGTCAAAAG
CTGCGCTGCCTAAGTAGGGGCACAACCTACCACTGCAAAGAGGAAGAGTGTGAAGGCGAC
TGA

So, EBI itself provides the corresponding coding sequences, which they name _DNA.fasta. I misinterpreted the file-name as genomic (unspliced) DNA.


text:
Where to find database for breast cancer phenotype?
Where can I find databases for breast cancer mutations (SNPs & indels) phenotype, with chromosome, position, ref, alt, gene, clinical significance.
So far, I have only found the Clinvar database.
I need to extract data for breast cancer, in order to train a machine learning model that will be able to predict the pathogenicity of variants (benign, pathogenic etc).
So I did the extraction of data from clinvar vcf file (chromosome, ref, alt, clinical significance), but the clinvar data for breast cancer is very small.

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.

Do you mean this (a product of the BRCA Challenge project):

https://brcaexchange.org/

yes it helps thank you

in addition,  Maryem,  • the following may need to note • To gain biological insights, we frequently associate mutational signatures with factors such as driver mutations, germ line variation, epigenetic modifications, and environmental or therapeutic exposures. However, these remain associations until causation is proven.

[link] (https://www.researchgate.net/publication/353501050_Mutational_signatures_emerging_concepts_caveats_and_clinical_applications)

https://drive.google.com/file/d/18JPIjuPZPM2hnj5O18T4HeQnNil3fU9u/view?usp=sharing

https://drive.google.com/file/d/1ARmO9q9JooamZ_WgNMRuwIjvxKqtMaj5/view?usp=sharing

https://ugc.futurelearn.com/uploads/files/ca/71/ca712acd-e23a-4b3b-b228-45ee17ba1d9e/1_10_HGVS_nomenclature_KTB_AF_CL_KTB_AF2_JD_AF3_KTB_FINAL.pdf


text:
how to calculate the expression values per gene for all cells
I'm currently working with a seurat object and I'd like to calculate the expression values per gene for all cells within a particular cluster.
I've gone ahead and subsetted the cluster of interest.
tumor <- subset(x = data, idents = c("Tumor"), invert = FALSE)

How would I calculate the expression values per gene per cell?
FindAllMarkers

Finds markers (differentially expressed genes) for each of the identity classes in a dataset per ident.

No Seurat function is needed. In the seurat object "cells", you can access the raw counts or normalized counts using cells@assays$RNA@counts (raw counts) or cells@assays$SCT@data (SCT normalized data)


text:
Multiple sequence alignment with gaps
I'm working on a problem where the elements of multiple arrays need to be aligned with each other. While doing some research on algorithms that could possibly be applicable to this task, I stumbled on sequence alignment. I believe this is very similar to what I'm trying to achieve.
In short, I need to have equal elements aligned with each other across multiple arrays. However, the difference with my problem compared to sequence alignment seems to be that the latter doesn't tolerate gaps. On the other hand, gaps would be perfectly acceptable for me. So, for example, given input:
[A, ,B,A, ]
[ ,C, ,D,A]
[B, ,C,B,A]
A desired output would be:
[A,B, , , ,A]
[ , ,C, ,D,A]
[ ,B,C,B, ,A]
But there cannot be any errors. That is, every column must contain elements of the same value.
Would it be possible to achieve the desired result with sequence alignment algorithms?


text:
Is it possible to cluster a sc-RNAseq expression matrix with weighting or embedding?
I have an RNAseq mouse retina gene expression matrix(24659,49300) to cluster effectively. Here is a sample(7,4):
gene            cell1   cell2   cell3
KITL            0       0       1
TMTC3           3       0       0
CEP290          1       3       0
4930430F08RIK   2       1       2
1700017N19RIK   0       0       0
MGAT4C          0       0       4

Is there any chance to apply a weighting schema or embedding to increase ARI, NMI or DBI? Dataset is Macosko, 2015: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE63472
I want to try a new approach to cluster rnase data. Question is effectiveness to use weighting schema or embedding. Is it worth to try? May be someone tried it before.

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.


text:
How can I get data of single cell RNA sequence with raw count?
I am using dataset GSE85241, but I can't find the read counts of the dataset. It only provides with RPKM values.
How can I find the read counts of a dataset?

Hi @KevisLin I think the answer below is the correct answer, haci is an expert. Please consider accepting the answer

If the authors did not provide the raw counts, and I don't think it is compulsory, unfortunately you will need to download the raw data, align and count.


text:
Filter bam with sambamba
I have a very large bam file and I want to filter it to keep only a handfull of positions as defined in a bed file. Can I do this with sambamba or do i need another tool?
Thanks

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 think it is possible:
sambamba view --regions regions.bed --with-header --format=bam reads.bam  > filtered.bam

Hi @Pablo could you kindly edit your question inserting this as a suggestion?

You need bedtools, specifically bedtools intersect command, not sambamba:
bedtools intersect -abam A.bam -b B.bed > A_in_N.bam


text:
Compare my VCF to gnomAD variants
I have a VCF with small variant calls against HG38 and I would like to determine which of those calls are present in the gnomAD database.
Is there an existing tool that can do this? Should I be looking at variant annotation tools? or is this something that I would look at BCFtools for?

You can use vcfanno to annotate your vcf with the gnomad vcf file. There's an example config file for annotating with gnomad (along with other databases such as dbSNP) here.

I haven't checked how fast this is but you could try using Hail (disclaimer: I'm a maintainer). You'll need to install the GCS connector first to read from Google Cloud Storage:
curl https://broad.io/install-gcs-connector | python3

import hail as hl

mt = hl.import_vcf('your.vcf')
mt = hl.split_multi_hts(mt)  # use the same biallelic representation that gnomAD uses
mt.write('data.mt')  # store in Hail's native format for faster queries
mt = hl.read_matrix_table('data.mt')
sites = hl.read_table(
    'gs://gcp-public-data--gnomad/release/3.1.2/ht/genomes/gnomad.genomes.v3.1.2.sites.ht'
)
mt = mt.anti_join_rows(sites)

That'll filter down to just those variants (defined by chromosome, position, and alleles) which are not present in gnomAD v3.

If you have sufficient available disk space, you could download the files from the gnomad downloads page and use a bedtools intersect command to find the overlapping variants.
However, the recommended platforms for gnomAD are to use hail utilities developed by them and others, due to the constraints of working with ~1TB of data.
I am not an expert but potentially someone else knows more about using the site hail data.


text:
Problem with setting fixed topology for MrBayes
I need to use partitioned DNA input (concatenated multiple genes) and fixed topology (species tree from ASTRAL) to get branch lengths in MrBayes. I have input NEXUS like:
#NEXUS
BEGIN TAXA;
  DIMENSIONS NTAX = 318;
  TAXLABELS
    Achillea_nobilis_SAMN11585351
    Anthemis_arvensis_SAMN11585352
    ...
    Haplopappus_sp_SAMN11585375
  ;
END
BEGIN TREES;
  TRANSLATE
    1       Achillea_nobilis_SAMN11585351,
    2       Anthemis_arvensis_SAMN11585352,
    ...
    318     Haplopappus_sp_SAMN11585375
  ;
  TREE treeshrink_exons_good_nobrl = [&R] (1,((2,3...);
END;
BEGIN DATA;
dimensions ntax=318 nchar=169024;
format datatype=dna missing=? gap=-;
MATRIX
Achillea_nobilis_SAMN11585351               TTATGCAA...
Anthemis_arvensis_SAMN11585352              TTATGCAA...
...
Haplopappus_sp_SAMN11585375                 TTATGCAA...
;
END;
BEGIN MRBAYES;
  charset At1g09700 = 1-213;
  charset At1g10740 = 214-456;
  ...
  charset At5g66280 = 168530-168757;
  charset At5g66680 = 168758-169024;
  outgroup Achillea_nobilis_SAMN11585351;
  partition by_gene=467:At1g09700,At1g10740,...,At5g66280,At5g66680;
  set partition=by_gene;
  prset topologypr=fixed(treeshrink_exons_good_nobrl);
  prset applyto=(all) ratepr=variable;
  lset applyto=(all) nst=mixed rates=gamma ngammacat=10;
  unlink tratio=(all) revmat=(all) omega=(all) statefreq=(all) shape=(all) pinvar=(all) correlation=(all) ratemultiplier=(all);
  mcmcp ngen=10000000 nruns=4 nchains=4 samplefreq=5000;
  ss;
  sump relburnin=yes burninfrac=0.25 nruns=4;
  sumt relburnin=yes burninfrac=0.25 nruns=4;
END;

I'm not 100% sure about the MRBAYES block, but generally I think it's correct. When I launch it (MrBayes 3.2.7a) I get an error which I'm unable to resolve:
$ mb sp_treeshrink_exons_good.mb.nex
...
   Executing file "sp_treeshrink_exons_good.mb.nex"
   UNIX line termination
   Longest line length = 169069
   Parsing file
   Expecting NEXUS formatted file
   Reading taxa block
      Allocated taxon set
      Defining new set of 318 taxa
   Exiting taxa block
   Reading trees block
      Successfully read tree 'treeshrink_exons_good_nobrl' # NOTE
   Exiting trees block
   Deleting user trees
   Deleting previously defined taxa
   Reading data block
      Allocated taxon set
      Allocated matrix
      Defining new matrix with 318 taxa and 169024 characters
      Data is Dna
      Missing data coded as ?
      Gaps coded as -
      Taxon   1 -> Achillea_nobilis_SAMN11585351
      Taxon   2 -> Anthemis_arvensis_SAMN11585352
...
      Taxon 318 -> Haplopappus_sp_SAMN11585375
      Successfully read matrix
      Setting default partition (does not divide up characters)
      Setting model defaults
      Seed (for generating default start values) = 1654610114
      Setting output file names to "sp_treeshrink_exons_good.mb.nex.run<i>.<p|t>"
   Exiting data block
   Reading mrbayes block
      Defining charset called 'At1g09700'
      Defining charset called 'At1g10740'
...
      Defining charset called 'At5g66280'
      Defining charset called 'At5g66680'
      Defining partition called 'by_gene'
      Setting by_gene as the partition, dividing characters into 467 parts.
      Setting model defaults
      Seed (for generating default start values) = 1901076251
      Setting Topologypr to Fixed for partition 1
      Setting Topologypr to Fixed for partition 2
...
      Setting Topologypr to Fixed for partition 466
      Setting Topologypr to Fixed for partition 467
      Could not find tree 'treeshrink_exons_good_nobrl' # NOTE
      Could not set fixed topology from user tree 'treeshrink_exons_good_nobrl' # NOTE
      Error when setting parameter "Topologypr" (2)
      The error occurred when reading char. 26-52 on line 1445
         in the file 'sp_treeshrink_exons_good.mb.nex'
...

So it correctly reads tree treeshrink_exons_good_nobrl and after while it complains it couldn't find tree treeshrink_exons_good_nobrl...? I can't find where this weird error comes from.

The answer is pretty simple. What you are requesting is fixing the tree prior at a tree you are specifying. The getter/setter is complaining it cannot fit the data to the specified tree because of 26 bp at 26-52 on line 1445 of the nexus file.
The first thing needed is to post output to:
sed -n '1445p' ./sp_treeshrink_exons_good.mb.nex

What I suspect you will find are non-nucleotide characters associated with "missing data" for example xxxxx, or 0000 what the nexus file wants is ????? for missing data. It is possible part of the degenerate nucleotide code is not being observed (long while since I used MB), e.g. RYB I think its unlikely.
The second approach, if the output of the sed command is pretty normal is to remove the prior setting the tree, i.e. the tree you are supplying is violating e.g. Jukes-Cantor correction.
Finally and a related issue is there's an alignment error, i.e. the alignment at this region is out of frame for these residues and Jukes-Cantor is violated so it refuses to proceed. This would be really unusual however, because I've used MB where there are clear pockets of saturation and it hasn't complained.
Personally, I'm impressed by the level of foresight in MB under written into this getter/setter (its checking your input). This algorithm is rarely used but is clearly well written.

Okay there is no line 1445. What I recommend is placing the alignment in either Sequtron or Clustal X and scrolling down at the position 26-52 and if not successful 'optical analysis' of the alignment.
Personally, I would use Biopython AlignIO for this task and Counter from the Collections library to interrogate the bug. However, thats coding stuff and optical analysis does the same thing. Essentially, you loop through each sequence and Counter will make a frequency count of everything in each sequence. You then say if you are not AGCT- tell me what taxa this is. If you code in Python thats a good follow-up question (separate question).

Thank You, @M__ for Your answer. The output is bit unexpected:
$ sed -n '1445p' ./sp_treeshrink_exons_good.mb.nex
  prset topologypr=fixed(treeshrink_exons_good_nobrl);

The data matrix ends at line 973 (sequences are on slngle lines). I don't find any weird characters in data matrix.

Thanks @Tilia I've updated the question. Its as far as I can go on this question I'm afraid.

OK, thank You for Your help, @M__

No worries. Please remember to  upvote (if it helps) or else a "resolved" if its resolved.

I don't think it's problem with alignment. It complains about line 1445 of input, which contains command to fix the topology, and character range fits the name treeshrink_exons_good_nobrl, so I think the issue is related to it.


text:
beginner RNA-seq Replicate papers
have a good R and statistical analysis background (also with machine learning). in addition, I'm a fresh biotechnology grad. I would like to try to replicate some Rna-seq analysis with R papers (with their provided data). Any SHORT (beginner-friendly) papers to recommend?

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

I like Gierlinksi, et al., Bioinformatics, 2015 for this type of problem.
I like this paper because the methods are descriptive and clear, they have lots of figures that you can reproduce along the way, they spend a lot of time discussing "normalization" and quality control, and the data is easy to download.
Even better, the authors spent a lot of time making sure all the samples were split across lanes and flowcells in a consistent way so that both biological and technical differences could be robustly identified.
They also compared a few different methods published at the time for performing differential expression analysis, so you can see how you can process the same data with different statistical models.
Unfortunately, they don't provide the explicit code they run for this paper, so you can't look at their code as a reference to guide you.
But the methods are descriptive enough that I don't think there should be too much difficulty.

I am a beginner too.  I will try the galaxy (training):

https://training.galaxyproject.org/training-material/topics/transcriptomics/tutorials/ref-based/tutorial.html

Here's https://drive.google.com/file/d/1Lwh6ofxi_bb92mB0HoRrMJObE1MmtcMi/view?usp=sharing
It is from the book of Methods of mathematical oncology
https://github.com/ncbi/sra-tools

Suzuki, T., Poignard, C., Chaplain, M., & Quaranta, V. (2021). Methods of mathematical oncology : Fusion of mathematics and biology, osaka, japan, october 26-28, 2020 (1st ed. 2021 ed.). Singapore: Springer Nature Singapore : Imprint: Springer.

Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.


text:
Should I expect / use a prokaryote-specific substitution model when building a phylogenetic tree?
I have a couple of prokaryotic genes that I'm trying to build phylogenetic trees from.
I inputted them into IQTree, which, as part of the workflow, runs through multiple substitution models until it finds the one that seems to fit the data best.
My issue: quite a few of the genes ended up with eukaryotic_specific substitution models (such as Q.plant or Q.insect). Is this likely to bias my data? Should I re-run the analysis and remove the substitution models that are meant for different groups?
Thank you for your time!

That is a complicated question @Laura and is certainly worrying. I do know insect mutation matrices and compositional bias but cannot think why that would be confused for bacteria. What is the %GC of the bacteria of interest? Secondly, what is the data, I presume it is nucleotide data.

Hello @M__ ! It's protein data, I forgot to mention. The GC is roughly 50-55%, with a few outliers (e.g. my outgroup has ~60% GC content).

Thanks only the protein bit matters ... and yes it is likely to make a different. Could you tell me whether it is using any other protein mutation matrix e.g. DayHoff, JTT or is Q.insect the matrix? It is fairly easy to generate a mutation matrix, rather than use an established one.

Thank you @M__    . The program is running the data on a lot of amino-acid exchange rate matrices (including JTT, Dayhoff, and a lot of variants of each accounting for rate heterogeneity across sites http://www.iqtree.org/doc/Substitution-Models). It seems to pick the best model on something called Bayesian Information Criterion? (it also lists the best model according to Akaike Information Criterion and Corrected Akaike Information Criterion but it doesn't look like they influence later steps as far as I can see). Overall (with the variants) it goes through about ~700 models each time.

... is there any chance the bacteria are insect symbionts @Laura?

IQTree is a fashionable method for speeding up bootstrapping and used on ultra large datasets. What IQTree might be doing is rigidly fixing a model and parameterisation is not estimate via the ML. This might be how it achieves its performance.
In context, IQTree can use an establised mixed model, which is one approach to protein phylogeny. As a personal note, I don't favour it and come back to the reason later. It MIGHT also generate its own mixed model. There are not 700 established (published) amino acid matrices plus established mixed models, combined there are no more than 50.
There are around 30 protein mutation matrices, Dayhoff models (I honestly don't who Dayhoff is), one for insects which might be for mtDNA ... mtArt, one for mtDNA, mtREV I think and so on. Please check carefully the number of models it is evaluating.
Its easy to see how 4 or 5 combinations (or possibly permutations) of 30 become 700 models if it is generating its own mixed model.
Alternatively, what it does (mixed model or not) BIC and AIC do is given a tree estimate (prior) what's model is the best probabilistic fit. In your case Q.insect is cropping up, whether this is mixed model or not ... dunno.
What I don't know is how it calculates the mixed model, is it creating a new singular matrix by averaging (which is probably closer to the algorithm), or is it sequentially passing the data through first one model, then another etc ... Personally, I avoid mixed models for precisely this reason, I dunno what the actual matrix is, but accept it can provide strongest AIC/BIC.
It suggests either:

BIC/AIC results in an outcome that isn't biologically realistic .... OR
there has been some weird parallel/convergent evolution - this happens in parasites for certain

If you are using a single matrix not mixed model and Q.insect is cropping up - I would be concerned. If it is part of a mixed model, I wouldn't worry.
Is it checking 30 models (plus published mixed models) or 700 matrices??
If the final point was correct (biologically real result), you check it via blast - simply exclude bacteria and see if you get a plant of insect sequence as the top hit.
Ultimately, the uncertainty is the reason I avoid mixed models.

To answer the question, I'm not aware of a bacterial amino acid matrix. HOWEVER, you can select a mixed model. It might just be screening around 50 models via AIC/BIC and within that Q.insect is the most probable.
Summary
The alternative is to look for support using a general amino acid matrix (for example JTT, LG or JTTDCMut).
Personally, I would be very cautious about proceeding with a matrix generated by insect genomes outside a mixed model for bacteria and in the absence of any biological support. FOR EXAMPLE, if these were insect symbionts - that makes sense.
Which ever approach is used a Bayesian approach is often required.


text:
How can I download from NCBI all the ITS genes and the related taxonomy?
I would like to download all the ITS1 and ITS2 genes from NCBI in a fasta file. And, I'd like to download even the related taxonomy of each sequence.
Thanks,
Marco

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 recommend to use Biopython:
from Bio import Entrez

Entrez.email = "your@somemail.com"
Entrez.api_key = "YOUR_API_KEY"

search_result = Entrez.read(Entrez.esearch(db="nucleotide", term="ITS1[Gene Name]"))
with Entrez.efetch(
    db="nucleotide", id=search_result["IdList"], rettype="fasta", retmode="text"
) as handle, open(filepath, "w") as file:
    file.write(handle.read())

For batch operations, it would be better to acquire a NCBI API key and combine GNU parallel and Entrez Direct .
If you already have the accession numbers of ITS1 and ITS2 genes, put the accession numbers in a text file like "accession-number.txt". Otherwise search the ITS genes you want to download and then get the accession number:
esearch -db nuccore -query 'you query term[Gene Name]' |
  efetch -format acc >accession-number.txt

To download all the sequences in fasta format.:
parallel -a accession-number.txt -j8 'efetch -db nuccore -id {} -format fasta >{}.fa'

Here -j8 is used to control the number of jobs run in parallel, it should be less than the limit of 10.
To get related taxonomy (for eg, taxonomy rank, scientific name, taxonomy id, genus and species) of each sequence:
body='BEGIN { OFS = "\t" } { print "{}", $0}'
paralle -a accession-number.txt -j8 "elink -db nuccore -id {} -target taxonomy -name nuccore_taxonomy |
  esummary |
  xtract -pattern DocumentSummary -element Rank ScientificName TaxId Genus Species |
  awk '$body'"

Thank you! I thought, that they do not allow to make more than 10 requests per second for an API key regardless the fact it is parallel or serial. Which is real restrictions?

I think the real restriction is 10 requests per second. You can try -j with a number larger than 10 to see if there are any errors.

To get taxonomy info without awk:  paralle -a accession-number.txt -j8 "elink -db nuccore -id {} -target taxonomy -name nuccore_taxonomy |   esummary |   xtract -pattern DocumentSummary -pfx '{}\t' -element Rank ScientificName TaxId Genus Species"


text:
Identification of unique chains in proteins by sequence similarity
I'm trying to find unique chains or proteins within their PDB files. Many proteins have multiple chains, but very often they are identical (say, the PDB file consists of 4 homodimers, for example, so there is one unique chain represented 8 times total). I am processing certain structural features of the chains, so processing two identical chains is redundant for my pruposes.
I am using Biopython and the way I'm currently identifying unique chains is looking through each PDB file, checking the sequences of the chains and comparing them to one another using alignments (specifically pairwise2.align.globalxx(seq1, seq2). I consider a chain non-unique if it shares >90% sequence similarity to another chain. (Most chains have 100+ residues.)
My question is, is this threshold appropriate? Should I go lower, maybe 80%, or the opposite, higher? For now, it seems to work fine, but I'll be processing a lot of PDB files and want to be on the safer side. If I accidentally process identical chains, that's not much of an issue, but it adds unnecessary time to the processing. Is there maybe a source/article to help with what sequence similarity I should be using? I've found multiple articles, but they usually talk about homology (i.e., common ancestry of proteins), but I'm only doing this within the protein itself and its subunits (chains), so this is of no use for me.

I commented yesterday w/ a link to SIFTS —sorry about the terseness: I meant to say why not use SIFTS table of PDB chain to Uniprot if you are just after finding which chains are the same. Alternatively, the PDB API data will tell you in a more convoluted way what is present. Neither will tell you about gaps in sequence and the latter will be weird for C-Alpha traces (as the precise AA is unknown)

Came across this looking for a solution myself last night but now I have created my own I will share it.
Rather then using Biopython (which is a Python package I use a lot), it is more computationally efficient to parse the .fasta files from the protein data bank (RCSB in this case) which as far as I can see, contain this information. For example, this is 2P16's .fasta file (with 2 unique chains):
>2P16_1|Chain A|Coagulation factor X (EC 3.4.21.6) (Stuart factor) (Stuart-Prower factor)|Homo sapiens (9606)
IVGGQECKDGECPWQALLINEENEGFCGGTILSEFYILTAAHCLYQAKRFKVRVGDRNTEQEEGGEAVHEVEVVIKHNRFTKETYDFDIAVLRLKTPITFRMNVAPACLPERDWAESTLMTQKTGIVSGFGRTHEKGRQSTRLKMLEVPYVDRNSCKLSSSFIITQNMFCAGYDTKQEDACQGDSGGPHVTRFKDTYFVTGIVSWGEGCARKGKYGIYTKVTAFLKWIDRSMKT
>2P16_2|Chain B[auth L]|Coagulation factor X (EC 3.4.21.6) (Stuart factor) (Stuart-Prower factor)|Homo sapiens (9606)
KLCSLDNGDCDQFCHEEQNSVVCSCARGYTLADNGKACIPTGPYPCGKQTLE

and 6WON's (with 8 of the same chains):
>6WON_1|Chains A, B, C, D, E, F, G, H|YohF|Salmonella typhimurium (strain LT2 / SGSC1412 / ATCC 700720) (99287)
QGMTKVAIVTASDSGIGKACALLLAQNGFDIGITWHSDERGAQETAKKAAQFGVRAETIHLDLSQLPEGAQAIEHLIQRLGRVDVLVNNAGAMTKSAFIDMPFTQWRQIFTVDVDGAFLCAQIAARHMIKQGEGGRIINITSVHEHTPLPQASAYTAAKHALGGLTKSMALELIEHHILVNAVAPGAIATPMNDMDDSDIEPGSEPSIPIARPGSTHEIASLVAWLCSEGASYTTGQSLIVDGGFMLANPQFNAK

So far this always seems to be the case and parsing the .fasta is a lot easier then using BLAST or sequence similarity such as you are or were.
For the implementation, whereas the .pdb, .cif files etc use the same download URL (https://files.rcsb.org/download/), the .fasta files use this link: https://www.rcsb.org/fasta/entry/. Parsing the downloaded .fasta files and only reading the lines that begin with >, stripping out unnecessary chains (in my specific case I just want the first unique chains in the PDB so A, B being the same and C, D being the same would return me a list with [A, C], you can modify this script to your own ends) and running the Python code gives this output:
First unique chains in 1d2z (2 unique chains): ['B', 'A']
First unique chains in 2p16 (2 unique chains): ['A', 'B']
First unique chains in 6won (1 unique chain):  ['A']

There was an issue with some entries having Chain or Chains so both had to be considered. Additionally, there are commas that can be removed from all but the last chain ID but for my purposes I do not need to strip it, you can put something in _fasta_helper() to this end before return. This is the Python code and hopefully it helps someone out in the future:
import re
from typing import Union
import urllib.request
from urllib.error import (HTTPError,
                          URLError,
                          ContentTooShortError as ContentError)

def rcsb_downloader(pdb: str, ext: str = None, down_path: str = None) -> str:
    """Downloads the PDB from the protein databank RCSB.

    Parameters
    ----------
    pdb : str
        Protein databank accession 4 letter code i.e. '6won'
    ext : str, default = None
        Extension of the file to download i.e. '.fasta', '.cif' etc. Note: if
        None argument then ext = '.pdb'

    Returns
    -------
    str | bool
        Returns str path if download successful, False if failure
    """

    if ext is None:
        ext = '.pdb'
    if down_path is None:
        download_path = ('/your/path/to/dir/' + pdb + ext)
    else:
        download_path = down_path + pdb + ext

    rcsb_path = 'https://files.rcsb.org/download/' + pdb + ext

    if ext == '.fasta': # fasta has a different web search then others
        rcsb_path = 'https://www.rcsb.org/fasta/entry/' + pdb

    try:
        urllib.request.urlretrieve(rcsb_path, download_path)
        return download_path
    # Try again in case there was a temporary issue with the first attempt
    except (ContentError, HTTPError, URLError):
        try:
            urllib.request.urlretrieve(rcsb_path, download_path)
            return download_path
        except (ContentError, HTTPError, URLError):
            print(f"PDB failed to download {pdb}")
            return False

    def _fasta_helper(self, line: str, single_unique_chain: bool = True) -> Union[str, list]:
        """Separate the first chain ID from the .fasta comment line."""

        # [1] some of the chains had [auth C] etc after so this is removed
        # [2] was originally a list comprehension that filtered based on the word
        # chain, however some of the fasta entries had "Chain" or "Chains" in their
        # description which broke this. The first arrow is removed, split at |, the
        # second entry (1st in Python) entry contains the chain information and the
        # first entry removed (which so far has always been chain, anticipate their
        # being edge cases with older PDB entries). This appears to hold up.

        clean_chain = [
            x.replace(",", "")
            for x in list(
                re.sub(r"\[.*?\]", "", line).lstrip(">").split("|")[1].split(" ")[1:]  # [1]  # [2]
            )
        ]

        if not single_unique_chain:
            return clean_chain
        return clean_chain[0]

def fasta_parse(path_: str) -> list:
    """Parse .fasta files and analyse lines beginning with '>'."""

    with open(path_, 'r', encoding='utf-8') as fasta:
        total_unique_chains = 0  # Use as a flag to test if more than one unique chain
        unique_chains       = []
        for line in fasta:
            if line[0] == '>' and total_unique_chains == 0:
                total_unique_chains += 1
                unique_chains.append(_fasta_helper(line))
                continue
            elif line[0] == '>' and total_unique_chains != 0:
                unique_chains.append(_fasta_helper(line))
    return unique_chains

pdb_diff_chain  = '2p16' # Chain ...
pdb_diff_chains = '1d2z' # Chains ...
pdb_same_chains = '6won' # Chains ...

download_p_diff  = rcsb_downloader(pdb_diff_chain, '.fasta')
download_p_diffs = rcsb_downloader(pdb_diff_chains, '.fasta')
download_p_same  = rcsb_downloader(pdb_same_chains, '.fasta')

unique_s  = fasta_parse(download_p_same)
unique_d  = fasta_parse(download_p_diffs)
unique_ds = fasta_parse(download_p_diff)

print(f'First unique chains in {pdb_diff_chains} (2 unique chains): {unique_ds}\n'
      f'First unique chains in {pdb_diff_chain} (2 unique chains): {unique_d}\n'
      f'First unique chains in {pdb_same_chains} (1 unique chain):  {unique_s}')

This is a great way, but unfortunately I have the files as PDBs, so I'd have to re-download everything as FASTA files instead. However, I really like your solution, thank you for it! It's a great way to do it and definitely better than processing sequences from PDBs (due to chain breakings and its processing by BioPython.Polypeptide module), so I'll keep this in mind when I need to do it again.


text:
How to make a network overlap of gene lists across 3 conditions like this image?
From 3 gene lists (lung, nasal, gut), how to make this gene overlap network plot?
This image is figure 4a from: https://www.nature.com/articles/s41591-020-01227-z

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'm the author of gravis, an open source package for interactive graph visualization in Python.

It can visualize directed and undirected graphs (=having edges with or without arrows) as well as multigraphs (=having more than one edge between a pair of nodes).
It allows to control visual properties such as node sizes, colors, shapes or edge widths and color with data.
It recognizes graph objects from various network analysis packages such as NetworkX, igraph, graph-tool and others.
The output is based on web technologies (HTML/JS/CSS), so it can be viewed in a browser or Jupyter notebook. There is also support for exporting static images (JPG, PNG, SVG).
It comes with many examples on the documentation site.

Here's a quick attempt to reproduce a part of the figure you've provided:
import gravis as gv
import networkx as nx

s1 = 40
c1 = 'white'

s2 = 20
c2 = '#ffd703'

s3 = 20
c3 = 'orange'

s4 = 20
c4 = 'red'

g = nx.Graph()
g.add_node('Lung', color=c1, size=s1, border_size=1)
g.add_node('Nasal', color=c1, size=s1, border_size=1)
g.add_node('Gut', color=c1, size=s1, border_size=1)

g.add_node('MDFI', color=c2, size=s2)
g.add_node('LOXL1', color=c2, size=s2)
g.add_node('CNBP', color=c2, size=s2)

g.add_node('C3', color=c3, size=s3)
g.add_node('CA2', color=c3, size=s3)
g.add_node('SEPP1', color=c3, size=s3)

g.add_node('TFF3', color=c4, size=s4)
g.add_node('PTMA', color=c4, size=s4)

g.add_edge('Lung', 'MDFI')
g.add_edge('Nasal', 'LOXL1')
g.add_edge('Gut', 'CNBP')

g.add_edge('Lung', 'C3')
g.add_edge('Nasal', 'C3')

g.add_edge('Lung', 'CA2')
g.add_edge('Gut', 'CA2')

g.add_edge('Nasal', 'SEPP1')
g.add_edge('Gut', 'SEPP1')

g.add_edge('Lung', 'TFF3')
g.add_edge('Nasal', 'TFF3')
g.add_edge('Gut', 'TFF3')

g.add_edge('Lung', 'PTMA')
g.add_edge('Nasal', 'PTMA')
g.add_edge('Gut', 'PTMA')

gv.d3(g)

Resulting output in a web browser:

I've moved the nodes manually a bit so that the layout resembles the provided image. Once everything looks as intended, a static image can be exported with a single button click on the right-hand side.
There are many other ways to produce static graph visualizations, e.g. programmatically with other Python libraries mentioned before, or interactively in standalone tools such as Gephi, Cytoscape or Tulip.

You could use OmniGraffle to make this figure.
Objects in OmniGraffle have so-called magnets to which you can attach lines, which can join two objects.
Once you have objects joined, you could group them and position them as desired.
As you drag an object, any attached lines will also change their position as needed. You can drag groupings of objects (e.g., conditions or other categories of objects), which does the same thing.
You can export an illustration to PDF or SVG format from OmniGraffle, which will let you do any further work in Adobe Illustrator or other vector illustration tools, or share the illustration with colleagues over the web, for instance.

This question is a bit difficult to fathom, but essentially I think the answer you are look for is called network theory. Its quite simple, but alas not fashionable. You'll find an introduction here and try and avoid 'directed graphs' because that gets complicated.
The alternative way is to use Google fusion tables, which incorporates network analysis, which is really easy to use and investigate. Its very interactive. You'll find more information here,
https://support.google.com/fusiontables/answer/2566732?hl=en
Thats it and nothing more, undirected network theory is pretty simple. Its does not tell you very much (in my opinion), but it makes pretty diagrams.
Please note there is an analysis called 'Graph theory' and is not needed here. Whilst this has some similarity to 'network theory' it is distinct and the two shouldn't be confused.


text:
Issue with iCTNet2 installation
Has anyone encountered this same installation issue of iCTNet2?
I encountered this issue for installation of iCTNet2 via Cytoscape App Manager and similarly via installation from website download/install page.
I tried both installations via Cytoscape and via website takes a long time, when both is still installing even 4 hours later, I cancelled the installation. Does it really take that long? Or I should just wait for more than 4hrs to install this?
Any recommendations please? Thank you.


text:
Assemble Anaeroplasma species genome from metagenomic PacBio data
I have a fasta file containing reads generated by PacBio HiFi whole genome sequencing of a feces sample from mouse.
I would like to use this dataset to generate an assembled circularized genome for an Anaeroplasma organism known to be in this sample.
I am familiar with Linux and R. I would very much appreciate your advice on the needed steps and the required tools or R packages.

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

Many tools prefer PB reads in their original BAM format. What kind of PB reads? HiFi?

Yes HiFi reads. I have the reads in fastq and fasta fromat.

Try hifiasm-meta (paper here). Whether it works depends on the abundance and the complexity. You may also try HiCanu.

I would suggest the following workflow, understanding you to have metagenomic HiFi sequencing that you know to include your species of interest at high abundance:

metagenomic assembly of HiFi with metaFlyE
Annotate your resulting long contigs (>500Kbp) using GTDB-tk
Look for your Anaeroplasma among long contigs
Assess whether your genome is closed.
If the genome is closed, you are done.
If the genome is not closed, you need to close it. Assess its completeness using e.g. checkM.
Map your HiFi reads to the genome, see if you can close the remaining gap with e.g. SSpace-LongRead.

I'd suggest using this paper as an example for further workflow details.

Thanks for your answer and the references

@Amroon, you are very welcome. By the way, it's traditional on SE sites to upvote answers that you find helpful (there are now a couple of reasonable answers to your question), and accept the answer that you see as helping you solve the problem. For more info on responding to answers on SE, see here: https://stackoverflow.com/help/someone-answers


text:
How to identify most important Amino Acid residues in a given protein?
Continuing from the title, in my context, important means that changing the set of AAs (either or one more many) changes the domain and/or fold and/or function and/or family etc.
I  noticed that there are two types of methods Protbert (protein language models, sequence-based methods) and/or Alphafold (structure (tertiary) based methods); one needs both methods to "reallY" identify the essential set of AAs in a given protein?
Thanks

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

Sorry. This question is ill posed as it completely depends on context —namely is simply a question of what are the main roles AA residues can take. Core residues when mutated may cause and increase in Gibbs free energy of folding (i.e. less energy is released when folded from a magically linear protein form), these are key to Structure as seen in PDB and predicted by Alphafold. Catalytic residues are key to activity. Channel gate residues are key to permeability. etc etc etc etc

Do you have many examples of your protein from different species? Can you look at alignments of the various homologs?


text:
What is the basic difference between a protein and a ligand?
Is ligand a type of protein or vice versa?

A ligand is any molecule that can bind to a receptor. A receptor is always a protein. A ligand can be a protein, but also something else, be it a combination of protein with carbohydrates and-or fats, or just molecules such as cyclic AMP, any chemical compound, structures expressed on the cell surface, pretty much any structure you can think of.

Not only receptor, it may be any protein (as haemoglobine). The term in molecular biology has not a strict definition unlike inorganic chemistry. But it is rare when we look at proteins as ligands (usually little oligopeptides). Commonly it is named as hetero-oligomeric complexes.


text:
Could someone please help me understand the outputs for HapFlk
Could someone please explain the output file for HapFlk?
My input file was a ped file which has three different populations. Two with same number of snps by imputation and third one is an outgroup population which has the same number of snps but some of them are marked as 0 since it's missing and unimputed.
I ran HapFlk with only the --file and --outgroup option so the output would be Flk not the LD based HapFlk and outgroup would be leftout after building the tree
When I checked the flk result after running,
the columns were rs, chr, pos, pzero, flk, pvalue which I could understand but some of the positions have

nan for pvalue
nan for flk and pvalue
1.0 for pzero and nan for flk and pvalue
0.5 for pzero and 0.0 for flk and 1.0 for pvalue

I'm not sure how to interpret these positions cause the positions aren't missing in input ped file.
Please help me

crossposted https://www.biostars.org/p/9526800/

I'm sorry , it's not allowed to crosspost?

It’s not against the rules, but it’s typically seen as polite to comment that you’ve posted elsewhere, to avoid people duplicating answering effort.

Regardless, you are likelier to get a good response to this at Bioinformatics SE :).

I'll keep that in mind thanks


text:
statistical test and visualization of genes present in an experiment
I have a list of genes (n = 120) that are involved in breast cancer. And I have a list of differentially expressed genes (n = 70) that are present in the breast cancer gene list. My question is: Are differentially expressed genes significantly more likely to be breast cancer related vs non-differentially expressed genes?
which statistical test will be suitable to address this question and how I can visualize its results in R.

I would suggest using an existing tool to evaluate this issue, such as the GO enrichment analysis toolkit: http://geneontology.org/docs/go-enrichment-analysis/. Some of the details are finicky, there are a lot of people who have thought very hard about the issue of annotating genes and GO is what people have settled on, rather than a bespoke list of BC genes. If you must visualize this (rather than simply report the numbers or have a table of enriched categories, which I'd prefer), I'd use a Venn diagram, which you can google for R.

This question can be addressed by gene set enrichment analysis. You should use Gene Ontology and KEGG Pathways to find out which genes from your gene set are enriched in GO terms and KEGG pathway terms. You can use clusterProfiler package in R for this anlysis and for visualization also. Here is the link to the guide:
https://yulab-smu.top/biomedical-knowledge-mining-book/enrichment-overview.html


text:
Running hubCheck utility from inside K8s's Pod
I'm trying to run hubCheck utility* from inside Kubernetes Pod but for some reason I get this Network error:
# I ran hubCheck hubCheck from inside the pod
$ ./hubCheck -noTracks -verbose=2 ftp://ftp.ebi.ac.uk/pub/databases/Rfam/12.0/genome_browser_hub/hub.txt
### kent source version 431 ###
Found 2 problems:
Couldn't open ftp://ftp.ebi.ac.uk/pub/databases/Rfam/12.0/genome_browser_hub/hub.txt
TCP non-blocking connect() to ftp.ebi.ac.uk IP 193.62.193.138 timed-out in select() after 10000 milliseconds - Cancelling!

However, when I run hubCheck locally it works fine.
$ ./hubCheck -noTracks -verbose=2 ftp://ftp.ebi.ac.uk/pub/databases/Rfam/12.0/genome_browser_hub/hub.txt
### kent source version 430 ###
1 tracks in sacCer3
1 tracks in mm10
1 tracks in ce10
1 tracks in galGal4
1 tracks in ci2
1 tracks in danRer7
1 tracks in dm6
1 tracks in hg38

The first thing I checked is if the pod has access to the internet, which was the case as you can see below, I can use curl to download hubCheck:
# curl -o tools/hubCheck http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/hubCheck
  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100 33.2M  100 33.2M    0     0  10.0M      0  0:00:03  0:00:03 --:--:-- 10.0M

I thought this may be caused by a firewall limiting the access to external URL, but I don't think it's the case because other URLs like the one below is working fine both outside and inside K8s:
$ ./hubCheck -noTracks -verbose=2 http://younglab.wi.mit.edu/mESC_ChIAPET_hub/hub.txt
### kent source version 431 ###
4 tracks in mm9

Any help on this would be greatly appreciated!

hub Check is a database which facilitates online UCSC Genome Browser display

There are not many Kubernetes users here, but you never know. This question is perhaps better answered via UCSC's online support

@M__, Thank you, I think I'll do that, I asked on SO but the question didn't get enough attention


text:
How can I get latest .vcf files with annotation data?
I'd like to perform the annotation (1000genomes, COSMIC etc.) of my variants using SnpSift and SnpEff, however so far all I get are vcf files for separate chromosomes:
http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/
I'd like to annotate my variants and then filter out common SNP's. What in your opinion is best annotation source and where can I get the file?

Do you mean Somatic Sequence Variants?,
if so, Table 1 inside the following full-text link:
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5707196/
And additional information from my search-->
Here's
https://docs.google.com/document/d/1hVtgMYAtvDoUcigKx8OLDmongDRzDQAd/edit?usp=sharing&ouid=102408774767674945852&rtpof=true&sd=true

https://drive.google.com/file/d/18JPIjuPZPM2hnj5O18T4HeQnNil3fU9u/view?usp=sharing

https://drive.google.com/file/d/1RfIn2pIyr_zCyK18rQCX_J1Vr0OQ9Q4W/view?usp=sharing

https://drive.google.com/file/d/1l6M5xlKulDGQ2aFLz21TuNYY9s2LpFAd/view?usp=sharing

https://drive.google.com/file/d/1s86a0gtDE-15h0r3UORfMux3pT8Ns1ls/view?usp=sharing

According to your answer I tried to get the files which might be useful regarding annotation:

thanks a lot,  how can you do all of that within 24 hrs?  Congratulations!


text:
Can I use samtools addreplacerg to replace multiple RG entries at the same time?
I have a bam file that contains two @RG lines:
@RG     ID:Foo.Bar_L1.Foo   PU:Bar_L1.Foo SM:Foo        PL:Illumina     LB:Library.fa
@RG     ID:Foo.Bar_L2.Foo   PU:Bar_L2.Foo SM:Foo        PL:Illumina     LB:Library.fa

I wish to replace it with one @RG line:
@RG     ID:Foo SM:Foo        PL:Illumina     LB:Library.fa

What is the best way to do this?
samtools addreplacerg looks useful, but it looks like I can only affect one @RG per call, and I have no desire to write out the whole bam file twice

Not sure about samtools, but picard AddOrReplaceReadGroups can clear existing read groups and write a new one.


text:
Heatmap of motifs present in protein sequences
I have a file containing multiple protein amino acid sequences. I am interested in visually presenting the location of amino acid motifs PQG, QQG and AQG in the sequence in the form of a heatmap in R. Every location of the motifs should be represented as a differently coloured barcode. Example sequence follows,
>cat archae.fasta
>A0A075FKY9
MRRITLPKGIEISGKSPFGTSSGVYEVVISISENKPSDASIASNSFNSLWKDNFHLRVSG
DHFSQVLGSGSNPIPDNVFDLGSVWIVIQDQFSPVHTSFQFNISPSNQPASKPEEKPRRV
ETSTVKRTITKPGRLGPPGGKGPRGPSGYPGIKGDKGGRGPTGDKGDKGDKGIPGPQGEK
GKTGQAGDKGDKGITGPPGEKGPRGSTGPPGDKGDKGVQGKIGDKGLTGTTGPVGDKGTQ
GPQGPPGGKGLTGIPGPQGDKGTKGPLGPVGERGPTGKPGESGLQGPQGIQGAQGERGHS
GPPGPQGEKGLIGDKGEIGERGTRGPPGPPGEKGSQGGMSDEGKRLIKELLELLASKNII
STEEQIKLTSYLY
>A0A075FMB0
MPGDKGTSGIRGVPGDKGDKGPQGPPGDKGLTGPVGVPGEKGPTGTSGVQGEKGIQGTPG
PLGEKGVTGPAGDKGPIGPPGPAGSKGITGPAGPLGDKGTIGPPGPLGDKGSKGPEGPTG
DKGPQGPQGPAGAKGLTGVPGPQGEKGEKGPLGPIGEKGTTGATGPPGDKGPQGPQGTQG
ERGTTGPSGEKGPQGPQGIQGPQGERGPTGPIGSIGETGPIGGQGPVGPAGPRGPPGPPG
EKGPSGGMSSEQKAIFKELLEILTTKEIISTEEQIKLMSYLY
>A0A075G2F8
MPKGIEISGKSPFGTSSGVYEVKISVSQSKPSNDSITNNSFNSLWKDNFHLRVSEGYFSE
ILGSDSNPIPDNVFDLGSVWIVIQDQFSPVHTSFEFKISQDSKSIPTPKEKPRRVDTSSV
KKVRTIPTRQGPTGSKGIPGERGYPGTKGEKGGKGPTGDKGDKGDKGVPGPQGEKGKTGS
TGDKGDKGITGSPGEKGPRGPTGPPGDKGDKGVQGKIGDKGLTGTTGPIGDKGPQGPQGP
PGDKGLTGIPGPQGDKGGRGTIGPSGERGPTGRPGESGLQGPQGIQGPQGERGHSGPPGP
QGEKGIIGDKGEIGERGPRGPPGPPGEKGSQGGMSDDGKILFKELLQLLASKNIISTEEQ
IKLTSYLY
>A0A075GLR7
MPKGIEISGKSPFGTSSGVYEVRISISQDKPSDDGITNNSFNSLWKDNFHLRVSEGHFSE
ILGTDSNPIPADVFGLGSIWIIIQDQFSPVHTSFELHISQDNKLLSVPQGKPKPVNATKV
QKVRAQPGRPGPQGEKGGKGPAGYPGTRGEKGGRGITGDKGDKGDKGIPGPQGEKGIKGP
AGDKGDKGITGPVGPSGPKGSTGPPGDKGDKGVQGKSGDKGLTGVTGPLGDKGTQGPQGP
PGPKGLTGVPGPQGDKGMRGSIGPAGERGPTGKPGDPGLPGPQGIQGPQGERGHTGPQGP
PGEKGLVGDKGEIGDRGPRGPPGPPGEKGSQGGMSEESKRLVKELLELLASKNIISTEEQ
IKLTSYLY

As an example, please see a two-colour heatmap. I need something similar but with three colours corresponding to each of the three motifs.

I looked at the heatmaply package but the biggest problem is that my matrix is not numeric. Is there way to make heatmap using text strings in R?
Motif is defined as any three character string that matches the general pattern either PQG, QQG and AQG in the uploaded sequence. y-axis is the fasta code appended after > (A0A075FKY9 for example). x-axis would depend on how the data is converted into a matrix.

Firstly, the problem with your data is that it ain't aligned. Without that your motif positions are not relative from one protein to another. Basically you need to align it and there's no wriggle room.
Secondly, its easy you simply parse the motifs you want and their location within the string - after it is aligned. In Python you would use finditer within a comprehension. This will go into a dataframe and subject to a heatmap package, I'm not a 'heatmap' person but it would be easy to graph where the x axis is the position where each motif is found and a qualitative score given for each type of motif (1-3 here). I don't know the equivalent in R.
There is an issue whether an indel will remove a motif and you would need to check for that by running the same command against an unaligned protein and seeing if the frequency matches.
My advice is to ask a new question with some code. I dunno whether someone has coded this already, - its easy to do. Personally, I don't think there's sufficient requirement, but I could be wrong.


text:
Elegant Way to Extract Gene Symbols from FASTA Headers
I've downloaded the RefSeq transcript sequences (FASTA format) for the GRC38 patch 14 build of the human genome and would now like to extract the gene symbols (only) in order from the headers. I'm confident I can do this using grep and regular expressions--the trick will be coming up with an appropriate pattern. (The gene symbol appears to be surrounded by the last pair of parentheses in a header, but I'm not sure this is absolutely true.) However, I'm wondering if there is a more elegant way to do this--e.g., via the NCBI EDirect utilities or SeqKit--that I'm overlooking. (I have a list of the accession.version numbers.) Any thoughts on this would be greatly appreciated.

It would be helpful to see a sample of the file or a direct link to the file in question. Is there a reason that you want to do this from the FASTA header (which is probably not so strictly curated) rather than from a more traditional annotation file format such as GFF3?

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.

We can't help you parse data if you don't show us what you need to parse. Please [edit] your question and show us a few lines (preferably just fasta header lines) of your file and the output you would want from that example.

We can't be sure without seeing your file, but if what you really want is the text between the last parentheses on lines that start with >, then you can do this with GNU grep (the default on Linux) with:
grep -oP '^>.*\(\K[^)]+' file.fa

The -o means "print only the matching portion of the line" and the -P activates Perl Compatible Regular Expressions (PCREs). The regex looks for a > at the beginning of the line (^>), then as many characters as possible until the last parenthesis (.*\(), and then forgets everything matched up to this point with the \K. This means that the -o will only see things after the \K, a PCRE feature. So finally, we look for the longest stretch of one or more non-) characters.
If you don't have GNU grep, you can do the same basic thing with perl or sed:
sed -n 's/^>.*(\([^)]*\).*/\1/p' file.fa
perl -lne '/^>.*\(([^)]+)/ && print $1' file.fa

Super helpful. Thank you!

You're very welcome, @MarkPauley. If one of the answers here solved your issue, please take a moment and accept it by clicking on the checkmark on the left. That is the best way to express your thanks on the Stack Exchange sites.

There are probably better ways of doing this with Python or another tool, but here's a quick method that works in bash (tested on Linux) assuming that the Gene Symbols are contained in parentheses, i.e. "(Gene Symbol)":
pat='\(([^)]+)\)'
while read line; do
    if [[ "$line" =~ $pat ]]; then
        echo ${BASH_REMATCH[1]};
    fi;
done < human.1.rna.fna

Thanks. Unfortunately, there can be more than one parenthesized expression in the header, so this pattern is too permissive. In every record I've looked at, the gene symbol is in the last parenthesized expression, but I don't know if that is guaranteed to be true.

Could you maybe use another file provided by the RefSeq ftp such as this one? It may have more genes than in your FASTA file, I'm assuming not fewer, but you could always do a quick check using the regex output as a comparison by sorting both files, and performing a diff command in Linux.

You really, really don't want to do this in the shell. It is going to be orders of magnitude slower than anything else and is very hard to get right. For example, your entire loop can be done faster and simpler with grep '\([^)]*\) human.1.rna.fna. The loop only makes it slower and more fragile.

Relevant: Why is using a shell loop to process text considered bad practice?


text:
gtf2bed only retains gene information
I want to convert an hg38 GTF file to bed format using gtf2bed:
gtf2bed Homo_sapiens.GRCh38.106.gtf > Homo_sapiens.GRCh38.106.gtf.bed

It outputs a bed file but it only contains the coordinates for the whole genes and has lost the locations of the exons, UTRs, etc.
While I am running gtf2bed, I get thousands of warning like this:
WARNING:root:None   21  33276625    33276753    +   three_prime_utr
WARNING:root:WARNING End is after gene end (33279918 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  33279752    33279918    +   three_prime_utr
WARNING:root:WARNING End is after gene end (33283212 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  33283094    33283212    +   three_prime_utr
WARNING:root:WARNING End is after gene end (37101343 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  37100814    37101343    +   transcript
WARNING:root:WARNING End is after gene end (37101343 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  37100814    37101343    +   exon
WARNING:root:WARNING End is after gene end (37237744 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  37221419    37237744    +   transcript
WARNING:root:WARNING End is after gene end (37221597 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  37221419    37221597    +   exon
WARNING:root:WARNING End is after gene end (37237744 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  37237404    37237744    +   exon

How can I create a bed file that also contains the coordinates of the other genomic features?

I cannot reproduce this. Be sure to use the latest version of BEDOPS and to use the correct command:

wget http://ftp.ensembl.org/pub/release-106/gtf/homo_sapiens/Homo_sapiens.GRCh38.106.chr.gtf.gz

gtf2bed < <(zcat Homo_sapiens.GRCh38.106.chr.gtf.gz) > Homo_sapiens.GRCh38.106.chr.bed

$ head Homo_sapiens.GRCh38.106.chr.bed
1       1211339 1214153 ENSG00000186827 .       -       ensembl_havana  gene    .       gene_id "ENSG00000186827"; gene_version "11"; gene_name "TNFRSF4"; gene_source "ensembl_havana"; gene_biotype "protein_coding";
1       1211339 1214153 ENSG00000186827 .       -       ensembl_havana  transcript      .       gene_id "ENSG00000186827"; gene_version "11"; transcript_id "ENST00000379236"; transcript_version "4"; gene_name "TNFRSF4"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "TNFRSF4-201"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS11"; tag "basic"; transcript_support_level "1 (assigned to previous version 3)";
1       1213982 1214153 ENSG00000186827 .       -       ensembl_havana  exon    .       gene_id "ENSG00000186827"; gene_version "11"; transcript_id "ENST00000379236"; transcript_version "4"; exon_number "1"; gene_name "TNFRSF4"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "TNFRSF4-201"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS11"; exon_id "ENSE00001832731"; exon_version "2"; tag "basic"; transcript_support_level "1 (assigned to previous version 3)";
1       1213982 1214127 ENSG00000186827 .       -       ensembl_havana  CDS     0       gene_id "ENSG00000186827"; gene_version "11"; transcript_id "ENST00000379236"; transcript_version "4"; exon_number "1"; gene_name "TNFRSF4"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "TNFRSF4-201"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS11"; protein_id "ENSP00000368538"; protein_version "3"; tag "basic"; transcript_support_level "1 (assigned to previous version 3)";
1       1214124 1214127 ENSG00000186827 .       -       ensembl_havana  start_codon     0       gene_id "ENSG00000186827"; gene_version "11"; transcript_id "ENST00000379236"; transcript_version "4"; exon_number "1"; gene_name "TNFRSF4"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "TNFRSF4-201"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS11"; tag "basic"; transcript_support_level "1 (assigned to previous version 3)";

Finishes without any of the mentioned warnings.

Thanks! I did not know I should have run it as follows: gtf2bed < Homo_sapiens.GRCh38.106.gtf > Homo_sapiens.GRCh38.106.gtf.bed


text:
Converting genetic matrix to STRUCTURE
I want to convert a genetic matrix (rows are samples, columns are alleles) into a STRUCTURE file (output .str) within R. Is there a preferred way to achieve this? I've found genind_to_structure in the graph4lg R package, but I wanted to see if there are other approaches.
An example of what I mean by a "genetic matrix":
# Build a matrix of alleles
temp <- lapply(1:30, function(i) sample(1:9, 4, replace = TRUE))
temp <- sapply(temp, paste, collapse="")
temp <- matrix(temp, nrow = 10, dimnames = list(paste("ind",1:10), paste("loc",1:3)))
temp

# Convert matrix to genind format (for adegenet package)
test <- df2genind(temp, ploidy = 4, sep = "")
test

# How to convert to STRUCTURE?

If anyone has input, please chime in! Thanks

Structure is a bit old now. There is fastStructure (2014) http://rajanil.github.io/fastStructure/ which is Python. Happy to do this if there's a Pythonic Bio.popgen equivalent geared for fastStructure

I’m not totally sure what you’re asking, since the current method seems to be working? Is there a reason why you need another approach?

I ultimately found this function, which I used for the purpose of converting a genind to a STRUCTURE file. The graph4lg genind_2_structure command would regularly error on me (present working directory couldn't be passed as the output argument).


text:
DeconRNASeq: Extract gene names from returned mixing proportions
This question was also asked on Biostars
I am using the Bioconductor package "DeconRNASeq" to perform tissue deconvolution. Let's say I run the following code (this is from the manual):
BiocManager::install(c("DeconRNASeq"))
library(DeconRNASeq)

data(multi_tissue)
datasets <- x.data[,2:11]
signatures <- x.signature.filtered.optimal[,2:6]
proportions <- fraction

res <- DeconRNASeq(datasets, signatures, proportions, checksig=FALSE, known.prop = TRUE, use.scale = TRUE, fig = TRUE).
res.out <- as.data.frame(res[["out.all"]])
str(res.out)

And receive the following results:
'data.frame':   10 obs. of  5 variables:
 $ brain : num  0.0468 0.0718 0.0685 0.0424 0.6471 ...
 $ muscle: num  0 0.1172 0.6081 0.0489 0.1344 ...
 $ lung  : num  0.0939 0.0829 0.1265 0.6147 0.0709 ...
 $ liver : num  0.0947 0.6187 0.1355 0.1897 0.0812 ...
 $ heart : num  0.7645 0.1093 0.0615 0.1043 0.0665 ...

Is it possible for me to obtain the subset of genes that were successfully deconvolved? For instance, the first sample had 0.0468 proportion from "brain". Is it possible for me to see which genes specifically were composing this 0.0468 proportion?
Note: This R package does not appear to have a Q&A forum. If there is a recommended place to post this question other than Bioinformatics StackExchange, I would be grateful to hear suggestions.
Thank you.