text:
Under what circumstances can the SEQ field in a SAM file be a *?
I saw in the SAM format specifications that the SEQ field (10th column) can be a "*" if the sequence is not stored, instead of being the sequence of the mapped read. Under what circumstances is this expected to occur? I don't understand why it's allowed in the spec.
For a concrete example, I have pacbio data mapped to a reference sequence, and some of the records in the SAM file have a sequence of *, and I don't know why or if that means anything. I don't understand why they exist, since the data is ostensibly present in the FASTQ file I use as input to the read mapper, and I don't understand what purpose they could have, since it seems to me that these records couldn't be interpreted by variant calling software. At best it would seem they can be used to infer coverage, but it's coverage by sequences of an uncertain nature?

Related: How to output all sequences with bwa mem, not *?

Sounds like I should check if these are exclusively secondary alignments

This could be a duplicate question as @acvill is indicating


text:
How to rescale a function?
When we have a function which takes values in [a,b] and which returns values in [c,d], what transformation can we do so that the function takes values in [0,1] and returns values in [0,1] without the curve changing in both cases?

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

Please provide a concrete example that doesn't involve something like "function(x) return x;", which is what you're implying.

Assuming you want a linear transformation (which is what "without the curve changing" sounds like to me), to rescale from range c..d to range 0..1 subtract the minimum then divide by the difference, i.e.
$val_t = (val - c) / (d-c)$
The inverse of this can be applied to work the other way, i.e.
$inp_t = val \times (b-a) + a$
For a function $f(x)$, with ranges as you have stated, the full transformation of input and output would be:
$g(x_t) = \frac{f(x_t \times (b-a) + a) - c}{d-c}$

Very cool response

source: Stack Exchange

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

sort: votes


text:
Comparing differential expression across samples - is batch effect correction needed?
I have a bulk RNA-seq dataset made up of control and treatment conditions for a range of cell lines. This dataset was generated in two batches, such that the cell lines are split between batches but both the treatment and control for each cell line are within the same batch. As the cell lines are very different, I'm looking to do DE analysis for each cell line and then compare these DE responses between cell lines.
I'm wondering whether batch effect correction is useful here since I'm only ever comparing DE values and not the raw counts themselves. Also, the baseline expression values for the cells lines are so different that they already constitute a strong biological batch effect that is perfectly correlated with the technical batch effect, which batch correction wouldn't be able to account for.

There is no need for batch correction because all your comparisons are within batch. As long as you are making comparisons for each line, and treatment and control for the same line are in the same batch, then any batch effects are already accounted. for.

If your treatment and control are prepared in different batches, it's impossible to distinguish between treatment-specific effects, and batch-specific effects. DESeq2 will complain when you attempt to carry out such correction in the statistical model for this reason.
Baseline differences in expression are a large reason for carrying out batch correction in the first place; it makes differential expression analyses more powerful, not less. I'd recommend consulting with a statistician (prior to carrying out an experiment / sequencing) regarding the best way to split samples into different batches so that batch effects can be discovered and accounted for. As you mention, there are multiple effects that should be considered (e.g. cell lines, sequencing run).

Could you expand a bit on why batch correction helps with baseline differences? For example, if I have multiple cell lines and treatments from a single experiment, sequenced in the same run, are you saying that I should still do batch correction to adjust for the baseline expression? Doesn't that run the risk of altering real biological differences?

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


text:
What is the difference between flat and stereo stitch id?
For my thesis, I'm working on side effect prediction in a machine learning framework and I am working on data available on SIDER (http://sideeffects.embl.de/download/). If you open the README file, the first description you'll see is:

meddra_all_se.tsv.gz
1 & 2: STITCH compound ids (flat/stereo, see above)
3: UMLS concept id as it was found on the label
4: MedDRA concept type (LLT = lowest level term, PT = preferred term; in a few cases the term is >neither LLT nor PT)
5: UMLS concept id for MedDRA term
6: side effect name

I can not really understand what is the difference between flat and stereo compound ID; is it a chemical descriptor? Is it a notation issue? Moreover, I thought that the two IDS for the same compound could yield to similar results on PubChem but they don't. For example, in the dataset the STITCH ids for the first compound are CID100000085 and CID000010917; on PubChem the first is an aromatic compound and the second is not, they are widely different.
Thanks to everyone who will take the time to help me.

Hey, welcome to Bioinformatics. I see you asked the same question on another SO! Please try to avoid posting the question on multiple communities.

As for your questions, I may be wrong but I think those terms refer to the chemistry 3D representation in space. I think Stereo should stand for Stereoisomers. You can read more about it here. As for flat, it is how you would represent a compound on a 2D space.
The same compound arranged differently in 3D space can be different molecules and thus yield different results.


text:
plot the relation between RMSD and free energy
how can i plot the relation between rmsd and free energy like the image i attached as i need to find the most stable conformation structure for my cycle peptide after finishing REMD simulation

I look at that graph and I note the high ∆G at 0 Ångström and wonder if the reference pose for the RMSD was not energy minimised correctly... so maybe not the best reference. http://prody.csb.pitt.edu/ is a helpful toolkit for this kind of analysis. Is Python an option and what is your reference? (If any user knowledgeable about peptides sees this please step in as I doubt I can help —NMR works by magic, right?)

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

Thanks for the comment .. i'm trying to use the tool but couldn't .. in which part please i can plot RMSD with energy surface ?

You would need to calculate the free energy along a trajectory and the rmsd along that trajectory and then cross-correlate/eliminate variables.  There may be a utility which does this but I am not ware of it.
https://manual.gromacs.org/current/user-guide/cmdline.html
There are two utilities of interest for doing just that.  See gmx_energy and gmx-rmsd.  Then you can plot one against the other by elimination of the time/space variables. Check out other calculations you can do with GROMACS:
https://manual.gromacs.org/current/user-guide/cmdline.html.


text:
How to combine multiple .fasta files of primary assembly from Ensembl into one for sequence alignment?
I have some marmoset snRNA reads that I want to align with the reference transcriptome using cellranger. The primary assembly for marmoset is available here, which is broken down into 22 parts. However, cellranger mkref only accepts one .fa file to generate the transcriptome.
I tried concatentaing all the extracted .fa files into one (After removing the first line which is a comment), but I get the following error while trying to make a reference transcriptome using cellranger:
Creating new reference folder at /data/users/abhishek/ensembl/common_marmoset_ref/combined/common_marmoset
 ...done

 Writing genome FASTA file into reference folder...
 ...done

 Indexing genome FASTA file...
 Traceback (most recent call last):
   File "/data/public-data/cellranger-6.1.1/bin/rna/mkref", line 137, in <module>
     main()
   File "/data/public-data/cellranger-6.1.1/bin/rna/mkref", line 128, in main
     reference_builder.build_gex_reference()
   File "/data/public-data/cellranger-6.1.1/lib/python/cellranger/reference.py", line 630, in build_gex_reference
     self.process_fasta()
   File "/data/public-data/cellranger-6.1.1/lib/python/cellranger/reference.py", line 529, in process_fasta
     subprocess.check_call(["samtools", "faidx", self.fasta_path])
   File "/data/public-data/cellranger-6.1.1/external/anaconda/lib/python3.8/subprocess.py", line 364, in check_call
     raise CalledProcessError(retcode, cmd)
 subprocess.CalledProcessError: Command '['samtools', 'faidx', '/data/users/abhishek/ensembl/common_marmoset_ref/combined/common_marmoset/fasta/genome.fa']' returned non-zero exit     status 1.

There are also files which contain only one or two characters. I suspect that these aren't being handled properly with the direct concatenation of all the files. Some examples are below:
==> Callithrix_jacchus.mCalJac1.pat.X.dna.primary_assembly.1.fa <==
CCCATTATCCCTGCTTCTGCTAATCTTTCAAAATCTCACCACCATAGTC

==> Callithrix_jacchus.mCalJac1.pat.X.dna.primary_assembly.20.fa <==
TTTTT

==> Callithrix_jacchus.mCalJac1.pat.X.dna.primary_assembly.21.fa <==
TT

==> Callithrix_jacchus.mCalJac1.pat.X.dna.primary_assembly.22.fa <==
AAAATTGTTTCTTGTTT

==> Callithrix_jacchus.mCalJac1.pat.X.dna.primary_assembly.2.fa <==
C

==> Callithrix_jacchus.mCalJac1.pat.X.dna.primary_assembly.3.fa <==
TTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGG

How do I combine all these .fa files to generate the reference transcriptome for downstream analysis?

When you say "After removing the first line which is a comment", do you mean you removed the lines which start with a >? Those are the headers and are required to have a valid fasta file. They contain the name or identifier of the contig/chromosome/transcript. They separate the contigs in the fasta if there are multiple and will be needed for alignment in order to know which contigs your reads mapped to.

Why do not you run "samtools faidx" alone on the fasta file any see if you are getting any useful error message. Looks like your program is calling samtools for indexing. You can tell this from the last line of the error.

And your fasta file does not look correct.


text:
BEAGLE do not see and/or not use OpenCL
I am conducting phylogeny analyses using BEAST 2 as software tool. It is said that it can be accelerated by BEAGLE-library. I have a personal laptop with Windows & a cluster with CentOS.
I installed both BEAST 2 & BEAGLE onto them. This gives me a boost by 3.6 times. It was said too that using OpenCL with BEAGLE gives me more speed-up. Nonetheless, my calculations takes from twelve hours to 3 days so I got really interested in this. Installing it on Windows were very simple. When I ran beast -beagle_info, I saw:
                        BEAST v2.6.6, 2002-2021
             Bayesian Evolutionary Analysis Sampling Trees
                       Designed and developed by
 Remco Bouckaert, Alexei J. Drummond, Andrew Rambaut & Marc A. Suchard

                   Centre for Computational Evolution
                         University of Auckland
                       r.bouckaert@auckland.ac.nz
                        alexei@cs.auckland.ac.nz

                   Institute of Evolutionary Biology
                        University of Edinburgh
                           a.rambaut@ed.ac.uk

                    David Geffen School of Medicine
                 University of California, Los Angeles
                           msuchard@ucla.edu

                      Downloads, Help & Resources:
                           http://beast2.org/

  Source code distributed under the GNU Lesser General Public License:
                   http://github.com/CompEvol/beast2

                           BEAST developers:
   Alex Alekseyenko, Trevor Bedford, Erik Bloomquist, Joseph Heled,
 Sebastian Hoehna, Denise Kuehnert, Philippe Lemey, Wai Lok Sibon Li,
Gerton Lunter, Sidney Markowitz, Vladimir Minin, Michael Defoin Platel,
          Oliver Pybus, Tim Vaughan, Chieh-Hsi Wu, Walter Xie

                               Thanks to:
          Roald Forsberg, Beth Shapiro and Korbinian Strimmer

--- BEAGLE RESOURCES ---

0 : CPU
    Flags: PRECISION_SINGLE PRECISION_DOUBLE COMPUTATION_SYNCH EIGEN_REAL EIGEN_COMPLEX SCALING_MANUAL SCALING_AUTO SCALING_ALWAYS SCALERS_RAW SCALERS_LOG VECTOR_SSE VECTOR_NONE THREADING_NONE PROCESSOR_CPU FRAMEWORK_CPU

