text:
How can I create a combined barplot like this image?
I have data from different countries from two groups of samples (control and test). And, calculated three different indices alpha, beta and gamma. I have a large dataframe containing all the data. Here is a reproducible example of my dataframe (df):
df <- structure(list(
    SampleID = c("SRR490_profile", "SRR467268_profile",
                 "SRR467265_profile", "SRR467264_profile",
                 "SRR467263_profile", "SRR467258_profile",
                 "SRR467250_profile", "SRR467244_profile",
                 "SRR467242_profile", "SRR467227_profile",
                 "ERR710432_profile", "ERR710430_profile",
                 "ERR710429_profile", "SRR232646_profile",
                 "SRR232622_profile", "SRR232619_profile",
                 "SRR232617_profile", "SRR232616_profile",
                 "SRR232614_profile", "SRR232612_profile",
                 "HJ4386_profile", "HJ4382_profile",
                 "HJ4380_profile", "HJ4379_profile",
                 "HJ4377_profile", "HJ4375_profile",
                 "HJ4374_profile", "HJ4373_profile",
                 "HJ4372_profile", "HJ4370_profile",
                 "KLM3246730_profile", "KLM3246729_profile",
                 "KLM3246725_profile", "KLM3246722_profile",
                 "KLM3246721_profile", "KLM3246714_profile",
                 "KLM3246713_profile", "KLM3246712_profile",
                 "KLM3246711_profile", "KLM3246705_profile"),
    type = c("test", "control", "control", "control",
             "control", "test", "test", "test", "control",
             "control", "control", "test", "control",
             "control", "control", "control", "control",
             "test", "control", "test", "control", "control",
             "control", "control", "control", "test",
             "control", "test", "test", "control", "control",
             "test", "test", "test", "control", "test",
             "control", "test", "test", "test"),
    variable = c("gamma", "alpha", "beta", "gamma", "alpha", "beta",
                 "gamma", "alpha", "beta", "gamma", "alpha", "alpha",
                 "beta", "gamma", "alpha", "beta", "gamma", "alpha",
                 "beta", "gamma", "alpha", "alpha", "alpha", "beta",
                 "beta", "gamma", "alpha", "beta", "gamma", "alpha",
                 "beta", "gamma", "alpha", "beta", "gamma", "alpha",
                 "beta", "gamma", "alpha", "beta"),
    value = c(95, 90, 109, 112, 128, 109, 90, 87, 90, 99, 129, 101, 113,
              100, 111, 102, 99, 105, 85, 94, 133, 117, 116, 83, 113, 126,
              90, 128, 107, 121, 88, 109, 82, 42, 67, 100, 88, 101, 127,
              115),
    StudyID = c("swiss", "swiss", "swiss", "swiss", "swiss", "swiss",
                "swiss", "swiss", "swiss", "swiss", "australian",
                "australian", "australian", "australian", "australian",
                "australian", "australian", "australian", "australian",
                "australian", "thai", "thai", "thai", "thai", "thai",
                "thai", "thai", "thai", "thai", "thai", "english",
                "english", "english", "english", "english", "english",
                "english", "english", "english", "english")),
    row.names = c(NA, -40L),
    class = c("tbl_df", "tbl", "data.frame"))

Now, I want an image for each index containing barplots for two groups (control and test) and Y-axis showing the StudyID names (i.e., country names) like the following (this image shows plots for two indices alpha and beta):

Is this possible from my data? Can anyone give me a solution?
Thanks in advance

The ggplot2 function stat_summary or stat_summary_bin should get you most of the way there, something like this [not tested]:
library(tidyverse)
library(ggplot2)

df %>%
  ggplot() +
  aes(x=value, y=StudyID, colour=type) +
  stat_summary() +
  facet_wrap(~variable)

See ?stat_summary for more information.

If you're not wedded to a dot and whiskers plot as the one show, I'd probably plot this dataset as a boxplot with ggplot:
ggplot(df,aes(x=value,y=StudyID,colour=type)) +
geom_boxplot() +
facet_wrap(~variable)


text:
Is SBOL-VOWS operational?
Synthetic Biology Open Language's Visual Ontology Webservice docs have example URLs but they all return a github pages 404.
EDIT:
The SBOL-VOWS has site has been fixed but the section on the SBOL website still has links with a 404.

would you mind adding the link in the docs that returns 404 github page? And perhaps mentioning what did you expect to see there instead? It just so people will be more easily understand/find your question if then have the same problem.

@KamilSJaron ah my bad. I just realized that SBOL website itself still has broken links!

Many thanks for raising this issue. The web service is now working. Only the glyph related web service calls were affected due to restructuring folders at the SBOL-Visual GitHub repository. I now revived the glyphs using the sbol-visual-ontology GitHub repository. Over the summer, we will develop the web service further so that you can request different versions of glyphs. Currently, the web service works with the SBOL Visual v2.3 glyphs and we will provide support for SBOL Visual v3 glyphs soon.
The web service is available at:
https://vows.sbolstandard.org/
Please let me know if there is anything I can clarify.

Hello - I've just noticed the links on the SBOL website are still broken.

source: Stack Exchange

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

sort: votes


text:
polidy VCF file for Canvas somatic-wgs cnv calling
I am trying to run Illumina CANVAS cnv caller for Somatic-WGS. There is an option "--ploidy-vcf" which is mandatory to supply, but don't know what exactly that mean. I had supplied the CNV.vcf file from the germline WGS (matched normal run on Canvas smallPedigree), but it has given an error that says "Reference copy number > 2 is not supported" at the last stage of run.
I would like to know your experience when you used CANVAS on somatic WGS samples. What ploidy-file is appropriate here? Should it be only chrX and Y in it? Your suggestions please.

It seems you need to supply vcf files specifying the ploidy of the sex chromosomes.
See https://github.com/Illumina/canvas/issues/89 for explanation and example vcf files.

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