1 : NVIDIA GeForce 940MX
    Global memory (MB): 2048
    Clock speed (Ghz): 1.19
    Number of cores: 384
    Flags: PRECISION_SINGLE PRECISION_DOUBLE COMPUTATION_SYNCH COMPUTATION_ASYNCH EIGEN_REAL EIGEN_COMPLEX SCALING_MANUAL SCALING_AUTO SCALING_ALWAYS SCALERS_RAW SCALERS_LOG VECTOR_NONE THREADING_NONE PROCESSOR_GPU FRAMEWORK_CUDA

2 : Intel(R) HD Graphics 620 (OpenCL 2.1 )
    Global memory (MB): 3219
    Clock speed (Ghz): 1.00
    Number of compute units: 24
    Flags: PRECISION_SINGLE PRECISION_DOUBLE COMPUTATION_SYNCH COMPUTATION_ASYNCH EIGEN_REAL EIGEN_COMPLEX SCALING_MANUAL SCALING_AUTO SCALING_ALWAYS SCALERS_RAW SCALERS_LOG VECTOR_NONE THREADING_NONE PROCESSOR_GPU FRAMEWORK_OPENCL

3 : Intel(R) Core(TM) i5-7200U CPU @ 2.50GHz (OpenCL 2.1 (Build 10))
    Global memory (MB): 8067
    Clock speed (Ghz): 2.50
    Number of compute units: 4
    Flags: PRECISION_SINGLE PRECISION_DOUBLE COMPUTATION_SYNCH COMPUTATION_ASYNCH EIGEN_REAL EIGEN_COMPLEX SCALING_MANUAL SCALING_AUTO SCALING_ALWAYS SCALERS_RAW SCALERS_LOG VECTOR_NONE THREADING_NONE PROCESSOR_CPU FRAMEWORK_OPENCL

However, it still uses only resourse #0, even through forced beast -beagle_order 3 -threads 1 filename.xml. Although, in comparison to Linux case, it was more successful.
I have downloaded intel_sdk_for_opencl_applications_2020.3.494.tar.gz from here. Then I unzipped it & run install.sh. Installation was successful. After that, I have decided to go to build folder in beagle-lib & run cmake -DCMAKE_INSTALL_PREFIX:PATH=$HOME .. in order to go on to run make install but I get the next message:
-- JAVA_HOME=
-- JNI_INCLUDE_DIRS=/usr/lib/jvm/java/include;/usr/lib/jvm/java/include/linux;/usr/lib/jvm/java/include
-- JNI_LIBRARIES=/usr/lib/jvm/java/lib/libjawt.so;/usr/lib/jvm/java/lib/server/libjvm.so
-- Not using libtools for plugins
-- Could NOT find OpenCL (missing: OpenCL_LIBRARY OpenCL_INCLUDE_DIR)
CUDA_TOOLKIT_ROOT_DIR not found or specified
-- Could NOT find CUDA (missing: CUDA_TOOLKIT_ROOT_DIR CUDA_NVCC_EXECUTABLE CUDA_INCLUDE_DIRS CUDA_CUDART_LIBRARY)
-- Configuring done
-- Generating done
-- Build files have been written to: /export/home/rinkman/beagle-lib/build

It was solved through manual specification of directories cmake -DOpenCL_INCLUDE_DIR=~/intel/system_studio_2020/opencl/SDK/include/ -DOpenCL_LIBRARY=~/intel/system_studio_2020/opencl/SDK/lib64/libOpenCL.so.1.2 -DCMAKE_INSTALL_PREFIX:PATH=$HOME ... The compiling was successful after that. make test was executed without any problems. I detected the file libhmsbeagle-opencl.so.40.0.0 in ~/lib folder. But beast -beafle_info did not show OpenCL. As there are older versions of BEAST & BEAGLE on the root, I tried to change very roughly startup script to run only my local BEAGLE but nothing changed:
#!/bin/sh

if [ -z "$BEAST" ]; then
    ## resolve links - $0 may be a link to application
    PRG="$0"

    # need this for relative symlinks
    while [ -h "$PRG" ] ; do
        ls=`ls -ld "$PRG"`
        link=`expr "$ls" : '.*-> \(.*\)$'`
        if expr "$link" : '/.*' > /dev/null; then
        PRG="$link"
        else
        PRG="`dirname "$PRG"`/$link"
        fi
    done

    # make it fully qualified
    saveddir=`pwd`
    BEAST0=`dirname "$PRG"`/..
    BEAST=`cd "$BEAST0" && pwd`
    cd "$saveddir"
fi

BEAST_LIB="$BEAST/lib"

if [ -z "$JAVA_HOME" ]; then
  JAVA=java
else
  JAVA="$JAVA_HOME"/bin/java
fi

# use BEAGLE_LIB if the BEAGLE library is not in a standard location
if [ -n "$BEAGLE_LIB" ]; then
    if [ -n "$BEAST_EXTRA_LIBS" ]; then
        BEAST_EXTRA_LIBS="$BEAST_EXTRA_LIBS:$BEAGLE_LIB"
    else
        BEAST_EXTRA_LIBS="$BEAGLE_LIB"
    fi
fi

# Explicitly add /usr/local/lib to library search path to ensure
# beast continues to find beagle when installed here.  (This is
# necessary due to known problems with certain JREs.)
if [ -z "$LD_LIBRARY_PATH" ]; then
    export LD_LIBRARY_PATH=/usr/local/lib
else
    export LD_LIBRARY_PATH="$LD_LIBRARY_PATH":/usr/local/lib
fi

BEAST_EXTRA_LIBS="$HOME/lib/"

# use BEAST_EXTRA_LIBS variable to load BEAGLE and other libraries from non-default locations
# this assumes that the library path contains all these libraries (or are set through LD_LIBRARY_PATH)
#if [ -n "$BEAST_EXTRA_LIBS" ]; then
#  "$JAVA" -Dlauncher.wait.for.exit=true -Xms256m -Xmx1g -Djava.library.path="$BEAST_EXTRA_LIBS" -Duser.language=en -cp "$BEAST_LIB/launcher.jar" beast.app.beastapp.BeastLauncher $*
#else
#  "$JAVA" -Dlauncher.wait.for.exit=true -Xms256m -Xmx1g -Duser.language=en -cp "$BEAST_LIB/launcher.jar" beast.app.beastapp.BeastLauncher $*
#fi
"$JAVA" -Dlauncher.wait.for.exit=true -Xms256m -Xmx1g -Djava.library.path="$BEAST_EXTRA_LIBS" -Duser.language=en -cp "$BEAST_LIB/launcher.jar" beast.app.beastapp.BeastLauncher $*

Please, could anyone help with this?
P.S. I am a noviciate in unix-systems.

@M__ I installed BEAGLE successfully per se & it does work but I would like to run it with OpenCL.

@M__ Do you mean beast-users?

I better explain this style of parameter rich MCMC takes ages and the issues are how speed this up.

Having given this some thought my suggestion is to install miniconda or anaconda and then install the port to OpenCL, then install beast and beagle.
Anaconda is here: https://www.anaconda.com/products/individual
Then install openCL via conda:
 conda create -n beasty python=3.9
 conda activate beasty
 conda install -c conda-forge pyopencl
 conda install -c bioconda beagle
 conda install -c bioconda beast2

My idea was you'd access OpenCL via python described here. This wouldn't work ... its complicated.
Traditionally OpenCL and Beagle didn't work ... on an Intel Mac

2 : Intel(R) HD Graphics 620 (OpenCL 2.1 )

However, if this is your chip set - I wouldn't give up. In any other scenario  I wouldn't pursue this - personally, below.
One alternative is NVIDA e.g. Cloud computing used via CUDA and Beagle. CUDA will easily install and Beagle stacks straight on top.
Any heavy duty MCMC is very amenable to traditional CPU parallelisation - if you think at it - it would be and thats why GPU - CUDA will work. Personally, looking at your resources I would parallelise under you cluster using a MINIMUM of 8-cores per calculation outside CUDA/OpenCL. If that has not already been done, that is the solution in context in my opinion. Unless you're a Google developer few will have access to a NVIDA cluster in any case sufficient CPU parallelisation should be able to achieve the runtime close to that of a good GPU (its gonna be a lot of CPUs though).
Conclusion
Honestly, my personal advise is that phylogenetics sump clusters no matter what the calculation is being used and that has always been the case. The best solution is to max out on available cores. When eukaryote genome guys complain about 1000 hour runtimes (obviously thats the total CPU parallelised), thats a routine calculation in phylogenetics. Buying in NVIDA doesn't necessarily solve this, because whilst MCMC is massively accelerated for maximum likelihood isn't accelerated to that extent by comparison.
Extra information In case the OP is not aware, the enormous advantage of Beast2 is the ability to resume an MCMC at the point it ended (or some checkpoint therein), you don't need to run the whole thing again. You should seek to exploit that, so a calculation started on a laptop can be finished on a cluster. I can't state sufficiently how powerful that is, it stands to reason because the random number stream is declared and every sample thereafter. Keeping all calculations is standard practice in any case, but especially with a Beast2 run and don't every discard it (hard disk is cheap right?). You might find you have to return to it and its a whole lot better than needing to do the whole thing again.
Its time to move on.

I have gotten `Found conflicts! Looking for incompatible packages.
This can take several minutes.  Press CTRL-C to abor
failed
UnsatisfiableError: The following specifications were found to be incompatible with each other:
Output in format: Requested package -> Available versions
Package libstdcxx-ng conflicts for:
python=3.10 -> libffi[version='>=3.4.2,<3.5.0a0'] -> libstdcxx-ng[version='>=7.3.0|>=7.5.0']`

`beast2 -> beagle-lib -> libstdcxx-ng[version='>=10.3.0|>=9.3.0|>=7.5.0|>=7.3.0']The following specifications were found to be incompatible with your system:

feature:/linux-64::__glibc==2.17=0
python=3.10 -> libgcc-ng[version='>=7.5.0'] -> __glibc[version='>=2.17']

Your installed version is: 2.17`
What does this mean? Which packages should I delete? Of what is 2.17 a version? Please could you help?

I get the same conflict with a spanking new environment :(

beast2 can be installed with beagle & opencl only through adding conda-forge & bioconda channels, & creating new environment from scratch as conda create -n env_name opencl beagle beast2. Nonetheless, conda run beast -beagle_info shows only 1 resource: CPU without OpenCL

There is a way to save states in BEAST 1. Should it be left here as a comment or as the answer correction?

That is certainly interesting @Vovin , no its more than a comment. A separate answer would be cool, because I'm not aware of this and it is important or I can ask a separate question.

In this case, I think I could made a Q&A according to this

Thats perfectly, look forward to this @Vovin


text:
Biopython reads my tree eternally long
I have a nexus tree (1332 taxa) with a lot of additional data. When I tried to read it through tree = Phylo.read(treepath, "nexus"), my kernel got eternally loaded.
If I abort the process, I get the following message:
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Input In [95], in <cell line: 1>()
----> 1 tree = Phylo.read(r"D:\Downloads_hard_disk\Job\NAs_Hepatitis_A\Results\experiments\test_right\test_right_ExpPop_3\test_right_ExpPop_3.biopy.tree", "nexus")

File ~\miniconda3\lib\site-packages\Bio\Phylo\_io.py:60, in read(file, format, **kwargs)
     58 try:
     59     tree_gen = parse(file, format, **kwargs)
---> 60     tree = next(tree_gen)
     61 except StopIteration:
     62     raise ValueError("There are no trees in this file.") from None

File ~\miniconda3\lib\site-packages\Bio\Phylo\_io.py:49, in parse(file, format, **kwargs)
     34 """Parse a file iteratively, and yield each of the trees it contains.
     35
     36 If a file only contains one tree, this still returns an iterable object that
   (...)
     46
     47 """
     48 with File.as_handle(file) as fp:
---> 49     yield from getattr(supported_formats[format], "parse")(fp, **kwargs)

File ~\miniconda3\lib\site-packages\Bio\Phylo\NexusIO.py:40, in parse(handle)
     32 def parse(handle):
     33     """Parse the trees in a Nexus file.
     34
     35     Uses the old Nexus.Trees parser to extract the trees, converts them back to
   (...)
     38     eventually change Nexus to use the new NewickIO parser directly.)
     39     """
---> 40     nex = Nexus.Nexus(handle)
     42     # NB: Once Nexus.Trees is modified to use Tree.Newick objects, do this:
     43     # return iter(nex.trees)
     44     # Until then, convert the Nexus.Trees.Tree object hierarchy:
     45     def node2clade(nxtree, node):

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:668, in Nexus.__init__(self, input)
    665 self.options["gapmode"] = "missing"
    667 if input:
--> 668     self.read(input)
    669 else:
    670     self.read(DEFAULTNEXUS)

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:718, in Nexus.read(self, input)
    716     break
    717 if title in KNOWN_NEXUS_BLOCKS:
--> 718     self._parse_nexus_block(title, contents)
    719 else:
    720     self._unknown_nexus_block(title, contents)

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:759, in Nexus._parse_nexus_block(self, title, contents)
    757 for line in block.commandlines:
    758     try:
--> 759         getattr(self, "_" + line.command)(line.options)
    760     except AttributeError:
    761         raise NexusError("Unknown command: %s " % line.command) from None

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:1181, in Nexus._tree(self, options)
   1179     elif special == "W":
   1180         weight = float(value)
-> 1181 tree = Tree(name=name, weight=weight, rooted=rooted, tree=opts.rest().strip())
   1182 # if there's an active translation table, translate
   1183 if self.translate:

File ~\miniconda3\lib\site-packages\Bio\Nexus\Trees.py:82, in Tree.__init__(self, tree, weight, rooted, name, data, values_are_support, max_support)
     80 # there's discrepancy whether newick allows semicolons et the end
     81 tree = tree.rstrip(";")
---> 82 subtree_info, base_info = self._parse(tree)
     83 root.data = self._add_nodedata(root.data, [[], base_info])
     84 self._add_subtree(parent_id=root.id, tree=subtree_info)

File ~\miniconda3\lib\site-packages\Bio\Nexus\Trees.py:123, in Tree._parse(self, tree)
    121 elif not incomment and tree[p] == ")":
    122     plevel -= 1
--> 123 elif tree[p:].startswith(NODECOMMENT_START):
    124     incomment = True
    125 elif incomment and tree[p] == NODECOMMENT_END:

KeyboardInterrupt:

I am interested what is the specific reason of this: either the ineffective algorithm or a bug in Biopython? Is there a way to read the tree as Bio.Phylo.Newick.Tree object? Or should I completely use another package? In the last case, could you recommend effective one?

The reason is the RAM has been exceeded for Biopython's capabilities and you are encountering a 'RAM-bottleneck'. I've listed the options from 1 or 3 and recommend option 3c.
The options are use:

machine with much greater RAM
acceleration coding techniques in Python (such parallelisation)
Beast tools written in Java

I assume this calculation is being performed on your laptop because in theory <1500 isn't much, but agreed there are ~1500 taxa.
Option 1 I personally don't think you are very far into the bottle-neck. Certainly a Cloud computer (e.g. GCP) could easily perform this calculation and maybe your the headnode of your local cluster. Usually the headnode has the most RAM.
Option 2 Code acceleration is worth keeping in mind because this achieves >> 1 log performance. Parallisation will solve this using an optimal chunk-size. However, large amounts of code can be needed. Code optimisation - Python is very sensitive to this because it uses C-libraries - for example, use map and reduce and avoid ANY regex in your Python code because IMO its a key weakness of Python. Pandas iterrows is famously slow. If the Biopython library is the problem code efficiency isn't an option (unless you want to rewrite Biopython). Label compression or header recoding techniques work, I'm not sure whether this has been done in Beast output already(?).
Option 3 The tools for dealing with this which are written and optimised to achieve rapid output are:

Tracer - does lots of analysis on MCMC output. Its really geared towards assessing convergence, but it can do a lot of stuff.
TreeAnnotator - makes a summary tree from the MCMC output you have

TreeStat2 <- I think this is what you are looking for.


text:
What does this accession NCBI code mean: 6MWN_B?
According to this article, accession codes should consist from a combination of uppercase letters following a combination of digits. If this is a RefSeq, it can have a prefix as a combination of uppercase letters with underscore.
But this accession code 6MWN_B does not consist with all aforementioned. What is it? I could not find an answer via Google.

In general, an accession is a value that points to (provides access to) a single database entry. Different databases (GenBank, RefSeq, Ensembl, SwissProt, and so on) have different conventions for accessions.
The accession in question points to an RNA sequence associated with a Protein Data Bank (PDB) record: "Crystal structure of hepatitis A virus IRES domain V in complex with Fab HAVx". In the PDB record, you can click the "sequence" tab and select chain "B" from the drop-down, and you will see an RNA sequence (ACGU) matching the DNA sequence (ACGT) in the NCBI record to which you linked.
I'm not sure what the rules or conventions for PDB accessions are, but they are clearly different from RefSeq accessions.

Thank you! I am just surprised that this lies in nuccore & did not get a new local accession code.


text:
Handling failed blocks within scattered WDL workflow
I'm trying to adapt an existing workflow to make it more robust to failures. Below is a simplified example of the workflow.
workflow main {
   File filelist="list_of_input_files.txt"
   Array[String] files=read_lines(filelist)
   scatter(file in files) {
      call task1 { input: in=file }
      if(defined(task1.out)) { call task2 { input: in=task1.out } }
      if(defined(task2.out)) { call task2 { input: in=task2.out }
   }
}

task task1 {
   File in
   String out1 = basename(in,'.txt') + '_1.txt'
   command { if(...) { ./script1.sh ${in} > ${out1} } }
   output { File? out = out1 }
   runtime { ... }
}

task task2 {
   File in
   String out2 = basename(in,'.txt') + '_2.txt'
   command { if(...) {./script2.sh ${in} > ${out2} } }
   output { File? out = out2 }
   runtime { ... }
}

task task3 {
   File in
   String out3 = basename(in,'.txt') + '_3.txt'
   command { if(...) {./script3.sh ${in} > ${out3} } }
   output { File? out = out3 }
   runtime {...}
}

with my list_of_input_files.txt looking like
firstfile.txt
secondfile.txt
thirdfile.txt
...

Obviously this is oversimplified, but I am trying to catch failed blocks by making the output of each task optional (using File?) and then having an if(defined(taskX.out)) condition to call the next task. My reasoning is that if something fails (represented by the if(...) within the command block of each task), then the out file will never be created in that task and so won't exist as input for the next task and thus if(defined(taskX.out)) will resolve as false and skip the next task (and therefore all the subsequent tasks).
With current testing, the first file in the list satisfies the if(...) conditions of all the tasks and succeeds in giving me the final output file firstfile_1_2_3.txt. However, I am getting a Required file output 'secondfile_1.txt' does not exist' error message because it fails the if(...) condition in task1 and so the output file is never created.
From what I've read in the WDL v1.1 spec, optional File outputs are allowed and are undefined if the files they point to don't exist. However, after much confused googling, I've realised it's probably because I'm using Cromwell, which is stuck with WDL v1.0.
My question is whether this should be working in WDL v1.0 and, if not, is my logic correct and should this (psuedo-)code work with WDL v1.1? Is there a way to achieve what I'm trying to do using pure WDL v1.0?

I stayed clear of WDL because seemed to over complicate things. Can re-write this as a bash script but would need to study the pipeline a bit more.

I think the thing to do is contact Broad Institute directly because this appears to be a version issue pertaining the upgrade from 1.0 to 1.1

Your WDL code looks more like draft-2 WDL than 1.0 or 1.1 WDL. In WDL 1.0/1.1 there is a separate input section in tasks and workflows, and the version must be specified on the first line. In any case, optional file outputs should be supported in all versions.

@M__ Thanks - my confusion is based around the different specs. It seems that the features I'm trying to use (optional files) don't seem to be supported in the spec that Cromwell uses.

@DavyCats Thanks for the comment! I think I've made the mistake to assume that Cromwell incorporates the latest drafts, but it actually uses WDL 1.0 and so these features are not supported. Despite the WDL 1.1 spec being > 1 year old, I can't see any Cromwell implementation based around it yet.

Optional outputs should work in WDL 1.0. Here is an example that runs without issue through cromwell (version 74).
Note that for the second call in the scatter the input file needs to be converted to from an optional to an non optional value. Since this only gets run if this value is defined, this conversion can safely be done using the select_first function (which returns the first defined value in a given list of optionals).
version 1.0

workflow test {
    input {
        File filelist = "files.txt"
    }
    Array[String] files = read_lines(filelist)

    scatter (file in  files) {
         call blah as b1 {input: in=file}
         if (defined(b1.out)) {
             call blah as b2 {input: in=select_first([b1.out])}
         }
    }
}

task blah {
    input {
        File in
    }
    String outPath = basename(in,'.txt') + '_blah.txt'

    command {
        if $(grep -q "yay" ~{in})
          then
             echo yay > ~{outPath}
        fi
    }

    output {
        File? out = outPath
    }
}

files.txt looks like this:
test1.txt
test2.txt

With test1.txt containing yay and test2.txt containing nay.


text:
Detect mutation context in a read of a sam file
After sequencing (Illumina) of some DNA, I generated a sam file through alignment of a fastq file (using Bowtie2). Instead of using variant calling programs, I want to know the specific variances, compared to the reference genome, in each specific read.
For example, I'm only interested in C to A mutations. I would like to include reads that contain MORE than 5 C to A mutations. Afterwards I would like to analyze the specific sequence context of that mutation and for doing that I need to know the 3 nucleotides before and after the C to A mutation, which gives a context of NNNCNNN.
I looked at samtools mpileup or bcftools, but I can't find the exact thing that I want to analyze. Any ideas?
Thank you!

If you want to look at a few, specific reads manually, you can do this on VarSome. We have recently released a new feature that allows you to enter single NGS or Sanger reads and get the list of variants detected in that read. For example, here:

In this example, there was one variant identified when comparing to the human hg38 genome assembly:

CYP24A1(NM_000782.5):c.1186C>T which is at chr20:54158136 G⇒A (note that the CYP24A1 gene is on the - strand which is why the C to T variant in HGVS notation is a G to A in the genome). You can then click on the variant link(s) to investigate the variant's effect and annotation.

Obviously, this isn't practical for many reads, but it can be handy if you just want to look at a few manually. More importantly, NGS is very error prone (which is why we sequence at high depth), and you should absolutely not assume that just because a variant is present in a single read it was also actually present in the source DNA you sequenced.

Disclaimer: I work for the company that builds VarSome.


text:
Ligand and protein preparation in Vina docking
I am an undergraduate novice, from CS major, intending to run a virtual-screening run using Autodock Vina. A few questions in this regard:

How to choose the best charge model for the ligand?

By default (p. 2), Vina installs Gasteiger-Marsili (GM) partial charges prior to docking. However, GM, being a conformationally-independent approach, is a poor choice when compared to other conformationally-dependent approaches.
In that case, can I choose B3LYP/6-311G/AIM which is the best-performing EEM-parametrization-set, yet published?

What are the upsides and downsides?

How to choose the best charge model for the protein?

Once again, Vina installs Gasteiger-Marsili (GM) partial charges prior to docking and there exist better models.
Is there any better open-source alternative than MOPAC-PM7?

Vina removes non-polar hydrogens and lone pairs from the ligand as well as protein.

Is there any upside to this operation other than decreasing computational complexity?
If time etc. is not a factor, are docking outputs going to be affected to any significant degree in absence of this modification?

@M__ Narrowed down my previous question

Ideally StackExchange questions should be related to a single problem. Can you please give more details about your project / problem? If this is for an assignment (or similar canned problem), please add the tag 'homework'. If you're interested in general discussion rather than a solution to a specific problem, I recommend you instead post to https://reddit.com/r/bioinformatics

Your opening statement forces me to ask: have you done a re-docking test with vanilla settings and what was wrong with it? Or are you well past that step?
Gasteiger partial charges are more than adequate. The drude particle mod in some FF addresses the polarisability problem, but is uncommon. Some people run QM charge calculations in Gaussian, for better partial charges, but most folk don't. The major issues in docking stem from protein —charged active site sidechains, especially involving cysteine and histidine, loops that ought to have been flipped to close up on the complex etc


text:
Wrapping R scripts in WDL and run in docker
I am trying to run an R script as a task in WDL in a docker image. Since this is my first time working with these, I am stuck in something very simple (and in the initial stage) and am hoping somebody here can point me in the right direction.
Here is my dockerfile.
FROM rocker/r-base:latest
RUN Rscript ...r-docker/install_packages.R

Here is my WDL script
version 1.0

workflow Rpipe {
   input {
     String name
   }
  call PrintNameR{
        input:
        name = name
}
  output {
     File outfile = PrintNameR.output_R
  }
}

task PrintNameR {
      input {
           String name
}

command<<<
    R --no-save --args name <<Rscript
    args <- commandArgs(trailingOnly = TRUE)
    for (i in args){
    print(i)}
    Rscript
  >>>

runtime {
            docker: "firsttrydocker"
       }
output {
     File output_R = stdout()
  }
}

I validated, generated an input json file - edited it to include the string I want "Smith", and ran it.
java -jar ../cromwell/womtool-77.jar validate rtry1.wdl
java -jar ..cromwell/womtool-77.jar inputs rtry1.wdl > rtry1.json
java -jar ..cromwell/cromwell-77.jar run rtry1.wdl --inputs rtry1.json

I get an output:
> args <- commandArgs(trailingOnly = TRUE)
> for (i in args){
+ print(i)}
[1] "name"

I get "name", I wanted to get "Smith" because that's what my input is.
What am I doing wrong? Clearly, it seems that I am not passing the WDL task input to the WDL task command -in the Rscript. What is the proper way to do it?

In order to reference a WDL input within a command block you need to use an expression placeholder. Expression placeholders have the form ~{*wdl expression*}, where the inner part can be an valid WDL expression (ie reference to an input declaration)
When an engine sees a String or a command block, it first looks to see if there are any of these expression placeholders within it. If there are, it will evaluate the inner expression, and replace the entire placeholder ~{...} with the resulting value.
You can see how this can be applied below
String first_name = "John"
String last_name = "Smith"
Int age = 45

String full_name = "~{first_name} ~{last_name} # equals "John Smith"
String name_and_age = "John is ~{age}" # equals "John is 45

# Any valid WDL expression can go within an expression placeholder
String complex = "~{first_name + last_name + age}" # equals "JohnSmith45"

A key thing to note is that WDL command templating and expression evaluation is completely separate from the actual execution of the command. An execution engine will evaluate the command template, swap the expression placeholders with their evaluated values, and  THEN, execute the script.
This is relevant for your use case since you are using heredoc to directly define the rscript you are runnning. If we briefly look at a python example, you can see what I mean:
task hello_world {
  input {
   String name = "john"
  }

  command <<<
   python <<CODE
     # this is valid, because once the placeholder evaluates
     # this command becomes name = "john", which is valid python syntax
     # We did not have to pass ~{name} as an arg to the python command
     name = "~{name}"
     print(name)
   CODE
  >>>
}

Now getting to your problem specifically. You are trying to reference the name input value without using any sort of expression placeholder. To the execution engine, the name value looks like any other string. To swap the name in the command for the value passed into the task, you simply need to wrap the variable name in an expression placeholder: ~{name}.
task PrintNameR {
  input {
    String name
  }

  command<<<
    # I would wrap this in quotes to prevent any side effects of
    # shell expansion
    R --no-save --args '~{name}' <<Rscript
    args <- commandArgs(trailingOnly = TRUE)
    for (i in args){
    print(i)}
    Rscript
  >>>

  runtime {
            docker: "firsttrydocker"
  }

  output {
     File output_R = stdout()
  }
}
```

I figured I had to pass ~{sep=" " name} as an argument
command<<<
R --no-save --args ~{sep=" " name} <<Rscript
args <- commandArgs(trailingOnly = TRUE)
for (i in args){
print(i)}
Rscript


text:
Low mapping reads for RNA sequence
I am new to RNA sequencing analysis. I have RNA samples for bacteria for different treatment. I did contaminant filtration, remove adapters and checked fastQC report: attached below for one sample (fwd and rev) .  After mapping with the reference genome using bbsplit but the mapping is only 0.62%.
Any help/suggestions highly appreciated.
Thanks

Well, when you looked at the sequences of some of your reads, what species did they BLAST to?

In addition to swbarnes2 comment, I'd suggest a more systematic approach.
Kraken2 is a tool for metagenomic analysis of reads, though it can be hard to set up.
mash screen is another method for metagenomic analysis, much easier to set up but noisier.
If you are dealing with contamination, treating your data as metagenomic might help.
Also, as always, ensure that your workflow and reference assembly are correct!!!


text:
LD Clumping vs. Conditional Analysis for independent associations
I am working with some GWAS summary statistics. My ultimate goal is to extract variants to use as instrument variables for Mendelian Randomization analysis, which requires me to find SNPs that are independently associated with the exposure variable of interest.
My current workflow involves extracting SNPs from summary statistics with P-value < 5E-08, then performing LD Clumping with a stringent R^2 cutoff (R^2<0.05) within a locus. Is this sufficient to ensure that my variants are independently associated, or should I also plan on running conditional analysis(e.g. GCTA-COJO) in addition to/in lieu of clumping for my variants?

Please clarify your specific problem or provide additional details to highlight exactly what you need. Input/output examples are really helpful. As it's currently written, it's hard to tell exactly what you're asking.

LD Clumping with $r^2$ filtering seems like a good start, but bear in mind that population structure has a significant effect on association statistics.
I recommend carrying out separate PCAs (or some other dimensional reduction) on your variable groups (e.g. risk/non-risk, or high/medium/low) to make sure that their genomic population structure is similar. Alternatively / equivalently, carry out a whole-group reduction, and confirm that the coordinate centres for each variable group are similar.
If the population structures are not similar (which is very likely), then it's important to [also] exclude structure-specific SNPs from the analysis, as these will confound the results of your analysis. Use a clustering method to identify people with similar genetic structure, and treat the different clusters as different trait groupings for identifying those SNPs.


text:
How can you put a bash array into a wdl variable?
Given two separate wdl arrays, my command section will do some modifications to each array and place the output into two separate bash arrays. Finally, I want to output the two separate bash arrays in the output section of the wdl task to be used in other tasks within the workflow.  The below code is an example of what I would like to do:
task AddA {
    input {
        Array[String] wdllistX = [ "True", "False"]
        Array[String] wdllistY = [ "Good", "Food"]
    }

    command <<<
        ###We add the letter 'a' at the end of each string in each wdl list and finally get the
      bash arrays:

        bashX = ("Truea" Falsea")
        bashY = ("Gooda" "Fooda")
    >>>
    output {
        Array[String] finX = bashX
        Array[String] finY = bashY
    }
}

The above example produces an error. It seems that I might have to convert the bash arrays back into wdl arrays in order to declare them as outputs in the output section, however, I am not sure how to do this. May someone help me with this? Moreover, a general explanation of how you can capture the output of specific things from the command section into the output section would be helpful aswell.

When a WDL task is run, the execution engine takes the command block and evaluates it as a template, rendering any wdl expressions it contains into the final script that is executed, BEFORE it is actually executed.
The actual execution environment is generally isolated from the execution engine or from other tasks. Information is passed between tasks by reading from or passing around file based outputs and not from accessing the bash environment directly, which the engine does not have access to.
So, in order to achieve what you are trying to do, you first need to get the values into a file format, and then you need to read those values using an engine function provided by WDL. Generally, all outputs need to be in a file if you want the engine to be able to interact with them. Please note that stdout and stderr are considered "files" and are accessible with stdout() and stderr() engine functions.
Option 1: Construct json files and use read_json
task AddA {
    input {
        Array[String] wdllistX = [ "True", "False"]
        Array[String] wdllistY = [ "Good", "Food"]
    }

    command <<<
        ###We add the letter 'a' at the end of each string in each wdl list and finally get the
      bash arrays:
        echo '["Truea","Falsea"]' > x.json
        echo '["Gooda","Fooda"]' > y.json
    >>>
    output {
        Array[String] finX = read_json("x.json")
        Array[String] finY = read_json("y.json")
    }
}

Option 2: use read_lines
task AddA {
    input {
        Array[String] wdllistX = [ "True", "False"]
        Array[String] wdllistY = [ "Good", "Food"]
    }

    command <<<
        ###We add the letter 'a' at the end of each string in each wdl list and finally get the
      bash arrays:

        echo "Truea" >> x.txt
        echo "Falsea" >> x.txt
        echo "Gooda" >> y.txt
        echo "Fooda" >> y.txt
    >>>
    output {
        Array[String] finX = read_lines("x.txt")
        Array[String] finY = read_lines("y.txt")
    }
}


text:
How to configure Cromwell to work with slurm?
I have a cluster that is set up with slurm and I seek to use Cromwell to dispatch jobs to the slurm scheduler. I followed the description explained on Cromwells Documentation Website. And created a configuration file called cromwellconfig.conf, however when I run it using the command java -Dconfig.file=/lccs/home/logio9/Configure_Cromwell/cromwellconfig.conf \ -jar cromwell-77.jar run test-workflow.wdl I get an error parsing the config file. I am confused as to where the error is, bellow is the contents of the configuration file:
include required(classpath("application"))

backend {
  default = SLURM
  providers {
    SLURM {
      actor-factory = "cromwell.backend.impl.sfs.config.ConfigBackendLifecycleActorFactory"
      config {
        runtime-attributes = """
        Int runtime_minutes = 600
        Int cpus = 2
        Int requested_memory_mb_per_core = 8000
        String queue = "short"
        String jobname = "testwork"
        String cwd = "/lccs/home/logio9/Configure_Cromwell"
        String err = "errors.output"
        String out = "ouput.out"
        """

        submit = """
            sbatch -J ${jobname} -D ${cwd} -o ${out} -e ${err} -t ${runtime_minutes} -p ${queue} \
            ${"-c " + cpus} \
        --mem-per-cpu ${requested_memory_mb_per_core} \
            --wrap "/bin/bash ${script}"
    """
    kill = "scancel ${job_id}"
        check-alive = "squeue -j ${job_id}"
        job-id-regex = "Submitted batch job (\\d+).*"
      }
    }
  }
}

This is the test wdl script I made, it copies a very large directory into another directory.
version 1.0

workflow Greetings {
    call SayHello {
    }
}

task SayHello {
    command <<<
        cp -R /lccs/home/logio9/Configure_Cromwell/tf1 /lccs/home/logio9/Configure_Cromwell/tf2

    >>>
    output {
        Array[String] out = read_lines(stdout())
    }
}

If anyone has had experience with this, please do share a bit on how exactly this configuration works.

Could you provide what is exactly the error that you are getting?

It says: java.lang.RuntimeException: Error parsing generated wdl: and
Caused by: wdl.draft2.parser.WdlParser$SyntaxError: ERROR: Sibling nodes have conflicting names:
Declaration defined here (line 5, col 1):
String cwd
^

Cromwell already defines the cwd, err and out runtime attributes by itself (it also defines a job_name and script). So when you define these variables in your config it causes a name conflict, since they already exist.
Cromwell will make a directory inside of the "cromwell-executions" directory (created when you run cromwell) for each job it runs. It will write the script, various stdout and stderr files,  and the rc file to this directory. The intend is also for the job to run in that directory. Hence cromwell provides these runtime attributes to point to the appropriate locations within this directory so they can be used in your custom submission command.
You can try just removing these variables in the runtime-attributes section, while keeping them in the submit section.

Thank you for your explanation, your solution works. The only issue that I am having is when I run Cromwell with the command above, Cromwell still seems to be running. If I do a CTRL + C, it stops Cromwell and quits the slurm job. Is there a way to simply dispatch the job to slurm and have Cromwell return/exit? Or perhaps have Cromwell run in the background?

@LogiosJames Cromwell is designed to run multiple consecutive jobs which rely on each other's outputs. So, it needs to keep running in order to check up on the status of the jobs it has dispatched and then submit the next jobs when their inputs have become available. You could just submit the cromwell process to the cluster in the traditional way, though, and have it run on an execution node (eg. sbatch [...] --wrap cromwell run worklow.wdl). Or use something like byobu/tmux and run cromwell in one window and do your other work in another.

Ok, good to note. I normally access and work with Cromwell on my cluster via an ssh connection, if I am running my workflow as I described above and my ssh connection to my cluster closes, would Cromwell continue executing, or would it exit/stop upon the connection closing?

@LogiosJames If you just run cromwell directly in your ssh session it will stop when you disconnect. If you submit cromwell to an execution node, it will keep running on that node. If you open a tmux or byobu session, those sessions will stay alive when you disconnect and you can reconnect with those sessions later using tmux attach-session -t 0 or byobu respectively.


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

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

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

import hail as hl

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

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

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


text:
Locate the saved template of minknow
I am using minknow on ubuntu. After saving the template (settings), where is that file locally stored in the system?
Also what is its format/extension?

It would help to understand a little better what you are actually doing. I can find a github page for MinKNOW and a page with a ton of tutorials. Where have you looked to try to find this information? Have you consulted the reference manual?

Hi! Thanks for responding. I want to share those files across a few systems, automate and load the exact same settings.
There is no content regarding the local storage for the templates.


text:
Is there a generally accepted method used to impute missing DNA methylation data (probe beta-values)?
I'm looking at DNA methylation (DNAm) data such as TCGA (e.g., BRCA, KIRP, KIRC, etc.). Currently trying to build use my model to predict DNAm age on test sets, but many of the data sets  are missing key probe values used in my clock/model.
I have inspected approximately 30 GEO data series and the following TCGA cancer samples: KIRP, KIRC, LUAD, LUSC, BRCA, THCA. They are each missing key probes I am using.
Is there a common way to impute the missing values in R?

Could you please clarify what "my clock/model" means, and what specifically is missing? How do you know that probes are missing? If this is from methylation chip data, have you confirmed that the probes exist in the probes used for that chip?

Hi, could you show what kind of data do you have and how you identify these missing probes? This might help having some code to build upon and answer your question. I think there isn't any general way to impute missing values, so knowing the data and its sparsity might help proposing imputation methods.


text:
DimPlot - How to highlight cells with identity colors?
I made this wrapper for DimPlot:
highlight_gene_expression <- function(seurat, gene_counts, colors = NULL, cells = NULL){
  Idents(seurat) <- 'cell_type'
  p <- DimPlot(
    seurat,
    raster = FALSE,
    shuffle = TRUE,
    reduction = "umap",
    pt.size = 0.01,
    cells = cells,
    cols = colors,
  )
  p <- p + guides(color = guide_legend(override.aes = list(size = 4), ncol = 1))

  if (!is.null(cells)) { gene_counts <- gene_counts[cells,] }
  expression_level <- {gene_counts |> rowSums()} / {gene_counts |> rowSums() |> max()}
  p[[1]]$layers[[1]]$aes_params$alpha <- expression_level
  p
}

It takes a Seurat object and a dataframe of gene expression levels.  It returns a UMAP with the transparency (alpha) of each point determined by the gene expression level:
highlight_gene_expression(
  seurat, # a seurat object
  trgd_counts, # A dataframe of gene expression levels.  Not important to understand for this question.
  cells = trgd_cells, # A vector of cell barcodes
  colors = colors # c("Erythrocytes"="#556b2f", "HSPCs"="#483d8b", "Proliferating T Cells"="#3cb371", ...)
)

Taking it one step further, I highlight the desired cells instead of just coloring them:
highlight_gene_expression <- function(seurat, gene_counts, cells = rownames(seurat@meta.data), colors = NULL){
  Idents(seurat) <- 'cell_type'
  p <- DimPlot(
    seurat,
    raster = FALSE,
    shuffle = TRUE,
    reduction = "umap",
    cols.highlight = colors,
    pt.size = 0.01,
    sizes.highlight = 0.01,
    cells.highlight = seurat@meta.data[cells,] %>% split(.$cell_type) |> lapply(rownames),
  )
  p <- p + guides(color = guide_legend(override.aes = list(size = 4), ncol = 1))

  expression_level <- {gene_counts |> rowSums()} / {gene_counts |> rowSums() |> max()}
  p[[1]]$layers[[1]]$aes_params$alpha <- expression_level
  p
}

This produces a similar plot, but with unhighlighted points in gray.  I think it looks quite nice:
# Same params as before
highlight_gene_expression(seurat, trgd_counts, cells = trgd_cells, colors = colors)

The problem is that the colors are in the wrong places.  Dimplot accepts my color vector, but it applies it sequentially when used for highlighting.  This means colors are no longer tied to a specific cluster as they were in the previous example.  A cluster will have a different color now depending on which ones are actually highlighted.
This is not a bug.  Rather, it is expected behavior if you look at the documentation for DimPlot.  The cols param corresponds to an identity class.  cols.highlight does not.  But I want it to.  How can I make cols.highlight explicitly define colors for each cluster?


text:
How to manage memory contraints when analyzing a large number of gene count matrices? I keep running out of RAM with my current pipeline
I have several hundred scRNA-seq count matrices, each from a different sample. For my other dataset containg a few dozen samples, I simply merged everything together into one Seurat object, but that won’t work here as far as I can tell. When I try to merge them in the typical way, I get the error:
Error in .cbind2Csp(x, y) :
  Cholmod error 'problem too large' at file ../Core/cholmod_sparse.c, line 89

This basically means I ran out of RAM.  I even ran it on the cluster with 1 TB of RAM, and I got the same error.  This is unsurprising, as memory requirements increase exponentially in relation to matrix size.  It seems that this happens when the merge() function attempts to convert all the sparse matrices to dense matrices.  Here is my code:
seurat_object_list <- "path/to/matrices" |> list.dirs(recursive = FALSE) |> lapply(\(matrix_folder) {
    matrix_folder |> Seurat::Read10X() |> Seurat::CreateSeuratObject()
})
seurat <- merge(seurat_object_list[[1]], seurat_object_list[2 : length(seurat_object_list)])

The Seurat object list is created with no issue, but it always fails on the merge.  I tested this code with a smaller subset of the matrices, and it works just fine.  It seems that I just have too many cells to process in this manner.
At this point, It looks like I have a few options.

I could hack together some way to force the matrices into one Seurat object.  But I don't know if that would even be useful, since I would likely keep bumping my head on the same memory contraints whenever I try to do anything with the combined object, such as running SCTransform() or FindMarkers().
I could divide the samples into groups that are small enough to merge, and process them in the normal way.  But this imposes limitations on the types of comparisons/analyses I could perform.
I could keep the matrices separate and analyze them iteratively.  This would solve the memory problems, but I haven't found any info online about anyone analyzing count matrices this way.  I imagine that I would lose a lot of out-of-the-box functionality that comes with the standard single-cell R libraries.  For example, I'm not sure how I would normalize all the gene counts against each other as separate objects.  This option seems like it might require a lot of re-inventing the wheel and a lot of code.

I like option 3 the best, but I'm not sure where to begin.  I am still pretty new to bioninformatics, so I would love to hear some input from those who are more experienced.  How do you deal with memory constraints when dealing with a large number of cells/samples?  Are there any standard ways to deal with this problem?  Is there some tutorial/doc online that I failed to find?  Any advice you could offer is greatly appreciated.

I get the same error but not running out of memory. I don't think this is a memory problem. It's some internal limitation of dgCMatrix object.

Scarf has been developed to specifically solve this issue. Scarf handles  large-scale scRNA-Seq/scATAC-Seq/CITE-Seq datasets that will otherwise not fit into memory. Unfortunately, this is not available in R yet.
Here is the Github link:
https://github.com/parashardhapola/scarf
Documentation:
https://scarf.readthedocs.io/
Paper:
https://www.nature.com/articles/s41467-022-32097-3
Please feel free to create an issue on Github for you use case (if you cant already find a solution in the documentation).
/PD


text:
How to convert 10x matrix.mtx.gz files to hdf5 format?
I have a bunch of folders containing barcodes.tsv.gz, features.tsv.gz, and matrix.mtx.gz from the Cellranger Count output for a single-cell dataset that was sent to me from another lab. I need the .h5 version of these count matrices, but I do not have access to the original Cellranger Count output which contains the generated .h5 files.
Is there an R or Python function that can reliably convert these files into the proper hdf5 format? I haven't found anything online yet.
(x-post: https://www.reddit.com/r/bioinformatics/comments/veku5r/how_to_convert_10x_matrix_tsv_files_to_hdf5_format/)

Conversion into HDF5 from Python3 is straight forward. Could you please forward your CellRanger output example? The key is ... the key and I'd need to see the stuff you want to write and look at the keys in a CellRanger file (keys are everything in HDF5)

Answered on biostars (https://www.biostars.org/p/9527636/#9527638) as well as on reddit: (https://www.reddit.com/r/bioinformatics/comments/veku5r/how_to_convert_10x_matrix_tsv_files_to_hdf5_format/icqsfwn/)
The biostars answer is better as it needs just one package (DropletUtils) as opposed to the 2-3 required by the Seurat-allied packages suggested on reddit. Plus, Seurat can take forever to load a dataset.

Conversion into HDF5 from Python3 is straight forward, if this is a simple key: value dictionary format then it should be easy to mimic the CellRanger format. If it is a more complex data structure that would be tricky.  The key to HDF5 is 'the HDF5 key' which you need to find (you would call it the directory)
To find the 'key' one approach is do is download https://www.hdfgroup.org/downloads/hdfview/ and then forward a screenshot of the CellRanger HDF5 file here. The directories are what HDF5 call 'keys' and you need these to write to and call from HDF5. There is a CLI for this operation.
An example of the data you are writing to would be needed.If per chance the data is in a Python dataframe you're in business, because these easily IO to HDF5
All you do is put the data into a simple dictionary and assign the dictionary to a key which is written to HDF5. Parsing the data from e.g. an AlignIO/SeqIO object to the dictionary is more complicated. When the HDF5 is read, the key has to be supplied then the dictionary is extracted and assembled back into an object.


text:
open source tools for CNV analysis from illumina SNP array
What are the best open source tools to analyse copy number variation using Illumina SNP array data...preferably R based tools, but open to other languages as well.
I'm working with leukemia sample and I have Illumina scan data ..in IDAT/GTC file format and i wish to use tools that will utilise the BAF and LRR metrics.
Microarray platform used- Illumina infinium 850K cytoSNP (beadchip array)

What kind of analysis? What questions are you trying to answer? Should we assume human?

@terdon I am trying to do CNV analysis...to detect copy number gains/loss and copy neutral loss of heterozygosity.... I'm working on acute lymphoblastic leukemia samples

Yes, I mean what kind of analysis? The word "analysis" suggests you aren't just trying to detect CNVs, but need to analyse them. Maybe annotate them? Classify them according to the ACMG criteria? If you just want to detect them, that's a different issue.

Maybe you can try karyoploteR. Great tutorial here: https://bernatgel.github.io/karyoploter_tutorial/


text:
GRanges coverage error with discordant reads
I mapped my fastq files from 4SU-seq experiment to human genome hg19 with bowtie 2 and got this results:
    67716784 reads; of these:
  67716784 (100.00%) were paired; of these:
    29941143 (44.22%) aligned concordantly 0 times
    33244334 (49.09%) aligned concordantly exactly 1 time
    4531307 (6.69%) aligned concordantly >1 times
    ----
    29941143 pairs aligned concordantly 0 times; of these:
      9287621 (31.02%) aligned discordantly 1 time
    ----
    20653522 pairs aligned 0 times concordantly or discordantly; of these:
      41307044 mates make up the pairs; of these:
        32859156 (79.55%) aligned 0 times
        3142531 (7.61%) aligned exactly 1 time
        5305357 (12.84%) aligned >1 times
75.74% overall alignment rate

I also ran a second mapping with --no-discordant and --no-mixed. After that I read the bam files into R like this:
x1 <- BamFile(paste0(bam.path, "14751.OnlyConcordantReads.sorted.bam"), asMates = T)
x2 <- BamFile(paste0(bam.path, "14751.AllReads.sorted.bam"), asMates = T)
x1 <- readGAlignmentPairs(x1)
x2 <- readGAlignmentPairs(x2)
x1_gr <- granges(x1, use.names=T)
x2_gr <- granges(x2,use.names=T, on.discordant.seqnames="drop")

The parameter on.discordant.seqnames is necessary, otherwise an error occurs.
z <- GRanges("chr1",IRanges(236958581,237067281))
xcov <- coverage(x1_gr)
xnum <- as.numeric(xcov$chr1[ranges(z)])
plot(xnum)
ycov <- coverage(x2_gr)
ynum <- as.numeric(ycov$chr1[ranges(z)])
plot(ynum)

I get very unusual coverage values (all summed up to 154 - see image) when I use on.discordant.seqnames="drop". Why is that?

In addition, the coverage does not match my bed files.


text:
System specifications for NGS data analysis
I have a 3.7 GB whole genome data of a eukaryote, for which genome assembly, gene prediction, and annotation steps have to be performed. In some time I would also need to analyze the transcriptome data of the organism. Can someone please let me know what are the system requirements for carrying out the aforementioned NGS data analysis? Thank you!

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

"3.7 GB whole genome data" what kind of data. sequencing? what platform? what tech specs?


text:
Extract single frame with MDAnalysis
How can I do to extract a single frame as a pdb file from a trajectory traj.xtc file, by using MDAnalysis?
I've got a universe
u=mda.Universe('confout.gro','traj.xtc')
Then I selected a residue by doing select_atoms, and I would like to extract a single frame as a pdb file

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

Q to more involved folk in admin: why was Community bot was invoked for this, as it was a case of nobody bothering to answer because it is a boring 'please link me to command X in the documentation' request. I.e. https://docs.mdanalysis.org/stable/documentation_pages/coordinates/init.html#writers

Check out gmx_trjconv - this command is in GROMACS and accepts several options where you can specify the specific single frame that you want to output/write - out to the new pdb/gro/xtc file.  https://manual.gromacs.org/current/onlinehelp/gmx-trjconv.html


text:
how to retain reads with low mapping quality (MAPQ) scores when using samtools view -q
I have been using the -q option of samtools view to filter out reads whose mapping quality (MAPQ) scores are below a given threshold when mapping reads to a reference assembly with either bwa mem or minimap2. However, using samtools idxstats to count total mapped reads and unmapped reads indicates that these reads with lower MAPQ scores are actually being discarded altogether, since the number of unmapped reads reported (last row) is always zero when the BAM file is created with a -q argument, even though fewer reads are being mapped to the assembly.
For example:
minimap2 -ax map-ont assembly.fasta reads.fastq | samtools view -b | samtools sort --reference assembly.fasta > nanopore_mapped_no_q_thresh.bam
samtools idxstats nanopore_mapped_no_q_thresh.bam

gives
contig_1      4751250 9897    0
contig_2      37697   528     0
contig_3      9131    2366    0
contig_4      98997   701     0
contig_5      4539    365     0
contig_6      3933    3164    0
*       0       0       142

while
minimap2 -ax map-ont assembly.fasta reads.fastq | samtools view -b -q 2 | samtools sort --reference assembly.fasta > nanopore_mapped_q2.bam
samtools idxstats nanopore_mapped_q2.bam

gives
contig_1      4751250 9716    0
contig_2      37697   528     0
contig_3      9131    1415    0
contig_4      98997   700     0
contig_5      4539    365     0
contig_6      3933    3161    0
*       0       0       0

Is there an alternative way to filter mapped reads by their MAPQ scores so that the filtered reads are still retained in the BAM file?

Looking at the samtools view docs, I believe that you need to also set the -p flag if you want to retain "filtered" reads in the file and simply have them be unmapped:

The -p, option sets the UNMAP flag on filtered alignments then writes them to the output file.

Yes, this does the trick. This flag was not recognized by the older version of samtools I was using (1.11) but after upgrading to 1.15 it works. Thank you Maximilian!


text:
How to calculate distance matrix from fasta?
I’d like to perform the following, but I am having difficulty.
I want to convert amino acid sequences in multifasta files into matrices (or vectors) with Euclidean distance and Pearson correlation coefficient.
If you know a solution, please let me know. Thanking you in advance.
File1
>BDP27474.1 |surface glycoprotein [Severe acute respiratory syndrome coronavirus 2]
MFVFLVLLPL
>BDP27485.1 |surface glycoprotein [Severe acute respiratory syndrome coronavirus 2]
MFVFLVLLPL

File2
>BDD77411.1 |surface glycoprotein [Severe acute respiratory syndrome coronavirus 2]
MFVFLVLLPL
>BDD77515.1 |surface glycoprotein [Severe acute respiratory syndrome coronavirus 2]
MFVFFVLLPL

Code
library(tidyverse)
library(ape)
library(Biostrings)

naseq1 <- readAAStringSet("File1.fasta")
naseq2 <- readAAStringSet("File2.fasta")

dist.aa(naseq1)
dist.aa(naseq2)

##Error in numeric(n * (n - 1)/2) : invalid 'length' argument

Add info.
iOS 10.15, Rstudio cloud

There's no matrix for file1 because MFVFLVLLPL and MFVFLVLLPL are identical

You are probably looking for Biostrings::stringDist to compare both strings


text:
How to assess quality of WGCNA module identification in blockwisemodules()
I'm exploring WGCNA for bulkRNA sequencing analysis with human subjects.  I have healthy, myocarditis, and heart failure patients (which can be further broken into ischemic or nonischemic)
I have been spending quite a bit of time tinkering with the parameters in blockwisemodules() to identify modules/MEs, and I'm only able to assess the quality of this step with a rudimentary intuition.
For example, if I use
net = blockwiseModules(datExpr, power = 10,maxBlockSize = 20000,TOMType = "unsigned", minModuleSize = 20,reassignThreshold = 50, mergeCutHeight = 0.35,numericLabels = TRUE, pamStage = TRUE, pamRespectsDendro = TRUE)
I get what is so far the "best" module detection - many modules (although small), and about 60-70% of genes seem to fall into a module. Image below.
I find myself at a bit of a loss as to how to judge the quality of this step.  The module detection looks extremely messy compared with the tutorial data, which i wasn't surprised at really, but can anyone provide guidance here?
I can provide other images of what I get from this step that looks like much poorer detection of modules if it would help.
Last related question, I have 19 healthy, and about 41 patient samples, I wasn't sure if the best approach to module detection is to put all samples together, then look at what modules are associated with disease states which i could convert into numeric for that step, or would it be better as i have started - to detect modules in healthy, or patient groups in independent datasets?


text:
Which components must be included in the HDR plasmid?
I am trying to design a plasmid for our knock in experiment and wanted to clarify some technical issues.
I wanted to know exactly which components must be included in the plasmid.
I know that I have to include the 5' homology arm , insertion gene sequence, marker sequence and 3' homology arm.
But, I would like to know for certain the following details:

Should I take the whole sequence of the gene or do I need to remove the introns/UTR regions or any other unnecessary part of the sequence?
Which sequence should I choose for gene insertion , gene sequence or mRNA transcript sequence?
When designing plasmids, should I insert extra promoters and/or terminators with the insertion genes?

Additionally, can you recommend an online tool to design HDR Plasmid?
I have already tried using Alt R HDR design tool and Benchling but our target genome is not available there.

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

You would probably get better answers on the Biology forum. It would also help to mention whether you are using a commercial backbone for the plasmid or designing the whole construct from scratch.

Thank you @Cloudberry for your suggesion.


text:
How are BLAST record scores calculated in biopython
Can someone explain how the score is calculated for Bio.Blast.Record.Description?  For example, if Record.Description.num_alignments is >1, is Record.Description.score simply the sum of the alignment scores for the description's hit?  Can alignments overlap?  How does the score of a given alignment deal with the scores of overlapping HSPs in that alignment?
I'm getting the Blast.Record object via NCBIXML.parse(), and would like to rely on the description data as much as possible to remove the need for code that parses the Record.Alignments objects, but I'm a nervous including things when I don't know exactly how they were created.
Cheers!
Tim

The results you'd get are the same as with any other NCBI BLAST+ tools, python does not compute it.
Biopython can run a blast query qblast on the internet or locally. To make it run locally you need to install blast tools on your machine and biopython will launch it as a subprocess. Here you can find the pertinent information: https://biopython.org/docs/dev/Tutorial/chapter_blast.html#running-blast-locally
Filtering data is as much a burden as it is important. Depending on your question you want to answer. Some people prefer a low e-value. The score is however affected by the length of the sequence since the number of matches is bigger, more details here: https://www.nlm.nih.gov/ncbi/workshops/2023-08_BLAST_evol/blast_score.html


text:
Polygenic Risk score calculation using vcf files
I have about 750 vcf files and I am looking to generate a Polygenic Risk Score for each individual. I am planning to use either PLINK or LDPhred to generate a Risk Score. I have phenotype information of these individuals.
Can anyone verify if the polygenic risk score is calculated using individual vcf files or should I look into combining the vcf files?

What kind of polygenic risk score? For what diseases? What genes? Using which of the many approaches?

For the tools you mentioned, you should convert to BED format. Other tools will accept VCF. It should work either way, combined or not, though you may have to adjust your pipeline depending on which you prefer. I'd imagine a combined file would be easier to work with.


text:
Bar Graph of Expression Data from Seurat Object
I am looking for a way to present some single-cell RNA sequencing data on expression of a certain gene. My dataset has 2 variables, cell type and condition. I'm looking to create a grouped aligned scatter plot bar graph of the expression of the gene, where the y-axis is expression levels, while X axis is the different cell types, with each cell type having 2 bars (1 for healthy, 1 for diseased). My dataset has 3 healthy and 3 diseased samples, but all of the data is integrated into a Seurat object. To first create an aligned scatter plot bar graph, what I did was generate a DotPlot for the expression of gene X in each sample, split by cell-type. Then I realized that DotPlot() prints out a data table with values including avg.exp, pct.exp and avg.scaled.exp. I used the avg.exp metric to plot the data, but I was posting this question to see if there are other ways I can access and present expression data from a Seurat object. Is the way I've presented accurate? I've posted a sample bar graph of what I'm looking for, and the main question is how I can gather the expression levels across cell types AND diseased or healthy condition.
Thanks!

Regarding the other answer with the below code:
data.df <- object@meta.data[, c("celltypes", "condition")]
data.df["value"] <- object$RNA@data["geneX", ]
data.df <- data.df %>% mutate(group = paste(celltypes, condition, sep=":"))
data.df <- data.df %>% filter(value > 0)
data.df.2 <- data.df %>% group_by(group) %>% summarise(sd=sd(value),value=mean(value))
X <- str_split(data.df.2$'group', ":", 2, simplify=TRUE)

data.df.2$'celltypes' <- X[, 1]
data.df.2$'condition' <- X[, 2]

ggplot(data.df.2, aes(fill=condition, y=value, x=celltypes)) +
    geom_bar(color="black", position="dodge", stat="identity") +
    geom_errorbar(
        aes(ymin=value, ymax=value+sd), width=.2,
        position=position_dodge(.9)
    )

I find the plot matches more the results of Seurat differential expression if you DO NOT perform the step of removing the cells with 0 counts (data.df <- data.df %>% filter(value > 0)). To avoid having huge error bars you can use standard error of the mean instead of sd.

In single-cell data a bargraph with error bar is not usually shown since the gene expression does not have a normal distribution due to drop-outs. A violin-plot is commonly used. In your case, first set the active.ident to cell types
seuratObj <- SetIdent(seuratObj , value = "Cell.Types")

Then call the violin plot
VlnPlot(seuratObj , features = "geneX", split.by = "condition")

If you still specifically need bargraph, you can get the output of the violin plot as a dataframe and then draw the barplot.

Following the answer of @Ajay.
You can add + geom_boxplot() after the command to turn the violin plot to a box plot.
As for the barplot, you can do it with these:
# the object is your seurat object
data.df <- object@meta.data[, c("celltypes", "condition")]
data.df["value"] <- object$RNA@data["geneX", ]
data.df <- data.df %>% mutate(group = paste(celltypes, condition, sep=":"))

# remove cells that do not express the gene
data.df <- data.df %>% filter(value > 0)

data.df.2 <- data.df %>% group_by(group) %>% summarise(sd=sd(value),value=mean(value))
X <- str_split(data.df.2$'group', ":", 2, simplify=TRUE)
data.df.2$'celltypes' <- X[, 1]
data.df.2$'condition' <- X[, 2]

ggplot(data.df.2, aes(fill=condition, y=value, x=celltypes)) +
    geom_bar(color="black", position="dodge", stat="identity") +
    geom_errorbar(
        aes(ymin=value, ymax=value+sd), width=.2,
        position=position_dodge(.9)
    )

Dropout is very usual in single-cell data. If you really want to show the expression with barplot, I would suggest removing the zero values or using denoised data for better visualisation.


text:
Picard's CollectGcBiasMetrics gc-bias plot- help with interpretation
Setup
I ran picard's CollectGcBiasMetrics tool for my WES sample, with the intention of assessing the performance of the Exome capture panel in terms of the uniformity of read coverage at various %GC windows. I obtained the gc bias plot attached below.

Question
I am a bit puzzled by the normalised coverage values in the windows containing higher GC % and need help with interpretation - in particular, whether it indicates that something has gone awry with the sample, or if it is merely the way the normalised coverage is computed in picard. Why do the normalised coverage values seem to peak at around 50-55 GC%?
Digging into the gc_bias_metrics.txt file:

The #windows, #read starts and normalised coverage values are highlighted in the red, blue and yellow boxes respectively. Here, I notice that the # reads / # windows ratio is high for the tail end GC content, which probably inflates the normalised coverage value.
Is this something unexpected for a WES sample, i.e. should I be concerned?
Follow-up question
If the coverage values provided by picard isn't that useful, are there other useful tools for me to calculate the uniformity of read coverage across various GC%? I am looking at plots like Figure 3 in the twist biosciences data page.

Often coverage is higher in the relatively high GC range due to some artefactual biases from things like PCR, though the exome kits are getting better these days (and can be PCR-free): https://www.nature.com/articles/s41598-020-59026-y


text:
Which type of variant caller should I use in a WES normal cell line sample?
I have whole-exome sequencing data of an immortalised non-tumor (normal) cell line that I wish to assess for the presence/absence of APC/Wnt mutations. This is to double check that the cell line is sufficiently "clean", without mutations in key regions, for further experimental analysis.
I understand that without a matched normal sample, I won't be able to distinguish if a variant is a somatic variant or a rare inherited germline mutation. As I just want to check for mutations in those regions, it is not important to distinguish between these types.
Question: Should I use the GATK haplotypecaller or a somatic variant caller such as Mutect2 without a matched-normal sample?

I don't understand. If this is not a tumor sample, why would you want to run a somatic variant caller? Do you want to run this in a super sensitive way to detect everything that might conceivably be a variant?


text:
Transcriptome analysis
I am trying to assemble reads belonging to two different readlength. Is it a valid way since I am looking for common genes among the species I am assembling.

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

My work focuses on identifying common genes amongst different species of insects. I am doing transcriptome analysis to find these genes. However, I am assembling reads from different readlength using Trinity. I am unable to understand how to proceed next.


text:
RNASeq NMF clustering using what kind of expression unit?
I would need some help regarding certain points in a RNASeq gene expression analysis (sorry if there are stupid questions, its my first project and i have no capable supervision). My workflow is the following:

TMM normalization

obtaining CPMs

filtering CPMs that are under 1 in half the samples

obtaining log2(CPM) - after I perform this step I get Inf values, should I treat those values as NA and remove them or should I use CPM + 1 (I know that adding +1 introduces certain bias)?

apply Median Absolute Deviation and use different thresholds to subset data

apply NMF clustering and find cluster specific genes and samples

use the sample clusters to obtain DEGs (using DESeq2) and proceed with GSEA

Now for the questions.. Firstly, I tried to look for the best practices in scientific papers and on forums but they vary in a huge manner.
In papers, such as https://pubmed.ncbi.nlm.nih.gov/26909576/ the authors use log2(CPM) for NMF clustering analysis, while looking online at other suggestions, they all suggest using variance stabilizing transformation or z-scores etc. Now, what I'm trying to do is find subtype specific genes and significant DE genes based on those subtypes.
From what I understand, VST transformation is ideal but why do all the papers that I've found only use TPM/CPM for subsequent analysis such as NMF clustering (other papers only mention gene expression or I missed that)?
I've tried both ways, with log2(CPM) and with VST and I've found that VST clusters contain more genes and tend to be of equal size (both genes and samples clusters) while log2(CPM) clusters returns a cluster 3 times bigger than the second (both for genes and samples). Both VST and log2(CPM) clusters were selected by the best cophenetic coef and also samples had silhouette widths above 0.6.
Now, what should be the best approach considering that I'm expected to write a decent paper with this dataset.
Thank you in advance!
PS: all the above are done in R


text:
Error of duplicated rownames although there are no duplicates
I have a data frame that I want to switch its row names from EnsembleID (GENEID) to gene symbol (SYMBOL). When I try to switch, I get this error:
Error in `.rowNamesDF<-`(x, value = value) :
  duplicate 'row.names' are not allowed
In addition: Warning message:
 Error in `.rowNamesDF<-`(x, value = value) :
  duplicate 'row.names' are not allowed

The code I'm using works for other similar data frames. Only this one is giving me trouble. There are no duplicates in the row names, I checked, so the source of this error is strange.
This is the code:
GetSymbols <- function(df) {   # Switch from geneID format to Symbol format

  ensembl.genes = as.vector(rownames(df))
  geneIDs = ensembldb::select(EnsDb.Hsapiens.v86, keys= ensembl.genes, keytype = "GENEID", columns = c("SYMBOL","GENEID"))
  dupli = duplicated(geneIDs$GENEID)
  geneIDs = geneIDs[!dupli,]
  rownames(geneIDs) = geneIDs$GENEID
  rownames(df) = geneIDs$SYMBOL
  return(df)

}

#Calling the function:
 df = GetSymbols(df)

A sample of the dataframe:
structure(list(logFC = c(0.362414024114299, 0.270660335317552,
0.250197714648158, 0.355106754020497, 0.860436247608723, 0.331838240496886,
0.642000001979217, -0.506238094622956, 0.199130048525963, 0.226010606161764,
-0.300273572714138, 0.763447672404514, -0.282220097838676, 0.686284440644756,
0.209266959343866, 0.201216572039244, 0.467207891400503, 0.207476905368996,
0.548263917432251, 0.582095719576479, 0.167369145307915, -0.239886539805618,
-0.133088162158502, 0.432528267703316, -0.232843401197013), AveExpr = c(2.63663927914016,
6.23730774047258, 6.74130527553022, 6.51256660562955, 3.51386518575641,
0.962682352871946, 3.60110717340049, -1.07537043066235, 4.55672701732791,
3.96125159652725, 1.20473881401599, 1.72955742088622, 2.10495435034315,
2.95446094011986, 4.31205976523475, 4.85412008885978, 4.09171246258306,
3.94681319329962, 6.27193781389662, 0.936821016422818, 7.97270280545117,
2.4470079966345, 5.75991260844846, 3.84225682530045, 3.35515997418382
), t = c(5.15993300980101, 4.74644046948475, 4.67895886939685,
4.65270476519835, 4.39369295809775, 4.51083389626156, 4.20489542862464,
-4.49987927160777, 4.1096966271885, 4.03273843222089, -4.09143780518073,
4.07713970281633, -3.99157698713235, 3.97814493972727, 3.94402861982903,
3.92673590933521, 3.92436825669156, 3.90538427413924, 3.89364498564491,
3.93360332510695, 3.80233390398449, -3.76935782400013, -3.7013905266149,
3.68313224298669, -3.68114307289328), P.Value = c(4.50582660244125e-07,
3.21245260466636e-06, 4.37260426507432e-06, 4.92523952290176e-06,
1.54869545509172e-05, 9.2842399715271e-06, 3.45373717288344e-05,
9.74370585569957e-06, 5.12014322967072e-05, 7.00192255628614e-05,
5.51721193906694e-05, 5.84844087520547e-05, 8.26184479488647e-05,
8.71763749953859e-05, 9.9850615705057e-05, 0.000106924308803713,
0.000107929018751166, 0.000116314737750083, 0.000121805424073131,
0.000104060063751203, 0.000173710084683786, 0.000197140244520663,
0.000255164374909145, 0.000273300166932535, 0.000275347732934516
), adj.P.Val = c(0.00727510763230165, 0.0198807293341929, 0.0198807293341929,
0.0198807293341929, 0.0357217668827299, 0.0262203124576875, 0.06970505049172,
0.0262203124576875, 0.0858444785191522, 0.09421086799483, 0.0858444785191522,
0.0858444785191522, 0.0968123298197961, 0.0968123298197961, 0.0968123298197961,
0.0968123298197961, 0.0968123298197961, 0.098333518854239, 0.098333518854239,
0.0968123298197961, 0.130401412989609, 0.138392451653505, 0.150726778618736,
0.150726778618736, 0.150726778618736), B = c(5.59649208873277,
4.23844568248881, 3.95926511142815, 3.85011503575047, 2.70095619489428,
2.52151007226572, 2.02920447371096, 1.82077681362955, 1.73104697407381,
1.44381690648743, 1.32264026808754, 1.26349133643967, 1.17476645994761,
1.17423988785926, 1.13844103561967, 1.07813935334313, 1.06670248721183,
1.0004861785984, 0.949862540556422, 0.634079467108277, 0.620590447127368,
0.489071423128816, 0.29054426889472, 0.254833988789798, 0.246174538031329
)), row.names = c("ENSG00000188938", "ENSG00000182768", "ENSG00000174579",
"ENSG00000119523", "ENSG00000271798", "ENSG00000204923", "ENSG00000281910",
"ENSG00000271092", "ENSG00000153767", "ENSG00000116903", "ENSG00000168101",
"ENSG00000206834", "ENSG00000180035", "ENSG00000251733", "ENSG00000058729",
"ENSG00000132467", "ENSG00000212464", "ENSG00000115561", "ENSG00000230669",
"ENSG00000206612", "ENSG00000275272", "ENSG00000162139", "ENSG00000143889",
"ENSG00000207233", "ENSG00000119878"), class = "data.frame")

Thanks very much! This is extremely helpful when doing immune deconvolution. I used quanTIseq and it said my row names are duplicated while all of my row names are unique. However, after removing all the ids with the duplicated() function, the immune deconvolution worked! May I know what magic duplicated() do?

I found an answer shortly after posting this question:
ensembl.genes <- rownames(df)
geneIDs <- ensembldb::select(EnsDb.Hsapiens.v86, keys= ensembl.genes, keytype = "GENEID", columns = c("SYMBOL","GENEID"))
dupli = duplicated(geneIDs$SYMBOL)
geneIDs = geneIDs[!dupli,]
ind = intersect(rownames(df), geneIDs$GENEID)
df = df[ind,]
rownames(df) = geneIDs$SYMBOL

Can you please mark this as an accepted answer if it works for you?


text:
How can to validate the presence of a certain type of cells in a single cell dataset?
I have a single cell dataset which I consider to be a reference , let say, for an human organ, I have identified some new clusters that correspond to cell types. I also have WGCNA modules that correspond to a identified cluster.
Given this background, is there any possibilities to validate the presence of the cell types from my reference dataset inside another single cell dataset ?
I was thinking about GSEA, but all the gene in a set have the same importance.
Could I use eigen vectors from first principal component for a given module expressed in the corresponding cell type, as a weight. Then calculate inside another dataset a sort of "first principal component" using the previous eigen vectors. Then I would simply calculate the correlation between the first principal component and the calculated one to validate whether my cell type is present or not.
I was also thinking about a sort of weigthed GSEA, if anything like this ever existed.

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

I think the typical way to do this would be to identify a small number of marker genes whose presence and absence is specific for that cell type.


text:
The difference between sequence-independent and sequence-order independent structural alignment?
I am confused about the difference between sequence-independent and sequence-order independent structural alignment. TM-align for example claims to be sequence-independent; however, I've seen it is also sequence order independent from the papers. To me, they sound like the same thing, but obviously, they aren't. Does anyone know about anything the difference?

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

TM-align is an algorithm for sequence independent protein structure comparisons. For two protein structures of unknown equivalence, TM-align first generates optimized residue-to-residue alignment based on structural similarity using heuristic dynamic programming iterations. An optimal superposition of the two structures built on the detected alignment, as well as the TM-score value which scales the structural similarity, will be returned.

TM-score has the value in (0,1], where 1 indicates a perfect match between two structures. Following strict statistics of structures in the PDB, scores below 0.2 correspond to randomly chosen unrelated proteins while those higher than 0.5 assume generally the same fold in SCOP/CATH.   https://seq2fun.dcmb.med.umich.edu//TM-align/ https://seq2fun.dcmb.med.umich.edu//TM-align/

maybe  sequence order independent from the papers.  means tha using with input A,B or B,A returns same result ??


text:
PL and QUAL values on VCF file?
I want to filter my VCF file to include only the relevant information I need but I have some questions about the results I got.
This is what the first 10 lines of my VCF file looks like:
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  DATA_VARIANT
chrI    826 .   T   C   255 .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    1/1:344,30,0:57:10:0,10,0,0:0,5
chrI    831 .   C   T   71  .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    0/1:139,36,274:71:12:0,8,0,4:3,2
chrI    836 .   G   T   59  .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    0/1:140,49,373:59:16:0,0,12,4:4,1
chrI    848 .   C   T   255 .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    1/1:556,48,0:75:16:0,0,0,16:0,5
chrI    855 .   G   A   255 .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    0/1:399,48,139:59:16:12,0,4,0:1,4
chrI    858 .   ATT AT  255 .   CNV=1;TYPE=INDEL    GT:PL:GQ:DP:ADP:ACN 0/1:480,48,160:85:16:4,12:1,4
chrI    876 .   T   C   255 .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    1/1:865,75,0:102:25:0,25,0,0:0,5
chrI    904 .   G   A   255 .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    1/1:1261,112,0:138:37:37,0,0,0:0,5
chrI    920 .   A   T   67  .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    0/1:244,145,1426:67:48:41,0,0,7:4,1
chrI    925 .   T   C   255 .   CNV=1   GT:PL:GQ:DP:BSDP:ACN    1/1:1786,157,1:255:52:0,52,0,0:0,5

QUAL:
As I understand it, the QUAL is based off of Phred quality scores, which would mean that a QUAL value of anything greater than 20 would be fairly sufficient in terms of quality for that genotype. However, most of my QUAL scores are 255 and from what I've read that means there is no mapping quality available?

What does this mean for my data? I used the same tool (NGSEP) to map and call the variants on this strain so it should have mapping quality scores.
Do I need to filter out all the 255 scores? Or can they still be used?

PL: From what I've read about PL scores, they indicate the likelihood of the presence of a genotype (0/0, 0/1, 1/1). Based on the results, there should be a clear least likely genotype, "almost" likely genotype, and a most likely genotype indicated by a value of 0.

Some of my PL values do not include a 0 (i.e., the second variant call = 139,36,274). What does this mean for my data?

Any help would be great! Thank you!
Note: I used NGSEP as my variant caller

The VCF specification provides the definition for the QUAL field.
However, QUAL values are often capped by variant callers to a given value. I have seen 100 being used and according to this documentation on NGSEP it seems 255 has been chosen for this caller as the maximum value. So you are all good keeping these high QUAL variant sites and filtering only those below a certain value (e.g. 20)
PL is also defined in the VCF specification: phred-scaled genotype likelihoods rounded to the closest integer (and otherwise defined precisely as
the GL field).
Again, different variant callers may have different implementations for the specification. GATK for example normalizes the values so that the most likely genotype has a PL of 0. But this might not be followed in other callers. It seems that NGSEP simply does not apply this normalization or does it differently. I think this is a detail you can ask directly the authors here by opening an "issue".
My suggestion is that you rather use GQ (genotype quality) if you want to hard-filter your variants. This is more standard to do, as opposed to filtering using PLs. These values are also in phred-scale.


text:
need bam file for pilon
I just ran an assembly on yeast genomes using Flye and I want to polish those assemblies with Pilon but it requires a sorted BAM file.
How do I make a BAM file of the resulting assembled.fasta?

What reads do you have? Only long read or short read or whatever? You need to align the starting reads to your assembly and then supply that bam, as implied by the user172818 answer.

Run a short-read mapper. For example:
bwa index assembled.fasta
bwa mem -pt16 assembled.fasta read1.fq.gz read2.fq.gz \
  | samtools sort -m4G -@4 -o align.bam -

This is good if the user has Illumina, that isn't clear (though it is implied by using pilon).