text:
What is the right way of calculating a Phred score by hand?
i am trying to calculate mean Phred scores for my sequencing data, but i feel not very comfortable about it. There are actually two ways of calculating. (I just use an existing sample)
giving: 3 reads with a mean phred score of 20, 10 and 3 (meaning Error rates of 1%, 10% and 50%, respectively)
Way 1 to calculate: Just take the arithmetic mean of the q scores:
(20+10+3)/3 = 11
Way 2:  Take the arithmetic mean of the error rates:
(0,5+0,1+0,01)/3 ~= 0.2
And then calculate the q score by
-10*log10(0.2) = 6.99
(https://gigabaseorgigabyte.wordpress.com/2017/06/26/averaging-basecall-quality-scores-the-right-way/)
The Author of this Blog claims, that Way 2 is the correct version, but i still do not get why.
Is somebody able to answer this question?
Thanks a lot

Not sure if my explanation is any good but let's try...
First, let's convince ourselves that converting phred to probability is the right thing to do (i.e. your Way 2).
Imagine a read of length 10000 bp with half of the bases with Q=10 (or P=0.1) and the other half with Q=20 (or P=0.01). The number of wrong calls you expect is therefore 550:
n = 10000
(n * 0.01 * 0.5) + (n * 0.1 * 0.5) # -> 550 wrong calls

So the overall error rate is 550/10000 = 0.055 which corresponds to Q=12.6. If you follow your Way 2 you do get this answer:
((n * 0.5 * 0.01) + (n * 0.5 * 0.1))/n = 0.055
-10*log10(0.055) = 12.6

If you follow your Way 1 you get Q=15, which means an overoptimistic error rate of ~0.03 or ~316 wrong calls. Why is that?
Because phred is on a non-linear, logarithmic, scale: a unit change  towards the high values is less important than a unit change towards the lower end. For example, across 10000 bases a change from Q=30 to Q=31 means a change from 10 to ~8 wrong calls - a difference of 2 only. But a change from  Q=10 to Q=11 means 1000 to ~800 wrong calls: the difference is a substantial ~200 calls.
Having said this, the straight average of phred scores is not "wrong". It measures the overall phred score but it does not reflect the expected number of wrong calls, which is probably what you want in this case.

Thanks alot for this clear answer! This just means, that for calculation of the expected wrong calls it is necessary to have a look first on the sequencing technique. If it is expected to have low/medium output reads (like a not that good MinION sequencing) it is necessary to go with way 2. For higher quality reads (like most of Illumina reads), both methods can be considered, but just to be safe, way 2 is the right one to choose. Thanks alot for the explanation


text:
How to introduce point mutations in a protein structure?
I want to introduce point mutations to a protein PDB file and automatically predict the 3D coordinates of a new amino acid residue in the mutated protein.
I used interactive Schrodinger's Maestro tool to mutate residues and minimize them to avoid clashes, however this time i want to automate the task and use some open source tools.
Is there any free software for this task?

PDBFixer and pyMutate are simple, lightweight, free and open-source. Both don't do any optimizing, though, so you may have to run some molecular dynamics afterward to get what would be the right "orientation".
VMD would work for sure, I believe PyMOL too. They are both graphical interfaces but you can start them in text only mode (or just import pymol in a python script and use the Python API). They are free for particulars and academics, but given their renown you may already be using them.
If you're looking for something a tad more complete and that may adapt the side-chain to the surrounding environment I'd go with Modeller, which provides two example scripts to do point mutations (here and there).
If you don't want to just get the side-chain of this amino-acid right but actually strive to see some general backbone adaptation too then use PyRosetta. It can probably do what all the others do, too, but given it is 1 GB big it may be overkill. Any of the package coupled with molecular dynamics would also do the trick.

Dear @O.Laprevote, thank you for your valuable answer! I actually tried running Rosetta's fixbb tool to mutate the truncated pdb residue with all missing atoms but one, however apparently it just ignores such residues, while mutating the truncated residues was my very task.

Modeller did the job very well for me. I am accepting this answer!


text:
What does presence/absence of black lines mean in IGV?

IGV (Integrated Genome Viewer) is a popular open-source tool for viewing alignment files. In my BAM file in IGV, some deletions have black lines running through them, and others don't. What causes this? Short black lines also appear adjacent to some mutations and insertions.
I generated the CIGAR for the SAM file myself, and am wondering if I made an error. I didn't find anything in IGV's documentation about this.

What do you see when you hover over these black lines using your mouse? Hovering on a few deletions with/without the black lines might give you an idea on the differences between them.

It doesn't show anything besides the general read metadata. Only insertions show the inserted bases on hover. Insertions with/without short black lines through them are indistinguishable on hover.

Those look like some single-T reads (i.e., reads with a single base). What do you see when you hover over them? Do you get a more coherent bit of info? What does the CIGAR look like for one of those?

They look like short reads, but it's actually a deletion within a long read that for some reason isn't showing a black line. I'm wondering what metadata difference in the SAM/BAM could possible cause the omission of the black line...

Black line represents an indel.  Read was lengthened to fit consensus.
No black line (blank), no read/sequence.  Read likely ended or was truncated.

In general, this is true. However, for my dataset, these are all long reads which I know span the region I'm looking at. For some reason, the black lines (which should be there) are not present for some INDELs.


text:
How to execute several processes sequentially, without channels?
Using DSL2.
I have processes A, B, C which I want to run in sequence - B waits for A to finish, C waits for B to finish.
They are not communicating via channels so I can't do A.out.collect().
I have tried making dummy channels but this is problematic and a hack. The onComplete subscription doesn't seem to let me trigger another process.
What is the correct way to run a sequence of processes which are not communicating or triggered by channels?
P.S. I am aware that this is not really the intended use case of Nextflow (which is oriented around channel communication/triggers... but... I just need to get this done.)

Could you please elaborate on how you tried to do dummy channels? Maybe by giving an example. I would suggest to use watchPath but using this type wouldn't allow your pipeline to finish.

I would guess it is not possible without channels unless you are using beforeScript and afterScript. I came up with this solution:
#!/usr/bin/env nextflow

nextflow.enable.dsl=2

workflow {
  // C(B(A("")))

  next = ""
  next = A(next)
  next = B(next)
  next = C(next)
}

process A{
  echo true
  input: val _
  output: val ""

  """
  echo $task.process
  """
}

process B{
  echo true
  input: val _
  output: val ""

  """
  echo $task.process
  """
}

process C{
  echo true
  input: val _
  output: val ""

  """
  echo $task.process
  """
}

By using workflows you can change the execution order. But it is still a hack and I would prefer something like a next: keyword to specify the next process.

And for those not interested in DSL2:
process A {
  echo true

  output:
  val "" into B

  """
  echo "$task.process"
  """
}

process B {
  echo true

  input:
  val _ from B

  output:
  val "" into C

  """
  echo "$task.process"
  """
}

process C {
  echo true

  input:
  val _ from C

  """
  echo "$task.process"
  """
}

Also, I am quite curious to see other solutions.


text:
Making Seurat object from a processed file
I have downloaded log2(TPM/10+1) values of 11,548 genes and 9609 cells from GSE146026 (10x) in tsv format as the raw data is not available
I see patient IDs, cell barcodes, genes, even assigned clusters are here
> dput(head(a[1:5,1:5]))
structure(list(Cell_ID = c("10x_barcode", "patient", "time",
"sample_ID", "clst"), X10x_1 = c("10x_3288_t1_AAACATACCTTCCG-1",
"5", "1", "3288.1", "1"), X10x_2 = c("10x_3288_t1_AAACATACTCCTAT-1",
"5", "1", "3288.1", "1"), X10x_3 = c("10x_3288_t1_AAACATTGAACTGC-1",
"5", "1", "3288.1", "1"), X10x_4 = c("10x_3288_t1_AAACATTGCTGACA-1",
"5", "1", "3288.1", "2")), row.names = c(NA, 5L), class = "data.frame")
>

I want to make a Seurat object of that but really I don't know how to do that
I have tried this
    CreateSeuratObject(getwd(),gene.column = 2,cell.column = 1,unique.features = TRUE,strip.suffix = FALSE)

But get the following error:
Error in CreateAssayObject(counts = counts, min.cells = min.cells, min.features = min.features) :
  No cell names (colnames) names present in the input matrix

Anybody have ever dealt with such a case to help me please?

The function you need is CreateSeuratObject() and not Read10X() as you start from TPMs.


text:
Is SNPRelate the goto package for PCA of SNP data in R?
I'm trying to evaluate a method that I'm not familiar with. They used SNPRelate. I'm wondering if this is the best / only choice... what other tools would you recommend?
Many thanks,

Does it have to be in R? If not, the standard things to make a PCA from genotype matrix is plink or eigenstrat.

That said, unless you are doing fancy stuff like PCA projection - PCA is just PCA and whilst methods might differ in terms of computational compexity, the resulting eigenvectors will be mostly very similar.

SNPRelate claims nearly 10x faster than plink... but I don't know.


text:
Convert SNP data from GEO into PED and MAP files?
I can easily grab SNP-Chip data from NCBI's GEO database for a given Illumina SNP-Chip array, e.g. from here:
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52147
GEO provides data as SOFT formatted family file(s), MINiML formatted family file(s) or Series Matrix File(s).
I'm looking for a way to convert these into PED and MAP files for use in PLINK...
Or does PLINK read these formats natively?

I see. Maybe PLINK cannot handle this. I see on GEO that series downloads are suitable for opening in spreadsheet.

I don't want to manually play with 3000 sheets...

Sadly the 'contact us' email at GEO didn't help at all.
I ended up probably re-inventing the wheel with this script:
https://gist.github.com/CholoTook/60968e3ab6d90cb8fd19be55a25592f1
YMMV


text:
Retain only part of a file name and fasta header in fasta directory
I have a MainDirectory/hundreds_of_subdirectories/thousands_of_fastas structure. For each fasta in the bottom level, I am trying to change the file name as well as the first line in the file to retain only some parts of the file.
So that the following text, which makes up both the file name and the first name of the file:
947-length-1150-cov-1000|contig:JAECWU010000213.1|slice:817050-818200|uce:uce-1452|match:817550-817700|orient:{'+'}|probes:4.unaligned.fasta

becomes just:
uce-1452.unaligned.fasta

For editing the first line in the file, I have tried
sed '1 s/^[^|uce-]*(|uce-[0-9]).*/\1/' hundreds_of_subdirectories/*

but this isn't the ticket..
Then messing around with some kind of loop...
for i in *\*.unaligned.fasta
do
  sed -E 's/^[^|uce-]*(|uce-[0-9]).*/\1/'
done

For the file name, it was suggested that I use rename instead of sed, but with sed I have tried something like:
for x in hundreds_of_subdirectories/thousands_of_fastas*unaligned.fasta
do
  echo $x | sed -r 's/^[^|uce-]*(|uce-[0-9]).*/mv' -v "\0" "|uce-[0-9].*/\1/.unaligned.fasta"/
done

Am I even barking up the right tree? Thanks!

What operating system are you using?

I think you're working too hard to capture the uce: expression.
Try
sed -Ee '1s/^.+uce:(uce-[0-9]+).+/\1.unaligned.fasta/' file_glob

to replace the first line content. You may want to add the in-place flag after testing the expression on a few files.
For the mv, use command substitution, which runs the command in a sub-shell and uses the output value as the parameter:
for f in file_glob
do
    mv -- "$f" "$(echo \"$f\" | sed '...')"
done

Note that your sed will not do what the OP wants: it will change the first line to >uce-1452 and not uce-1452.unaligned.fasta.

Thanks, @terdon. I was nudging OP in the right direction. The .unaligned.fasta is trivial to add once OP switches to a simpler solving mindset. I'll add it in to make the solution more complete though.

Force of habit: that protects from the (admittedly unlikely) case where a fasta file's name starts with -. Without the --, mv will treat that as an option and fail. Try  mv -badName goodName and contrast with mv -- -badName goodName.

Thanks, I did some googling and saw why the -- prevents filenames from manifesting nightmares in worst case scenarios, which is why I deleted my earlier comment that asked you why you recommend mv --.

Of course, if I were being careful, I should also have added -- to the echo in case the file name begins with -n. But I forgot and it really is vanishingly unlikely that this would happen with actual fasta file names.

@terdon - Only <() and >() are process substitution. $() is not that. It is a different invocation but definitely not process substitution as it does not replace a file. EDIT: I think it's called command substitution. "Subshell" is closer than "process substitution" though, so please edit judiciously.

Argh! You're absolutely right. I meant command substitution. Sorry, I should indeed be more careful. I just wanted to remove the "subshell" since that is misleading.

Let us continue this discussion in chat.

Brilliant, now to try the above. Thank you!

Thanks for the suggestions above @terdon. However, both suggestions (if I run from the main directory, take the .fasta files out of the subdirectories, correctly renamed, but moved from their subdirectory. This causes them to overwrite each other, as each subdirectory may have fastas with the same uce-number assigned. I ran: for file in */.fasta;do
perl -i.bak -pe 's/^\s>.uce:(uce-\d+).(.unaligned.fasta)/$1$2/' "$file"
rename 's/^\s>.uce:(uce-\d+).(.unaligned.fasta)/$1$2/' "$file"
done

If you are using bash, you can enable the globstar option (shopt -s globstar) which allows ** to match recursively. That way, **/*fasta will find all files whose name ends with fastq in this directory and all of its subdirectories. So that will let you easily iterate over all the target names. Next, if you have the perl rename command (called rename on Debian and Ubuntu and similar systems, prename or perl-rename on others), you can do it quite easily.
Note that perl-rename has the -n option that makes it print out the renaming it would do without actually renaming anything. So I suggest you first run this to test what you will be doing:
shopt -s globstar
for file in **/*.fasta; do
   echo perl -i.bak -pe 's/^\s*>.*uce:(uce-\d+).*(\.unaligned.fasta)/$1$2/' "$file"
   rename -n 's/^.*uce:(uce-\d+).*(\.unaligned.fasta)/$1$2/' "$file"
done

That will echo the perl command, so you can try it manually first and make sure it works, and will print out the rename operations without actually doing anything. If you are satisfied it will work, run again without the echo or -n:
for file in **/*fasta; do
    perl -i.bak -pe 's/^\s*>.*uce:(uce-\d+).*(\.unaligned.fasta)/$1$2/' "$file
    rename 's/^.*uce:(uce-\d+).*(\.unaligned.fasta)/$1$2/' "$file"
done

Or, if you don't have perl-rename, you can use the shell's string manipulation for this:
for file in **/*.fasta; do
    perl -i.bak -pe 's/^\s*>.*uce:(uce-\d+).*(\.unaligned.fasta)/$1$2/' "$file
    tmpFileName="${file#*uce:}"
    mv -- "$file" "${tmpFileName/|*/.unaligned.fasta}"
done

Explanation

shopt -s globstar: this enables the globstar bash option for the current shell session. It only needs to be run once, its effects will last until you close the terminal or start a new shell session. As for what it does, see man bash:

globstar
If set, the pattern ** used in a pathname expansion context will match all files and zero or more directories and subdirectories. If the pattern is followed by a /, only directories and subdirectories match.

for file in **/*fasta; do ... ; done: iterate over all files (and directories) whose name ends in .fasta, saving each file name as $file.

perl -i.bak -pe 's/^\s*>.*(uce-\d+).*(\.unaligned.fasta)/$1$2/' "$file: this changes the fasta header. The -i.bak option tells perl to edit the file in place, but to save the original file with a .bak extension. This can help you revert any changes if there is a mistake. Once you have finished, you could run rm **/*fasta.bak to delete the backups. Alternatively, just remove the -i.bak if you don't want to keep backups.
The -p option tells perl to process its input line by line and print every line after applying whatever script you give it with -e.
The s/old/new/ is the substitution operator, same as sed. The regular expression is looking for 0 or more whitespace characters and then a > at the beginning of the line (^\s*>). It then allows for 0 or more characters (.*) until it finds the string uce: followed by uce- and one or more digits (uce:(uce-\d+)). The parentheses "capture" the pattern, storing it as $1. Next, we want 0 or more characters again until we reach .unaligned.fasta (.*(\.unaligned.fasta)). We capture the .unaligned.fasta as $2 (second set of parentheses) only to avoid typing it out again. Note that this assumes you will never have anything after the .unaligned.fasta. Finally, we replace everything matched above with just $1 (the uce-numbers) and $2 (.unaligned.fasta).

rename 's/^.*uce:(uce-\d+).*(\.unaligned.fasta)/$1$2/' "$file": the rename opration. This will match everything from the beginning of the file name until the last occurrence of uce:uce-NNNN, capture the uce-NNNN as $1 then skip the remaining characters until .unaligned.fasta and finally rename as just $1$2 (the capturing logic is the same as above).

tmpFileName="${file#*uce:}": The ${var#pattern} syntax means "print the value of $var, but remove the shortest match of pattern from the front of the variable's value". The ${var//pattern/replacement} syntax means "print the value of $var and replace the first occurrence of pattern with replacement. Here, the pattern is |* which means "anything after a |" (note that this is shell stuff and is using globs and not regular expressions). This is probably easier to understand with an example:
$ echo "$file"
947-length-1150-cov-1000|contig:JAECWU010000213.1|slice:817050-818200|uce:uce-1452|match:817550-817700|orient:{'+'}|probes:4.unaligned.fasta
$ tmpFileName="${file#*uce:}"
$ echo "$tmpFileName"
uce-1452|match:817550-817700|orient:{'+'}|probes:4.unaligned.fasta
$ echo "${tmpFileName/|*/.unaligned.fasta}"
uce-1452.unaligned.fasta

So, given tmpFileName="${file#*uce:}", "${tmpFileName/|*/.unaligned.fasta}" is the desired file name and this is what we pass to mv.


text:
use same output in two processes in nextflow dsl2
This is my workflow:
pre_align()
pre_align.out.single_fastqs.view()
get_fq_info(pre_align.out.single_fastqs)

align_bwa(get_fq_info.out.fq_info)

align_bowtie2(get_fq_info.out.fq_info)

where I want to use the same output from get_fq_info as input for either align_bwa and align_bowtie2`.
As it's shown here only align_bwa is executed. If I comment align_bwa THEN align_bowtie2 gets executed.
How can I fix this?
Thank you!

I have solved this problem in the past by duplicating the output channel of the desired process in order to feed the same output into multiple downstream processes.

As far as I know, DSL2 should be able to handle this. Have you specified an output qualifier for your process get_fq_info (e.g. path)? I ran into a similar issue when I forgot to specify it.
For example, this will result in the error process align_bowtie2 declares 1 input channel but 0 were specified:
process get_fq_info {
   input:
   path reads

   output:
   "output.txt"

   shell:
   """
   myprocess.sh ${reads} > output.txt
   """
}

And is probably solved by this:
process get_fq_info {
   input:
   path reads

   output:
   path "output.txt"

   shell:
   """
   myprocess.sh ${reads} > output.txt
   """

```


text:
How does FDRtool work?
I have a question about using FDRtool. In the below code (on RNA seq data whose p values were acquired using Deseq2), the FDRtool was first used and thereafter p.adjust using the benjamini hochberg method. I didn't write this, its taken from someone else's code found online, in addition it is also the way the it is written in the pipeline at my university. Therefore, I assumed that this is how fdrtool is used and in my eyes it looks like they both adjust the p value with the same method. So that's why I asked, as I want to understand why you would need both these lines, which I've seen many people do:
FDR.ddsRes <- fdrtool(ddsRes$stat, statistic= "normal", plot = T)
ddsRes[,"padj"]  <- p.adjust(FDR.ddsRes$pval, method = "BH")

Dont they both correct for false discovery rate? Why do you need to run p.adjust after using fdrtool?
Thanks!!

Something like this book chapter might help on understanding FDR. The packages or functions you call differ in their estimation of pi0 or the proportion of null hypothesis.
Basically, the p.adjust(..method = "BH") is a more conservative method with Benjamini-Hochberg, you are assuming the proportion of null features to be 1.
When you use fdrtool, you are estimating the proportion of null hypothesis and from there, the local FDR. This in theory gives you more power.
However, I have to note that in your code, you are not using fdrtool to estimate the fdr as explained above, because:
FDR.ddsRes <- fdrtool(ddsRes$stat, statistic= "normal", plot = FALSE)

With this, the raw p-values are stored in FDR.ddsRes$pval and the q value, which is the corrected p-value, is stored in FDR.ddsRes$qval. In your code you are basically taking the raw p-values, and applying Benjamin-Hochberg:
p.adjust(FDR.ddsRes$pval, method = "BH")

You should apply a two tailed test. From the wald stats, it should be:
derived_p = 2*pnorm(-abs(ddsRes$stat))

In your code, you are running it wrongly by passing the test statistic into fdrtools, we can compare the raw p-values:
head(data.frame(stat = ddsRes$stat,derived = derived_p,
DESeq2=ddsRes$pvalue,fdrtools=FDR.ddsRes$pval))

        stat    derived     DESeq2  fdrtools
1  1.2405796 0.21476110 0.21476110 0.1216358
2  1.2314421 0.21815755 0.21815755 0.1244054
3 -0.3292516 0.74196552 0.74196552 0.6811997
4 -0.4964574 0.61957174 0.61957174 0.5356140
5 -0.4101163 0.68172066 0.68172066 0.6088427
6 -1.8007943 0.07173531 0.07173531 0.0246428

To summarize, use the adjusted p-value returned by DESeq2. The code you have above is wrong unless there's a good reason to do a one-tail test.


text:
How to index residues by icode in biopython?
I am indexing a structure by the code below:
PDB_CODE = '1lve'
RES_CHAIN = 'A'
RES_NUM = 27

structure = PDBParser(is_pqr=True).get_structure(PDB_CODE, PDB_CODE+".pqr")
residue = structure[0][RES_CHAIN][RES_NUM]

However, the 1lve contains residues with the same RES_NUM but different icodes 'A' to 'F' so my code either ignores them or throws an error if I pass '27A'. Does biopython allow this kind of structure indexing directly?

Okay, I found the answer myself. It is possible to index not only by RES_NUM, but by its id (check with residue.get_id()), so the below example would work:
#structure[0][RES_CHAIN][het, resseq, icode]
structure[0][RES_CHAIN][' ', 27, 'A']


text:
RDAVIDWebService (R package) manual installation, or alternative tool
Recently, I had to re-install all my R packages (R version 4.1.0 (2021-05-18) -- "Camp Pontanezen") and it turns out that "RDAVIDWebService" is not available anymore on Bioconductor:

package ‘RDAVIDWebService’ is not available for Bioconductor version '3.13'

Is there an another way to install it (couldn't find the files on the Bioconductor site to install the package manually) or is there another R package that allows to use the DAVID tool?

That package was deprecated in Bioconductor version 3.13 after several attempts to contact the package author/maintainer.
Your options include reverting back to Bioconductor version 3.12 or using one of the other options for GO/pathway enrichment analyses, e.g. topGO, enrichR, or clusterProfiler.


text:
Ccle miRna data units
Does anybody know what units the mirna expression data on ccle is in? I saw the pipeline in the original paper but it doesn't mention units.
Original paper: https://www.nature.com/articles/s41586-019-1186-3?proof=t
Data: https://portals.broadinstitute.org/ccle/data


text:
Python package or CLI tool outputting mutations of sequence with respect to reference
I'm new to bioinformatics, and I'm looking for a tool that takes two FASTA files, one containing a reference, let's call it A.fasta, the other another sequence, let's call it B.fasta, and outputs just the mutations that B contains with respect to A.
I'm open as to the output format. Is there a standard for this?
Background: I'm trying to analyse Sars-CoV-2 sequences and identify mutations with respect to the reference. I've done a google search for a tool that does what I want to but couldn't find something straightaway. Is there no standard tool for this?

I would say parsimony, but thats just me

If it's just a comparison of two sequences that you're looking for, then diffseq from the emboss toolkit should be okay:
https://www.bioinformatics.nl/cgi-bin/emboss/diffseq
I think this tool requires the sequences to be oriented in the same direction. If you discover a lot more variants than you expect, try reverse-complementing one of the sequences.


text:
Construct the Overlap Graph of a Collection of k-mers
I am doing a course on graph algorithms in genome sequencing and the current assignment involves building an overlap graph from a bunch of k-mers. Formally the problem is as follows:
Input: A collection of Patterns of k-mers.
Output: The overlap graph in the form of an adjacency list.
I have a working algorithm that handles most cases except for self-loops. The problem is that I don't know when/why I should be adding self-loops. Any help on this would be greatly appreciated. Below are some examples illustrating my question and my current solution. Thank you!
Example 1:
Input:
ACT
CTT
TTT
Output:
TTT->TTT
ACT->CTT
CTT->TTT
Example 2:
Input:
CCCC
Output:
CCCC->CCCC
Example 3:
Input:
CT
TG
TG
TC
TT
TC
Output:
TC->CT
CT->TC,TG,TT
TT->TC,TG,TT
The goal here is to represent a collection of kmers as an overlap graph where unique kmer will be represented as a node in the graph. Directed edges will connect nodes whose last $k-1$ characters overlap with the other node's first $k-1$ characters. Therefore, repeated kmers that overlap with a other kmer should be reflected as self loops. That's what I understood, however I am looking for a more general rule/solution.
Here is another example to illustrate this.
Input:
CT
TT
TT
TT
TT
TT

Output:
CT->TT
TT->TT

Current code:
/*
 * Creates an overlap graph from a list of k-mers
 *
 * \param patterns      a vector of string k-mers
 *
 * \return a string containing an adjacency list representation of the
 * overlap graph as described in the problem specification
 *
 * Output Format. The overlap graph Overlap(Patterns), in the form of an
 * adjacency list. Specifically, each line of the output represents all
 * edges (u,v) leaving a given node u in the format “u -> v”, where u and v
 * are both k-mers. If a given node u has multiple edges leaving it (e.g. v and w),
 * the destination nodes are comma-separated in any order. For example, if
 * there exist nodes ACG, CGT, and CGG, the resulting line of the adjacency
 * list would be: ACG -> CGT,CGG
 */

#include <algorithm>
#include <iostream>
#include <map>
#include <vector>
#include <set>
#include <string>

std::string overlap_graph(const std::vector<std::string>& patterns) {
    // Overlap Graph of the k-mers collection patterns, in which a vertex exists for
    // each unique k-mer, and an edge (u,v) exists if the last k-1 characters of u are
    // equal to the first k-1 characters of v.

    if (patterns.size() == 1) return patterns[0] + "->" + patterns[0];
    // Build adjancency list O(n^2)
    std::map <std::string, std::vector<std::string>> adj;
    for (const auto& kmer_u : patterns) {
        // Add kmer as key if it doesn't exist in the map and find its neighbours.
        if (adj.find(kmer_u) == adj.end()) {
            adj[kmer_u] = {};
            for (const auto& kmer_v : patterns) {
                if ((kmer_u != kmer_v) &&
                    (kmer_u.substr(1, kmer_u.length() - 1) == kmer_v.substr(0, kmer_v.length() - 1)) &&
                    (std::find(adj[kmer_u].begin(), adj[kmer_u].end(), kmer_v) == adj[kmer_u].end())) {
                    adj[kmer_u].push_back(kmer_v);
                }
            }
        } // Otherwise continue.
    }

    // Convert to string O(n).
    std::string ans;
    for (auto const& key_value : adj) {
        std::string kmer_u{ key_value.first };
        std::vector<std::string> kmers_v{ key_value.second };
        if (!kmers_v.empty()) {
            ans.append(kmer_u + "->" + kmers_v[0]);
            for (std::size_t kmer_v{ 1 }; kmer_v < kmers_v.size(); kmer_v++) {
                ans.append(',' + kmers_v[kmer_v]);
            }
            ans.push_back('\n');
        }
    }
    return ans;
}

void tests() {
    /* Sample Dataset */
    std::vector<std::string> patterns{ {"AAG"}, {"AGA"}, {"ATT"}, {"CTA"}, {"CTC"}, {"GAT"},
                                       {"TAC"}, {"TCT"}, {"TCT"}, {"TTC"} };
    std::string answer{ "TTC->TCT\nCTA->TAC\nAAG->AGA\nCTC->TCT\nTCT->CTA,CTC\nAGA->GAT\nGAT->ATT\nATT->TTC\n" };
    std::cout << answer << std::endl;
    // std::cout << (overlap_graph(patterns) == answer) << std::endl;
    std::cout << overlap_graph(patterns) << std::endl;
}

void solution() {
    std::vector<std::string> patterns;

    std::string pattern;
    while (std::cin >> pattern) {
        patterns.push_back(pattern);
    }
    std::cout << overlap_graph(patterns) << std::endl;
}

int main() {
    //tests();
    solution();
}

I think this is actually a de-bruijn as opposed to an overlap graph. I believe an overlap graph has an edge whenever any positive length suffix of a node overlaps with the prefix of another node.

Also, since you are only concerned with length $k-1$ overlaps, you can find all of the edges in $O(n)$ time, since each node can only have 4 possible outgoing edges. I can extend my answer to include details on this if interested

@Throckmorton Thank you for the feedback - I really appreciate it. However, in the de Bruijn graph, the nodes are the $|k-1|$  unique k-mers whereas in the overlap graph the nodes are the unique $|k|$ k-mers. For instance, the $k=3$ de Bruijn graph for genome AGT will be AG->GT and the overlap graph is AGT->AGT (I am actually having troubles with this part - why did we need to add the self-loop here?). I agree with you that the de Bruijn graph can be solved in $O(n)$  as we don't need to check for duplicates, ie. for TTTT and $k=3$ the graph is TT->TT, TT.

You can define $k$ to be whatever you want. The definition of a de-Bruin graph is a graph just all $k$-mers as nodes and an edge whenever there is a $k-1$ overlap between two nodes. This is exactly what the docstring in your question states.

Also, I don’t think there should be a self edge for AGT, in fact inter example the node “ACT” doesn’t have a self edge either.

And the O(n) time isn’t due to duplicates, it’s because for a single node you can find all out edges in constant time.

You will want to create a self loop for a $k$-mer $s$ whenever $s_1,...,s_{k-1} = s_2,...,s_k$. This will only happen when all of the characters of the $k$-mer are identical, as $s_1=s_2, s_2=s_3,...,s_{k-1}=s_k$.

The usual way to avoid loops in linked lists is to keep a record of the things that have already been used / visited, and stop when an already visited thing has been seen again (or when there is nothing more to visit).


text:
Annotating a .vcf with centimorgan information
Some programs (e.g. shapeit4) automatically annotate an INFO tag into a .vcf file which gives the cumulative genetic distance in cM between each SNP:
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##fileDate=10/11/2020 - 21:29:39
##source=G2H
##contig=<ID=9>
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=AC,Number=1,Type=Integer,Description="Allele count">
##INFO=<ID=CM,Number=A,Type=Float,Description="Interpolated cM position">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Phased genotypes">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
##bcftools_viewVersion=1.9-207-g2299ab6+htslib-1.9-271-g6738132
##bcftools_viewCommand=view -s -1_-1 ukb_cal_chr9_v2.HumanOrigins.common_snps.HumanOriginsInds.fiftyPercentEthnicUKBInds.PHASED.bgz.vcf.gz; Date=Mon May 24 14:15:51 2021
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  -1_-1
9   257034  rs573860    A   G   .   .   AC=0;AF=0.388699;CM=0;AN=2  GT  0|0
9   270321  rs17720310  A   G   .   .   AC=0;AF=0.0820595;CM=0.007994;AN=2  GT  0|0

However, this program is really used for phasing not for .vcf annotation.
Ideally I'd like a simple command line program which takes input the .vcf and a genetic map and interpolates the genetic distance and adds a cumulative map. This could be done manually with a custom script, but quite frankly I can't be bothered to do that.
Does anyone know of a tool that exists which will do this?
I posted this to biostars a while ago and received one reply which linked to this post, which includes links to tables containing cM locations for different rs numbers.
EDIT: I am interested in doing this for human data (although the same would be the case for any diploid organism).
Sample vcf:
https://pastebin.com/raw/YSbuU0qM

I think a useful tool to address your question is QCTOOL v2. Particularly, the option you are looking for is -annotate-genetic-map. This tool has been developed originally for GWAS datasets. It has multiple inputs data format.
If you could provide input data as an example, I will try working on it to provide the necessary steps to be done.
Hope this answers your question.

Sorry took a while, have added the sample vcf

@user438383 could you please specify whether your VCF is of hg38 or hg19? Also, do provide a proper VCF format (it is missing some columns like FORMAT followed by another column like in your original post example ). I tried using the file you provided, it is giving an error of segmentation fault. When I tried to add a vcf header with the proper format aforementioned, it works.


text:
Important genes beyond PAM50 for breast cancer classification
Note: this question has also been asked on Biostars
I am currently trying to complete a Breast Cancer Classification task using Neural Networks. I have experimented with using my full dataset of gene expressions(1000+) and managed to improve my results by using only the PAM50 genes. I was wondering though what other genes I could add to the list to even better improve my results.
I have tried looking for papers about extending the PAM50, but could not find anything. I would greatly appreciate any advice.

There is a CSV table in this paper with 33 sets of genes.


text:
RNA_Seq data aligned used uniquely or multi mapped reads impact on result interpretation
I have some transcriptomic (Whole) sequencing data that I should analyse. I would like to do raw data alignment to a reference genome taking into account the multi mapped reads and uniquely mapped reads. I would like to know the consequences on my transcriptomic analysis results : 1. if I should only focus on uniqueMapping reads (mapped to one loci) and neglect MultiMapping read (mapped to several loci or variants)! 2. or take them all together
So my questions are :
1- Considering unique and multi-aligned reads, wouldn't this be a source of bias in my results interpretation?
2- By hiding the multi-mapped reads, could I lose information on the expression of important genes?
thank you in advance

There are several methods/tools that take into account multi-mappers and deal with them in reasonable ways (e.g. pseudo-alignment methods like salmon and more traditional alignment-based methods like RSEM). If you completely ignore multi-mappers, then yes, you will be losing information, which may or may not be valuable to you.
Both of those methods utilize the expectation-maximization algorithm (or a variant thereof) to determine the relative contributions of multi-mapping reads to the various mapping locations. With salmon, you can forego alignment completely and end up with very accurate gene quantifications in hand in like 20 minutes if you want.

I followed your advice and I didn't have a big difference between the results of considering single mapped and multimapped reads. Thanks you Jared Andrews


text:
Smith-Waterman traceback with large gaps
I am currently trying to understand how the traceback algorithm is supposed to work for the smith-waterman algorithm as my current understanding breaks down in case of a large alignment gap.
Assume the sequences
3GA12CA
GA56CA

The obvious best alignment would be (or some other variation thereof)
3GA12-CA
-GA--6CA

However, assuming an affine gap penalty function of x => x + 2 and a binary similarity score (i.e., if equal then 1 else -1) would lead to a scoring matrix
    3 G A 1 2 C A
  0 0 0 0 0 0 0 0
G 0 0 1 0 0 0 0 1
A 0 0 0 2 0 0 0 0
6 0 0 0 0 0 0 0 0
C 0 0 0 0 0 0 1 0
A 0 0 0 1 0 0 0 2

When doing the traceback I would start with the largest value towards the lower right corner of the matrix

(6,6)
(5,5)
Now I do not know how to find the next best alignment (in this example at (3,2)) as neither going up nor going to the left nor going diagonal is going to lead to a score greater 0 meaning I have no next best value to go on.

How is the traceback function supposed to continue here?
Additionally, I read that the traceback should start at the highest score in the matrix itself. However, if we would extend the first subsequence match length by one, this would mean we miss the right alignment option completely.

Before we get started, remember that mismatches are often preferred over gaps, so the "obvious" best alignment would treat one of the 12 in seq1 as a mismatch and the other as a gap, not introduce gaps in both sequences unless your mismatch penalty >> the gap open penalty (which makes no sense)

I'm still trying to understand your problem statement. Affine gap penalties make expanding gaps less expensive than opening new ones, they do not make gaps cheaper than mismatches. "Gaps" are biological sequence insertions/deletions, and those are serious changes, you do understand that, right?

This question is based on a misunderstanding. The Smith-Waterman algorithm calculates the best local segment. You have to start at the $max(H_{i,j})$ point and run the traceback as described. It will (most likely) return only a local segment that, however, has the maximum value in terms of alignment score.
To get a global alignment, another algorithm has to be used. Needleman-Wunsch would be an example.
The difference between local and global is that a local alignment only alignes a subsequence of the provided sequences whereas a global alignment gives the best alignment of the whole sequences.


text:
PubMed Central Open Access archives are corrupted
I'm trying to download PubMed Central Open Access articles from ftp://ftp.ncbi.nlm.nih.gov/pub/pmc/oa_bulk. However, I am not able to decompress the archives without errors.
tar always gives me the error:
gzip: data stream error
tar: Unexpected EOF in archive
tar: Error is not recoverable: exiting now

Does PubMed provide corrupted archives or does downloading corrupt the archive?

Have you tried to decompress it in the past few days since it's been updated? Picking a smaller file, tar -zxvf non_comm_use.A-B.xml.tar.gz works fine.

I was eventually able to do it on Linux. It didn't work on MacOS for me (not even with gnu-tar).

The download attempt might be caused by an attempt to download the files while these files are being updated by the NCBI server. Try downloading on other days and make sure you have a strong connectivity.


text:
Detecting SARS-CoV-2 variants from the mixed virus population
I have a fastq file from Nanopore sequencing data that contains reads from both the UK and South Africa variants of SARS-CoV-2.
The variants are identified by three key mutations in the receptor-binding domain of Spike.
Aim
Firstly, I want to determine the number of reads then belongs to UK and South Africa variant. Then, I will create a bam file from the fastq file using the SARS-CoV-2 reference genome.
Questions

How do I determine the reads belonging to each variants?
Are there any tools available to classify the reads based on the pre-defined mutation criteria?

I don't think that's possible in general. You can obviously map the reads to the references and split but that's if you trust the references and in that case why do you even need the reads? Just thinking out loud.

Based on a quick request from someone else, I wrote some code a couple of years ago for assigning read groups based on variants at particular locations; this was set up for HIV variant splitting, so I expect it should also work on SARS-CoV-2 nanopore reads as well:
https://gitlab.com/gringer/bioinfscripts/-/blob/master/samVarSplitter.pl
It processes a SAM stream [i.e. in text mode] using Perl, so might be a bit slow, but should be fine for a few hundred thousand sequences. You'd use it as follows:
samtools view -h mapped_reads.bam | \
  samVarSplitter.pl [-ref <refName>] [-pos <int>] | \
  samtools sort > grouped_reads.bam

Then the reads can be split into different files based on their group using samtools split:
Usage: samtools split [-u <unaccounted.bam>] [-h <unaccounted_header.sam>]
                      [-f <format_string>] [-v] <merged.bam>

Fantastic, good response.


text:
Bioconductor Conda release delay
How long does it take for the latest Bioconductor to end up on Conda? I think I saw a blog post about the challenges with this process, but I cannot find it. For example, Bioconductor 3.13 was released a week ago, but is still not available. For comparison, many CRAN packages tend to show up within a day.
Additionally, is there a way to keep track of that progress?

I started the roll-out on the 30th in the evening. There's not much of a delay for the October Bioconductor release, but for the summer release we first have to build the newest R version and then rebuild all conda CRAN packages. It's probably not widely known, but Bioconductor packages are only compatible with a single R version. For example, the ones being released now are not compatible with R 4.0 and will not be compatible with R 4.2, but will instead only ever work with R 4.1. This is regardless of what the metadata provided by Bioconductor package authors might say, the Bioconductor core team has requested this (there's actually not even a guarantee that packages built under R 4.1.0 will work with 4.1.1, but such issues are exceedingly rare). This is why there is a multi-week delay every summer.
On the plus side, at least we're not waiting for the first bug-fix release of R to start the migration anymore.
For Bioconda we've never gotten around to making a migration tracker. If you have a local clone of Bioconda, you can run ./scripts/bioconductor/rootNodes.py --sinceNDays 2 --rootNodes to get a count of the packages completed in the past 2 days and a list of the root nodes on the dependency DAG that still need to be handled. That's how I estimate the fraction finished (currently about 600 packages).
For tracking the status on conda-forge you can see their status page: https://conda-forge.org/status/#r410

Devon will probably add a more precise answer, until then I just link his Tweet from last week:
https://twitter.com/dpryan79/status/1394714988720308226
Edit: 31.5.21
Rollout started yesterday, should take a week to complete:
https://twitter.com/dpryan79/status/1399273322601459712

That's promising. Would be great if there is was a way to track progress without waiting for the next tweet.


text:
Long-form gene summary
The following
(thx to these answers)
retrieves the info, in particular a long summary (such as "The membrane-associated protein encoded by this gene is a member of the superfamily of ATP-binding cassette... etc."), for a few mouse genes.
library(dplyr)

org.Mm.eg.db::org.Mm.egSYMBOL2EG %>%
  as.data.frame() %>%
  pull(gene_id) %>%
  head() %>%  # !!!
  rentrez::entrez_summary(id = ., db = "gene") %>%
  { do.call("rbind", .) } %>%
  as.data.frame() %>%
  select(name, description, summary)

Note the head. The protocol/server does not allow large queries. How can I get the same information (incl. the long summary) for all the mouse genes at once?

Chunk your query, then use purrr::insistently. No need for a loop, this can be vectorised. I also recommend against a fixed-length sleep; using an exponential backoff is more polite.

Thanks, I'll certainly have occasion for this one.

You could wrap the call too entrez_summary in a loop.
For example to cut the data into chunks of five genes:
gene.vector # this would be the genes you are interested in
gene.chunks <- ceiling(length(gene.vector) / 5)
# initializing an empty list
summary.list <- vector(mode = "list", length = gene.chunks)
for(i in gene.chunks){
   # select the ith set of five genes from your gene vector
   goi <- gene.vector[(i - 1) * 5 + 1:5]
   # if there are less than five elements there will be NAs that need to be removed
   goi <- goi[!is.na(goi)]
   # adding your code
   summary.list[[i]] <- org.Mm.eg.db::org.Mm.egSYMBOL2EG %>%
     as.data.frame() %>%
     filter(gene_id %in% goi) %>% # here the current five genes are selected
     pull(gene_id) %>%
     rentrez::entrez_summary(id = ., db = "gene") %>%
     { do.call("rbind", .) } %>%
     as.data.frame() %>%
     select(name, description, summary)

   Sys.sleep(10) # wait for 10 seconds
}

There will also be a way to do this using dplyr group. In this case the loop is a little easier to debug.
I added Sys.sleep at the end to reduce the frequency with which you are making requests from entrez. If you are posting too many requests in a short time this will be like a DDOS attack. You (and potentially your whole IP range including other people at your institution) could get banned from the service for a couple of hours. Maybe try it with a subset first.
If you really want to get the information for all mouse genes this probably will be rather slow. Alternatively, you can go to the entrez webpage and try to download the information as a list. Unfortunately, I don't know exactly where you will find this information.

The following bash script downloads the summaries into an XML file.
#!/usr/bin/env bash

# This requires the ncbi tool `gene2xml`
# sudo apt install ncbi-tools-bin

DIR="download_cache"
mkdir -p "$DIR"

URL="ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/ASN_BINARY/Mammalia/Mus_musculus.ags.gz"
wget -S -O - "$URL" | gzip -d | gene2xml -b T -l T | gzip -9 > "$DIR/mm.xml.gz"

URL="ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/ASN_BINARY/Mammalia/Homo_sapiens.ags.gz"
wget -S -O - "$URL" | gzip -d | gene2xml -b T -l T | gzip -9 > "$DIR/hs.xml.gz"

These XML files can be converted to tables along these lines, viz.:

import xml.etree.ElementTree as ET
import gzip

from pathlib import Path

from tcga.utils import mkdir, unlist1, relpath

def parse(src):
    def find_text(e, path):
        e = e.findall(path)
        if e:
            return unlist1(e).text
        else:
            return ""

    yield ('gene_id', 'description', 'summary')

    with gzip.open(src, mode='r') as fd:
        e: ET.Element
        for (x, e) in ET.iterparse(fd, events=['end']):
            if e.tag == "Entrezgene":
                gene_id = find_text(e, "**/Gene-track_geneid")
                descptn = find_text(e, "**/Gene-ref_desc")
                summary = find_text(e, "Entrezgene_summary")

                yield (gene_id, descptn, summary)

                e.clear()

def main():
    for src in Path(__file__).parent.glob("download_cache/*.xml.gz"):
        trg = mkdir(Path(__file__).with_suffix('')) / Path(src.stem).with_suffix('.tsv.gz')
        if trg.is_file():
            print("Skipping existing file:", relpath(trg))
        else:
            print("Writing:", relpath(trg))
            with gzip.open(trg, mode='wt') as fd:
                for t in parse(src):
                    print(*t, sep='\t', file=fd)

if __name__ == '__main__':
    main()

Btw., a table of vertebrate gene homology is available here (it's just a tab-separated table), and can be used to impute summaries by homology.


text:
Remove repetitive region vcf file using repeatmasker bed file
I have a 1000 genomes vcf file for chromosome 14. And I want to remove variants in repetitive regions flagged by repeatMasker. I have a bed file for those repeats. I downloaded that file from UCSC. How can I eliminate/remove those sites in vcf file using repeat masker's bed file? Insights will be appreciated.

A couple fast alternatives:
$ bedops --chrom chr14 -n 1 <(vcf2bed < variants.vcf) <(sort-bed repeats.bed) > answer.bed
$ bedops --chrom chr14 -n 1 <(vcf2bed < variants.vcf) <(rmsk2bed < repeats.rmsk) > answer.bed

Using --chrom chr14 focuses set operations on that chromosome.

This worked for me. Thanks Alex Reynolds

vcftools with the --exclude-bed option or bedtools with the subtract command should do the trick.

how about using bedtools intersect -v?

probably better to use bcftools view -R .bedfile - vcftools isn't maintained any more and a fair bit slower than bcftools.


text:
where to find sample contig data
Where or How can I find contigs? I am trying to learn Bioinformatics and I want to make a reference-based gene search. I also want to align contigs with a given reference sequence. From NCBI other databases I can find only reads and complete genome sequences. Please tell me how can I find contig sequences of a species?

Most of the reference genomes on NCBI are made of contigs or scaffolds (which is just a fancy name for contigs glued together by a bunch of unknown nucleotides). Hence, by download the whole genome, you also download all the contigs.
If you search for a specific contig (presumably a known locus) of a specific species, you can search NCBI for the species name + locus / scaffold / contig name and if it was sequenced it will pop up in the "genes" section of NCBI. For example, here is a search for Pan troglodites contig NC_036898.1 (I typed literally Pan troglodites NC_036898.1 in the search bar). You can even just type a species name, and click on "genes" and see loads of all the individual contigs sequenced.
--- edit ---
There is a bunch of ways to generate something we could call contig. However, usually what we mean is "continuous sequence assembled from sequencing reads" and it's generated from whole-genome sequencing reads. We sequence a genome, assemble the reads and what we get is a set of contigs as a representation of a haploid reference genome. Some of them will be short fragments of nucleotides with unknown functions, some of them can be long and contain hundreds of genes. The important part is that it's a continuous sequence.
With a bit more relaxed definition of "contig", we could also accept as contigs outputs of targeted sequencing (including old-fashioned sanger). It's a continuous stretch of nucleotides, usually bordered by a pair of primers. I am not 100% sure, but my best guess that most of the "gene" sequences in NCBI will be sanger-based (but I should disclaim I never investigated too much where all the NCBI sequences come from).

This helps thank you but how can i find the contigs/scaffolds before they were assembled to a reference sequence/whole genome. In the genes section of NCBI, I can find list of annotated genes with their sequence, i thought contigs/scaffolds are longer than one gene and contain other features?

I tried to clarify the answer to address your questions. Is it better now?


text:
Masking sites in a vcf file
I need to mask all sites in a vcf file flagged by the 1000 Genomes Project as being unfit for population genetic analyses. The sites for all chromosomes are available at:
1000Genomes masked sites
From above mentioned link you can get both bed files and Fasta files.
I wanted to mask all sites using VCFtools but VCFtools is not longer being updated and therefore, can't me used on the my vcf file because you can only VCFtools upto VCF file format v.2. And file format for my file is v.3. Can someone suggest another tool that I can use. Insights will be appreciated.

I've tried many ways to do this with bcftools but don't think it's possible to exclude a bed file. What you want is
bedtools intersect -v -a sample.vcf -b mask.bed -wa -header
The -v option gives you only the parts of file a that are not in file b. -wa only keeps entries from file a, and -header preserves the header from file a.

If you are familiar with Python, check out the pyvcf.VcfFrame.filter_bed() method I wrote:
from fuc import pyvcf
vf = pyvcf.VcfFrame.from_file('your_vcf.vcf')
filtered_vf = vf.filter_bed('your_bed.bed', opposite=True)
filtered_vf.to_file('filtered.vcf')

Have a look at BCFtools, which I think does most of the stuff that VCFtools does. I notice there's a --mask option in BCFtools consensus, which may do what you want:
-m, --mask FILE
 BED file or TAB file with regions to be replaced with N
 (the default) or as specified by the next --mask-with
 option. See discussion of --regions-file in Common
 Options for file format details

So just to make sure that I am using the correct command. It'll be:  bcftools consensus --mask A.fa A.vcf.gz > masked_file.vcf. Is this correct?

@John Sorry, I don't know anything more than what I've read; I have no experience with using masking in BCFtools or VCFtools. Try it and see if it works!


text:
Can ortholog genes have different names among different species?
To be specific, we have different lineages of E. coli. I wondered if the ortholog genes among them should all have the same name in all lineages or can they also have different names?

Gene naming is pretty subjective and even a single gene can have many names. If you're lucky orthologs will be named the same in different lineages, but you can't rely on it. Better to use gene ids (e.g. from Ensembl) and an ortholog db to map the orthologs. However, if you're looking at different strains, this is more tricky since most db don't go below the species level

Thanks for your answer. I was wondering if two genes belong to two different isolates of E. coli (i.e., they are still the same species but different isolates), should we call them ortholog or paralog

Gene IDs are not perfect either. Those can change as well, even between different versions of the annotation for the same species.


text:
How do you convert Raw Alignment Score to Bit Score?
I'm coding a pipeline where I make a lot of pairwise alignments, and I end up with raw alignment scores. But, I really need to look at my results in terms of bit scores.
I know that the formula is:
$′ = (ƛ* − ()) / (2)$
But, I don't know what ƛ or K would be. I'm currently using BLOSUM62.
I've looked around online for how to look up (and even empirically deduce) lambda and k, and I all I've come across were convoluted explanations.
I also know that this question from almost 3 years ago attempts to address the same topic, but the best response only leads to an NCBI link that requires login credentials.
Any help would be much appreciated.

Answer from @terdon, converted from comment:
The link just points to the NCBI ftp site where you can download the source package and look at the code. There is no need for login credentials, you can log in as anonymous. It's just that the version has changed so that answer's URL no longer works. Try this instead:
wget https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.11.0+-src.tar.gz

Then, extract the files and look at the file c++/src/algo/blast/core/blast_stat.c which, again according to the answer you linked to, should have the lamda and k tables at line 589.


text:
using Bowtie to map miRNA-seq data to rfam and also reference genome
I am using Bowtie to remove non-coding RNAs (tRNA, snRNA, rRNA) by rfam and also maping our microRNA-seq data with mirbase using following code:
bowtie -m 10  -l 18 -n 1 -v 2 -e 80 --best --strata -a   genome.index.file    inputfastqfile   --un unaligned.fq

Unfortunately I dont have any references or index for that code.
Are there any other settings and options to improve this code?

Why rfam specifically? You can you rRNA, tRNA, and snRNA database to get their sequences (assuming human genome) and use a hierarchical sequence alignment strategy to remove all the ncRNAs and, finally, map the filtered sequences to miRNAs.


text:
Entrez (Biopython) esearch and efetch not returning sequence as expected
I'm trying to use Entrez (through Biopython) to download the sequence of a TMV replicase gene.
I have the following code:
# Lookup ID
search = Entrez.esearch(db='gene', term='Tobacco mosaic virus[Orgn] AND replicase')
read = Entrez.read(search)
idlist = read["IdList"]

# Fetch sequence
search = Entrez.efetch(db='gene', id=idlist[0], retmode='text', rettype='gb')
read = SeqIO.read(search, "genbank")
tmv_replicase = read.seq

But when I run it it throws the exception on the SeqIO.read line: ValueError: No records found in handle.
The IDs provided in idlist are: ['1494081', '1494082'].

Did you try to fetch from the database nucleotide instead of gene ?
# load modules
from Bio import SeqIO
from Bio import Entrez

# Lookup ID
search = Entrez.esearch(db='nucleotide', term='Tobacco mosaic virus[Orgn] AND replicase')
read = Entrez.read(search)
idlist = read["IdList"]

# Fetch sequence
search = Entrez.efetch(db='nucleotide', id=idlist[0], retmode='text', rettype='gb')
read = SeqIO.read(search, "genbank")
tmv_replicase = read.seq

It returns
>>> read
SeqRecord(seq=Seq('GTATTTTTACAACAATTACCAACAACAACAAACAACAAACAACATTACAATTAC...GGG'), id='MN912489.1', name='MN912489', description='Tobacco mosaic virus isolate TMV:BR:SC:02:01, complete genome', dbxrefs=[])
>>> tmv_replicase
Seq('GTATTTTTACAACAATTACCAACAACAACAAACAACAAACAACATTACAATTAC...GGG')

Thank you, this gets me what I wanted. Out of interest though, do you know why using the gene database for this wouldn't work? (I'm new to Entrez).

It is simply 2 differents databases. The ID into nucleotide is not the same into gene. Thus from the website NCBI it is possible to extract a fasta from a gene query it doesn't seem possible with biopython https://biopython.readthedocs.io/en/latest/chapter_entrez.html

If you want to filter your queries from the nucleotide database in order to get only the specific genes you can check my code as an example: https://gitlab.mbb.univ-montp2.fr/bev/target_genes_from_ncbi/-/blob/master/select_query_gene_sequences_ncbi.py


text:
rmsd between re-docked complex and co-crystallized complex
I would like to calculate the rmsd between re-docked complex and co-crystallized complex for docking validation. How do I calculate the rmsd using Pymol? I tried the following command
align re-docked_complex, co-crystallized_complex

But the above command gives the rmsd of two proteins. I would like to get the rmsd of re-docked and co-crystallized ligand. How can I achieve this using pymol?

Welcome.
The Protein-Ligand-complex is the complete protein with Ligand!
If you only want to compare the two ligands, you can delete the rest of the proteins, or specify the input so only the ligand is choosen. (select Ligand and save selection. then use the function on that selections)
But be aware that for Docking comparison you probably don't want to align! but rather calculate the RMSD directly (because the absolute position and not the relative position is of interest?)
Maybe use: "rms_cur" (https://pymolwiki.org/index.php/Rms_cur)

Yup, spot on. rsm_cur is the correct command. One can pass the ligand w/o deleting: rsm_cur redocked and resn LIG and not element H, original and resn LIG if the residue name is LIG say, while the not element H is against hydrogens that were added for docking. However, one big caveat worth mentioning is...

... is that the RMSD is segi, chain, resi, resn and name based (i.e.alter/sort operations may be required to make them match), has an isomorphism problem: Take BHT, if you flip it's "arms" 180° around the hydroxy–methyl axis, you end up with a perfect fit, but the atom names differ. So, one must make sure this is not an issue (or use RDKit).


text:
What are the units for the force output by nstfout, the velocity output by nstvout, and the coordinates output by nstxout in gromacs?
What are the units of force, velocity, and coordinates output by gromacs trajectory? In .mdp, they are specified by nstfout,nstvout,nstxout.

nstxout, nstvout and nstfout are the number of steps between coordinates, speed and forces written in the output trajectory (Output control documentation).
As for the units present in your trajectory files (say a .trr):
coordinates are in nm, speeds in nm.μs⁻¹ (so in mm.m⁻¹) and forces in kJ.mol⁻¹.nm⁻¹. These should be consistent through all GROMACS files (Definitions and Units documentation).
Note that .xtc files multiply coordinates by 1000 to compress them, so they are stored in pm (see GROMACS xtc documentation). Given MDAnalysis documentation on the subject, I believe .xtc files can only store coordinates, no speeds or forces.
I would be wary of softwares used to analyze your trajectory results (MDAnalysis, VMD, etc), as they may automatically do the conversion of these units, and for example show you coordinates in Å.

This was a great answer O.Laprevote! We have a few unanswered GROMACS questions here, are you able to help us out at all with those? I'm trying to help take care of the unanswered queue.


text:
How to make CIRCOS plot of VCF file?
I have 3 VCF files, I want to make CIRCOS plot just like the following:
https://www.nature.com/articles/s41598-020-70527-8/figures/4
I try to follow the tutorial but could not understand. The genome is of Zea mays and the reference files used in making VCF files are downloaded from Ensembl plant. How I can make karyotype and config files?


text:
How to input the opposite strand of DNA?
I am trying to write a program that will count the repeats of a given DNA strand:
I am using a Counter library, something like this code shown below
from collections import Counter

        s="XXXXXAGCCTGCCAAGCAAACTTCACTGGAGXXXXXTGTGCGTAGCATGCTAGTAACTGCAXXXXXTCTGAATCTTTCAGCTGCTTGTXXXXXTGGGCCTCTCACAAGGCAGAGTGTCTTCATGGGACTTTGATATTTATTTTTGTACAACCTAAGAGGAACAAATCCTTTGACACTGACAAATTGGCTTCCATATTTTATACCTTAATCATCTCCATGTTGAATTCATTGATCAACAGTTTAAGAAAAAAAGATGTAAAAATGCTTTTAGAAAGAGAGGCAAAGTTATGCACAATAACTTCTCATGAAGTCACAGTTTGTTAAAAGTTGCCTTAGTTCACAATAAATAATTATGTATGCTCTATAATTTCAGTGA"

def chunk_string(s, n):
    return [s[i:i+n] for i in range(len(s)-n+1)]

counter = Counter(chunk_string(s, 5))

Repeat supposed to be "declared" as an identical sequence on the same strand (like XXXXX here that occurs 4 times). However, I can also have this repeat (for ex. XXXXX in this particular test case) on the OPPOSITE strand.
How am I suppose to intake the opposite strand in case of using .fasta files? As I know for today - fasta is presented only in a single strand way. Or am I completely mixed the terms?

Yes, you're going to want to take the reverse complement of your string in order to analyze the other strand, that is the standard.
Also, it is worth pointing out that you may want to use a generator expression instead of a list expression for your chunk_string function:
def chunk_string(s, n):
    return (s[i:i+n] for i in range(len(s)-n+1))

This will result in a generator that yields kmers as the Counter constructor iterates through them, which means that you don't need to first construct and store the entire list of kmers first and then pass to Counter but instead can feed them directly to it. Not essential for small strings, of course, but a good practice and useful for larger strings. See this SO post for a better explanation.
You can also use Biopython's reverse_complement function in order to easily obtain the reverse complement.

What's going to be simplest is to generate the reverse complement of your short search sequence, and search for that too.

Hi @swbarnes2, this looks like the start of a good answer; any chance you could expand on this approach?


text:
Using Entrez.efetch() to retrive .fasta file from any NCBI database?
I tried to use these python codes to retrieve a .fasta file from gene database, but the output were the parameters of the report
from Bio import Entrez
Entrez.email = "send_helps@earth.com"
handle = Entrez.efetch(db = "gene",id = "682",rettype="fasta", retmode="xml")
print(Entrez.read(handle))

i also read the documentation on this site about efetch() and seems like this function can only retrieve certain kinds of data for each database. So i can't take fasta file from gene database !?
I don't want to use the NCBI web UI. What's wrong with my code? Is there a way i can guarantee to retreive a .fasta file through python from NCBI ? Am i misundertanding something severly ?
Sorry for bad English. thank you.

In order to download the FASTA sequence of transcript variants using Entrez eutilities, you will have to use Entrez.elink to get a list of transcript accessions for the gene of interest and then use Entrez.efetch to fetch the sequence(s). Just to give you an idea, you can use Entrez Direct for this as follows:
elink -db gene -id 682 -target nuccore -name gene_nuccore_refseqrna \
  | efetch -format fasta

You can find more Entrez Link descriptions here. Note, there is a gene_fasta option for efetch but it still returns only the transcript sequence (i.e., no intron sequence).
If you need the full gene sequence (including the intronic regions), you can use the commands mentioned here.
Alternatively, you can use NCBI Datasets for this. You can use the web interface, command line tool or the python API. This page with the how-to articles has the information you seek.


text:
Tajima-Nei Distance estimate with BioPerl
I have been trying to estimate Tajima-Nei distance for my data (link below). The output is s phylogenetics distance (pairwise) matrix which accommodates a detailed transition transversion and the Jukes-Cantor correction.
I´m following this protocol from BioPerl: https://metacpan.org/pod/Bio::Align::DNAStatistics
I have 314 sequences in a fasta file and another file with the list of IDs.
Fasta:
>AVP78031.1
----atgttgtttttcttgtttcttcagttcgccttagtaaactc---------------
------------------------------------ccagtgtgttaacttgacaggcag
a----------------accccactcaatcccaattat--actaattcttcacaaagagg
...

IDs:
AVP78031.1
...

I´m using this Perl script to calculate the Tajima-Nei Distance in a pairwise comparison (314 * 314):
use strict;
use warnings;
use Bio::AlignIO;
use Bio::Align::DNAStatistics;

my $file = $ARGV[0];
my $idfile = $ARGV[1];

if ($file eq "" ) {
  $file = "NT_MSA_S_protein.fasta";
} elsif ($idfile eq "" ) {
  $idfile = "NT_ID_S_protein.csv";
}

#### Considerando un archivo
my @contentIDS;

open (LIST, $idfile) or die;
while (my $l = <LIST>) {
  $l =~ s/\n//g; # delete newline
    $l =~ s/\r//g; # delete CR
  next if (length($l) < 1);
  push @contentIDS, $l;
}
close LIST;

#### .... IDs list
my $stats = Bio::Align::DNAStatistics->new();
my $alignin = Bio::AlignIO->new(-format => 'fasta', -file   => $file);  ### $file: MSA file
while (my $aln = $alignin->next_aln) {
  #print "reading...A\n"; ### DIAG
  my $matrix = $stats->distance(-align => $aln, -method => 'Tajima-Nei');
  #print "reading...B\n"; ### DIAG
  ### Obtaining values for each pair (DISTANCE!)
  WL1:
  foreach my $aaa (@contentIDS) { ### ID #1
    WL2:
    foreach my $baa (@contentIDS) { ### ID #2
    next (WL2) if ($aaa eq $baa);
  my $data =  $matrix->get_entry($aaa, $baa);
  #($data = 0) if ($data < 0);
    print "DISTANCE\t$aaa\t$baa\t$data\n";
    } # END WL2
  } # END WL1
}

exit;

#

This script work it when I tried with small data, however, when I try with my real data this is the error message
MSG: Must provide a DNA alignment to Bio::Align::DNAStatistics, you provided a protein
---------------------------------------------------
Can't locate object method "get_entry" via package "0" (perhaps you forgot to load "0"?) at Tajima-Nei_Distance_NV.pl lin$

This is weird because I review my data for ambiguous characters and the characters are in majority "atcg" and on some occasions "n", at least that there are other ambiguous characters (maybe) that represent protein sequence.
I really don´t understand the message because the fasta file is clearly a nucleotide sequence.
Anyone with experience using Bioperl and estimation Tajima-Nei distance?
Link: https://github.com/MauriAndresMU1313/Example_Tajima-Nei_Distance_Bioperl/tree/main
Any comment or help is welcome!!
Thank!

The problem is the code itself.
I suspect there are R (argenine) and Y (tyrosine) in the sequence data. These are common ambiguities but also represent amino acids. The distance estimate has not accommodated these nucleotide ambiguities but is doing a getter setter for amino acids.
It would be easier to perform the same calculation in Mega X than debugging the code.


text:
GROMACS Coordinates and velocity and force trajectories do not match accuracy
I am using GROMACS to perform MD simulations of proteins. So, I tried to record the simulation results nstxout, nstvout, and nstfout with the same step size (nstxout 1, nstvout 1, nstfout 1). Then, I calculated the coordinate (step n) of an arbitrary step of the trajectory (step n) from the velocity of step n - 1, step n, and calculated the velocity of n - 1/2, and used that velocity to calculate the coordinate (n-1) of step n using the velocity of n - 1/2. LeapFrog was used in the simulation.
How can I get the coordinates, velocity, and force to be consistent?

From GROMACS documentation about reproducibility:

The following factors affect the reproducibility of a simulation, and
thus its output:

Precision (mixed / double) with double giving “better” reproducibility.
Number of cores, due to different order in which forces are accumulated. For instance (a+b)+c is not necessarily binary identical
to a+(b+c) in floating-point arithmetic.
Type of processors. Even within the same processor family there can be slight differences.
Optimization level when compiling.
Optimizations at run time: e.g. the FFTW library that is typically used for fast Fourier transforms determines at startup which version
of their algorithms is fastest, and uses that for the remainder of the
calculations. Since the speed estimate is not deterministic, the
results may vary from run to run.
Dynamic linking to different versions of shared libraries (e.g. for FFTs)
Dynamic load balancing, since particles are redistributed to processors based on elapsed wallclock time, which will lead to (a+b)+c
!= a+(b+c) issues as above
Number of PME-only ranks (for parallel PME simulations)
MPI reductions typically do not guarantee the order of the operations, and so the absence of associativity for floating-point
arithmetic means the result of a reduction depends on the order
actually chosen
On GPUs, the reduction of e.g. non-bonded forces has a non-deterministic summation order, so any fast implementation is
non-reprodudible by design.

I would add that in your case of hard coding everything by hand you may have to look a bit into your barostat and thermostat (although their influence should be found on speed of the atoms). Maybe a sample code with coordinates and speed of 4 atoms would help to give a more precise answer.
You can also try having a look at gromacs source code. I believe OpenMM provides a bit more control and clarity on each step, if for some unfathomable reason you are not doing this only for curiosity sake but because you actually have to make a code that reproduces these results.


text:
How to increase the size of the shapes that denotes effect sizes of different datasets and metaanalysis?
I am doing metaanalysis with forestplot from ggforestplot package. I have 7 different datasets from 7 different countries. After metaanalysis, I have 7 Effect sizes and a Meta-analysis effect size, i.e., total of 8 effect sizes. Now, I have generated a forest plot using this code:
# Forestplot
my_plot <- forestplot(
  df = final,
  estimate = estimate,
  logodds = FALSE,
  colour = country=="metaanalysis",
  shape = country,
  title = "Associations to disease accross populations",
  xlab = "Effect Size",
  #xlim = c(-1.5, 1.5),
  xtickbreaks = c(-1.0, -0.5, -0.45, -0.4, 0, 0.4, 0.45,0.5, 1.0)
)+
  # You may also want to add a manual shape scale to mark meta-analysis with a
  # diamond shape
  ggplot2::scale_shape_manual(
    values = c(5L, 0L, 1L, 2L, 8L, 19L, 4L, 3L),

  )

But, the shapes denoting the countries and the metanalysis seem very small. How can I make them larger for better visualization? Can anyone please help me?
Representative data for replication:
final <- data.frame(
  stringsAsFactors = FALSE,
              name = c("RIBOSYN2-PWY__flavin_biosynthesis_I_(bacteria_and_plants)",
                       "PWY-6385__peptidoglycan_biosynthesis_III_(mycobacteria)",
                       "PWY-6151__S-adenosyl-L-methionine_cycle_I",
                       "PWY-6387__UDP-N-acetylmuramoyl-pentapeptide_biosynthesis_I_(meso-diaminopimelate_containing)",
                       "PWY-6700__queuosine_biosynthesis",
                       "PEPTIDOGLYCANSYN-PWY__peptidoglycan_biosynthesis_I_(meso-diaminopimelate_containing)",
                       "PWY-6386__UDP-N-acetylmuramoyl-pentapeptide_biosynthesis_II_(lysine-containing)",
                       "PWY-5667__CDP-diacylglycerol_biosynthesis_I",
                       "PWY0-1319__CDP-diacylglycerol_biosynthesis_II",
                       "PWY-6122__5-aminoimidazole_ribonucleotide_biosynthesis_II",
                       "PWY-6277__superpathway_of_5-aminoimidazole_ribonucleotide_biosynthesis",
                       "PWY-6163__chorismate_biosynthesis_from_3-dehydroquinate",
                       "TRPSYN-PWY__L-tryptophan_biosynthesis",
                       "PANTOSYN-PWY__pantothenate_and_coenzyme_A_biosynthesis_I",
                       "GLYCOCAT-PWY__glycogen_degradation_I_(bacterial)",
                       "P161-PWY__acetylene_degradation",
                       "FUC-RHAMCAT-PWY__superpathway_of_fucose_and_rhamnose_degradation",
                       "P108-PWY__pyruvate_fermentation_to_propanoate_I",
                       "PWY-7328__superpathway_of_UDP-glucose-derived_O-antigen_building_blocks_biosynthesis",
                       "TCA__TCA_cycle_I_(prokaryotic)",
                       "PWY-6901__superpathway_of_glucose_and_xylose_degradation","PWY66-400__glycolysis_VI_(metazoan)",
                       "FUCCAT-PWY__fucose_degradation",
                       "PRPP-PWY__superpathway_of_histidine,_purine,_and_pyrimidine_biosynthesis",
                       "PWY-7234__inosine-5'-phosphate_biosynthesis_III",
                       "TRPSYN-PWY__L-tryptophan_biosynthesis",
                       "PWY-6122__5-aminoimidazole_ribonucleotide_biosynthesis_II",
                       "PWY-6277__superpathway_of_5-aminoimidazole_ribonucleotide_biosynthesis",
                       "RIBOSYN2-PWY__flavin_biosynthesis_I_(bacteria_and_plants)",
                       "PWY-6163__chorismate_biosynthesis_from_3-dehydroquinate",
                       "PWY-5667__CDP-diacylglycerol_biosynthesis_I",
                       "PWY0-1319__CDP-diacylglycerol_biosynthesis_II",
                       "PANTOSYN-PWY__pantothenate_and_coenzyme_A_biosynthesis_I",
                       "RIBOSYN2-PWY__flavin_biosynthesis_I_(bacteria_and_plants)",
                       "PWY-6386__UDP-N-acetylmuramoyl-pentapeptide_biosynthesis_II_(lysine-containing)",
                       "RIBOSYN2-PWY__flavin_biosynthesis_I_(bacteria_and_plants)",
                       "PANTOSYN-PWY__pantothenate_and_coenzyme_A_biosynthesis_I",
                       "PWY-6385__peptidoglycan_biosynthesis_III_(mycobacteria)",
                       "PWY-6387__UDP-N-acetylmuramoyl-pentapeptide_biosynthesis_I_(meso-diaminopimelate_containing)","PWY-6700__queuosine_biosynthesis",
                       "PWY-6151__S-adenosyl-L-methionine_cycle_I",
                       "PWY0-1319__CDP-diacylglycerol_biosynthesis_II",
                       "PWY-5667__CDP-diacylglycerol_biosynthesis_I",
                       "PEPTIDOGLYCANSYN-PWY__peptidoglycan_biosynthesis_I_(meso-diaminopimelate_containing)",
                       "PWY-6387__UDP-N-acetylmuramoyl-pentapeptide_biosynthesis_I_(meso-diaminopimelate_containing)",
                       "PWY-6385__peptidoglycan_biosynthesis_III_(mycobacteria)",
                       "PANTOSYN-PWY__pantothenate_and_coenzyme_A_biosynthesis_I",
                       "PEPTIDOGLYCANSYN-PWY__peptidoglycan_biosynthesis_I_(meso-diaminopimelate_containing)","PWY-5667__CDP-diacylglycerol_biosynthesis_I",
                       "PWY0-1319__CDP-diacylglycerol_biosynthesis_II"),
          estimate = c(0.492903365186722,
                       0.490116714405755,0.470681520902659,0.465726944992238,
                       0.461786439528495,0.455469435286562,0.449174326308655,
                       0.433755737833203,0.433670791448977,0.412009484469875,
                       0.412009484469875,0.409849212159568,0.407818381473059,
                       0.406187518931948,-0.405573831841341,-0.407806403622317,
                       -0.418351108232386,-0.419107645158175,-0.420781760133428,
                       -0.441930843508614,-0.444723954813865,-0.451224756237618,
                       -0.504642907784859,-0.515843696383127,-0.523823484265933,
                       1.01203740887784,0.860911808112254,0.860911808112254,
                       0.782824395343453,0.779871038531055,0.763412247455786,
                       0.763412247455786,0.744965434413444,0.735086809796288,
                       0.733132943355509,0.731855318166865,0.730523378560539,
                       0.728356885227315,0.694487466356475,0.694335496605963,
                       0.68469489116557,0.680753784497535,0.680732580641012,
                       0.668909750729078,0.665107737526883,0.662352062210176,
                       0.643767077562848,0.640204416487105,0.639460661647805,
                       0.637911093993277),
                se = c(0.10858902830706,
                       0.101375498952004,0.101193730594653,0.101213265518358,
                       0.101165178275411,0.101149604926044,0.101114528367103,
                       0.101177617463096,0.101176196339295,0.10102421433221,
                       0.10102421433221,0.100910703206783,0.101046445558438,
                       0.12265509469116,0.10486856029419,0.103267024648877,0.100867566762125,
                       0.11467836793049,0.10790947678297,0.100997454994136,
                       0.101684644439116,0.103275688295495,0.101358559584742,
                       0.101485340026678,0.131050876787633,0.12672527384678,
                       0.0911987727585976,0.0911987727585976,0.0473284627336248,
                       0.120947206065826,0.120594420271743,0.120594420271743,
                       0.120207965256539,0.12000489747134,0.119965054342127,
                       0.0890575546896603,0.0468935746714197,
                       0.0890043416222531,0.119198789457309,0.0466102635462743,
                       0.119011209638725,0.0883056120686039,0.0883053113523379,
                       0.0881390968551948,0.0880862640207957,0.118593197976584,
                       0.0877953030586674,0.118192523540133,0.0462081105569776,
                       0.0461972348759043),
           country = c("metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","Swedish","Kazakh","Kazakh",
                       "Japanese","Swedish","Swedish","Swedish","Swedish",
                       "Swedish","Swedish","Kazakh","Japanese","Kazakh",
                       "Swedish","Japanese","Swedish","Kazakh","Kazakh",
                       "Kazakh","Kazakh","Swedish","Kazakh","Swedish","Japanese",
                       "Japanese")
)

Output:

Thanks

A quick look at the code for forestplot shows that geom_effect is what needs to be modified. If you add a value for fatten changes the size of the shapes in the the figure (but not in the legend), or a value for size will change the shape, line and legend sizes (add to not geom_effect not to the aes layer)
This may be easier said than done as forestplot also seems to call an unexposed function (function quo_is_null). A couple of simple hacks are possible; copy and make your own forestplot in the default environment and remove the call to the function (it seems to just check if pvalue is set), or change the environment of the "new" forestplot to ggforestplot
envorionment(forestplot) <- envornment(ggforestplot::forestplot)

There may be other less awkward methods, but this will work. Fatten changed to size 6 which is massive.
forestplot <-
function (df, name = name, estimate = estimate, se = se, pvalue = NULL,
          colour = NULL, shape = NULL, logodds = FALSE, psignif = 0.05,
          ci = 0.95, ...)
{
  stopifnot(is.data.frame(df))
  stopifnot(is.logical(logodds))
  name <- enquo(name)
  estimate <- enquo(estimate)
  se <- enquo(se)
  pvalue <- enquo(pvalue)
  colour <- enquo(colour)
  shape <- enquo(shape)
  args <- list(...)
  const <- stats::qnorm(1 - (1 - ci)/2)
  df <- df %>% dplyr::mutate(`:=`(!!name, factor(!!name, levels = !!name %>%
                                                   unique() %>% rev(), ordered = TRUE)), .xmin = !!estimate -
                               const * !!se, .xmax = !!estimate + const * !!se, .filled = TRUE,
                             .label = sprintf("%.2f", !!estimate))
  if (logodds) {
    df <- df %>% mutate(.xmin = exp(.data$.xmin), .xmax = exp(.data$.xmax),
                        `:=`(!!estimate, exp(!!estimate)))
  }
  if (!quo_is_null(pvalue)) {
    df <- df %>% dplyr::mutate(.filled = !!pvalue < !!psignif)
  }
  g <- ggplot2::ggplot(df, aes(x = !!estimate, y = !!name))
  if (logodds) {
    if ("xtickbreaks" %in% names(args)) {
      g <- g + scale_x_continuous(trans = "log10", breaks = args$xtickbreaks)
    }
    else {
      g <- g + scale_x_continuous(trans = "log10", breaks = scales::log_breaks(n = 7))
    }
  }
  g <- g + theme_forest() + scale_colour_ng_d() + scale_fill_ng_d() +
    geom_stripes() + geom_vline(xintercept = ifelse(test = logodds,
                                                    yes = 1, no = 0), linetype = "solid", size = 0.4, colour = "black")
  g <- g + geom_effect(ggplot2::aes(xmin = .data$.xmin, xmax = .data$.xmax,
                                    colour = !!colour, shape = !!shape, filled = .data$.filled),
                       position = ggstance::position_dodgev(height = 0.5),
                       fatten=6) +
    ggplot2::scale_shape_manual(values = c(21L, 22L, 23L,
                                           24L, 25L)) + guides(colour = guide_legend(reverse = TRUE),
                                                               shape = guide_legend(reverse = TRUE))
  if ("title" %in% names(args)) {
    g <- g + labs(title = args$title)
  }
  if ("subtitle" %in% names(args)) {
    g <- g + labs(subtitle = args$subtitle)
  }
  if ("caption" %in% names(args)) {
    g <- g + labs(caption = args$caption)
  }
  if ("xlab" %in% names(args)) {
    g <- g + labs(x = args$xlab)
  }
  if (!"ylab" %in% names(args)) {
    args$ylab <- ""
  }
  g <- g + labs(y = args$ylab)
  if ("xlim" %in% names(args)) {
    g <- g + coord_cartesian(xlim = args$xlim)
  }
  if ("ylim" %in% names(args)) {
    g <- g + ylim(args$ylim)
  }
  if ("xtickbreaks" %in% names(args) & !logodds) {
    g <- g + scale_x_continuous(breaks = args$xtickbreaks)
  }
  g
}
environment(forestplot) <- environment(ggforestplot::forestplot)

Thanks a lot, Greg!!! It was really helpful. I am really indebted to you as I was looking for it for last 3 days and was almost hopeless. You saved me!!! Many thanks!!!


text:
What are common ways to calculate gene length for TPM calculations?
I am aware of this similar question. But the accepted answer there answers how to calculate TPM given a mapping from gene name to gene length.
My question is, given an annotation file of genes (for example from Ensembl here), how do you calculate gene length for the purpose of calculating TPM on a gene count matrix coming from RNA seq data. That is, I understand the actualy calculation of TPM given gene lengths, but I do not understand what the gene lengths are.
From what I understand, to calculate gene lengths from annotation file for RNA seq data, you only consider the exons of each gene, but even so, many exons overlap. What is the standard method of dealing with this?
For example, I have read that using the length of the union of all exons of the gene is acceptable, but from what I understand TPM is meant to represent the "expected number of reads coming from gene X when randomly selecting 1e6 transcripts". So taking the union of all exons will be an estimate of this but it might under-estimate this.
EDIT: From what I can tell this question is not a duplicate of RNASeq: Normalization, stabilization, gene length and rlog
That question does discuss how to perform normalization that takes gene length into account, but it does not ask "what is gene length? How does one determine the length of each gene?" Possible answers could be like what I already mentioned: length of the union of all exons in the gene, or another possible answer is one from the linked question "length of the longest transcript isoform" but I feel that it is not really discussed in that answer/question.

I think the most common way is to use something like Kallisto which figures out what proportion of reads are likely assigned to each transcript variant, and calculated TPM  based on that.


text:
How to know if ID correspond to nucleotide sequence or protein with a script
I have a text file that contains a list of IDs (314 sequences):
AVP78031.1
AVP78042.1
ATO98108.1
ATO98120.1
ATO98132.1
...

My goal is to make a script (maybe using Python or Perl) to check in the list if all the IDs are nucleotide or protein sequences.
For example:
AVP78031.1 -> protein (this is a nucleotide sequence, I change nucleotide for protein to show an example).
AVP78042.1 -> nucleotide
ATO98108.1 -> nucleotide
ATO98120.1 -> nucleotide
ATO98132.1 -> nucleotide

Any idea to do a script?
Thank everybody!

If these are all GenBank or RefSeq accessions, you can use Entrez Direct for this as shown below:
$ cat accs.txt
ATO98108.1
ATO98120.1
ATO98132.1
AVP78031.1
AVP78042.1
$ cat accs.txt | epost -db nuccore | efetch -format acc
## no output because none of them are nucleotide accessions
$ cat accs.txt | epost -db protein -format acc | efetch -format acc
AVP78042.1
AVP78031.1
ATO98132.1
ATO98120.1
ATO98108.1

NOTE: This will work only if the accessions are currently live because epost does not find any suppressed accessions. For example:
$ cat accs.txt
NM_002826.3
NM_002826.4
NM_002826.5
$ cat accs.txt | epost -db nuccore -format acc | efetch -format acc
NM_002826.5

Here, all three accessions are valid nucleotide accessions but only the last one, NM_002826.5, is alive.
An alternate way is to use the accession prefixes defined here and come up with an appropriate regular expression query.

The alternate way (accession prefixes) is a great idea!


text:
Why Do I Get Different Results Each Time From AutoDock Vina?
I just tried to use Vina for calculating binding affinity.
I used the same protein, the same ligand, the same options, but there are different results.
Why do I get different results?
Is it the limit Point of docking tools?
Then why we use it?

Vina uses a non-deterministic docking algorithm. See the manual.

The docking algorithm is non-deterministic. Even though with this receptor-ligand pair, the minimum of the scoring function corresponds to the correct conformation, the docking algorithm sometimes fails to find it. Try several times and see for yourself. Note that the probability of failing to find the mininum may be different with a different system.

The manual states that you can set the random seed in order to obtain reproducibility.

I changed something, and now the docking results are different. Why?
Firstly, had you not changed anything, some results could have been different anyway, due to the non-deterministic nature of the search algorithm. Exact reproducibility can be assured by supplying the same random seed to both calculations, but only if all other inputs and parameters are the same as well. Even minor changes to the input can have an effect similar to a new random seed. What does make sense discussing are the statistical properties of the calculations: e.g. "with the new protonation state, Vina is much less likely to find the correct docked conformation".

The answer by @Throckmorton is 100% correct in that it's not deterministic. I just wanted to address a few things —so this is a long comment.
All too often the box is set as the whole protein, this is bound to make it trip up on local minima. Sure, with native ligands, binding pockets close to the active site increases the dissociation constant of the complex, but docking is not the way to calculate these. Generally this is done by tethered MD, where the ligand is pulled out.
Docking returns generally either a binding potential (∆∆G_bind) as a kcal/mol (Autodock) or kJ/mol, where negative is good, or a rank. This is not the typically binding affinity measurements, e.g. an IC50 or dissociation constant (K_D) —The Eyring equation converted in a relative form (ratio makes constants disappear) allows you do some brute conversions to a k-on for the site, but this is highly not recommended.
Docking is not perfect and, worse, gets abused a lot. However I'd flip the OP's statement around and say the limitation by far is wet lab side, where doing an ITC properly is hard and crystallising a protein with a ligand is even worse...