text:
replacement has 1 row, data has 0
I have tried to use UpsetR to visualize the input file which can be found here
> library("UpSetR")
> orthogroups_df<- read.table("orthogroups.GeneCount.tsv",  header=T, stringsAsFactors = F)
> #All species
> selected_species <- colnames(orthogroups_df)[2:(ncol(orthogroups_df) -1)]
> selected_species
 [1] "Atha" "Cann" "NQLD" "Natt" "Ngla" "Nlab" "Nsyl" "Ntab" "Ntom" "Slyc" "Stub" "Vvin"
> head(orthogroups_df)
  Orthogroup Atha Cann NQLD Natt Ngla Nlab Nsyl Ntab Ntom Slyc Stub Vvin Total
1  OG0000000    0    0  965    0    0    3    0    0    0    0    0    0   968
2  OG0000001    0    1    3    0    0  448    0    0    0    0    0    0   452
3  OG0000002    0    1  313    0    0  120    1    0    1    0    0    0   436
4  OG0000003    0   93   15   21   46   16   33   63   36   25   39   26   413
5  OG0000004    1   42    2   34  109    6    8  154   11    9    4    0   380
6  OG0000005    0    2   61    1   34   44   91   70   43   20    1    0   367
> ncol(orthogroups_df)
[1] 14
> orthogroups_df[orthogroups_df > 0] <- 1
> upset(orthogroups_df,
+       nsets = ncol(orthogroups_df),
+       sets = rev(c(selected_species)),
+       queries = list(list(query = intersects, params = list("NQLD", "Nlab", "Nsyl"), color = "#238c45", active = T),
+                 list(query = intersects, params = list("NQLD", "Nlab"), color = "#ffd977", active = T)))

Error in `$<-.data.frame`(`*tmp*`, "freq", value = 45L) :
  replacement has 1 row, data has 0

How is it possible to fix the above error?

Please provide some reproducible data first. Nobody is going to be able to help otherwise.

I suspect the error is caused by a 'lack' of intersection between "NQLD", "Nlab", and "Nsyl"? If you remove "Nlab" from the first query the command appears to complete as expected, e.g.
library(UpSetR)
orthogroups_df <- read.table("~/Desktop/orthogroups.GeneCount.tsv",
                            header=T, stringsAsFactors = F)
selected_species <- colnames(orthogroups_df)[2:(ncol(orthogroups_df) -1)]
selected_species
#>  [1] "Athaliana_447_Araport11.protein_primaryTranscriptOnly"
#>  [2] "Cann"
#>  [3] "NQLD"
#>  [4] "Natt"
#>  [5] "Ngla"
#>  [6] "Nlab"
#>  [7] "Nsyl"
#>  [8] "Ntab"
#>  [9] "Ntom"
#> [10] "Slyc"
#> [11] "Stub"
#> [12] "Vvin"
head(orthogroups_df)
#>   Orthogroup Athaliana_447_Araport11.protein_primaryTranscriptOnly Cann NQLD
#> 1  OG0000000                                                     0    0  965
#> 2  OG0000001                                                     0    1    3
#> 3  OG0000002                                                     0    1  313
#> 4  OG0000003                                                     0   93   15
#> 5  OG0000004                                                     1   42    2
#> 6  OG0000005                                                     0    2   61
#>   Natt Ngla Nlab Nsyl Ntab Ntom Slyc Stub Vvin Total
#> 1    0    0    3    0    0    0    0    0    0   968
#> 2    0    0  448    0    0    0    0    0    0   452
#> 3    0    0  120    1    0    1    0    0    0   436
#> 4   21   46   16   33   63   36   25   39   26   413
#> 5   34  109    6    8  154   11    9    4    0   380
#> 6    1   34   44   91   70   43   20    1    0   367
ncol(orthogroups_df)
#> [1] 14
orthogroups_df[orthogroups_df > 0] <- 1

upset(orthogroups_df,
      nsets = ncol(orthogroups_df),
      sets = rev(c(selected_species)),
      queries = list(list(query = intersects,
                          params = list("NQLD", "Nsyl"),
                          color = "#238c45", active = T),
                     list(query = intersects,
                          params = list("NQLD", "Nlab"),
                          color = "#ffd977", active = T)))

Created on 2022-11-28 with reprex v2.0.2
You can see in the plot that "NQLD", "Nlab", and "Nsyl" don't appear to have any values in common, i.e. the 'dots' corresponding to those variables aren't joined. So, to fix the error you need to select variables that intersect in your queries.


text:
Why do we need to find minimum energy in a protein chain?

High–quality protein backbone reconstruction from alpha carbons using Gaussian mixture models

The above research paper is about a software tool for reconstructing a protein's main chain model only from its alpha-carbon backbone (aka, C-alpha trace).
This tool provides two options: conversion (1) with energy minimization, and (2) without energy minimization.

Energy Minimization

The above paper says that:

...
The goal of energy Minimization is to find a set of coordinates representing the minimum energy conformation for the given structure.
...
In molecular modeling we tend to are particularly curious about minimum points on the energy surface. Minimum energy arrangements of the atoms correspond to stable states of the system; any movement off from a minimum provides a configuration with a better energy.
...

I haven't understood the quote from the second paper.
Why do we need to know the minimum energy conformation of a protein chain? I.e., why do we want to do energy minimization?

After translation, proteins fold to achieve a specific 3D structure which allows them to perform their specific functions. The structure of a protein computationally inferred through energy minimization often corresponds quite well to that specific, real-life folded 3D structure.
If you are interested in the exact reason why proteins tend to fold to a conformation of minimal energy (even though this is an oversimplification of what's really happening) you would have to dive a bit deeper into the thermodynamics of proteins: https://en.wikipedia.org/wiki/Protein_dynamics
TL;DR: Energy minimization is often used to get an idea of the protein's real-life conformation (i.e. structure) through computational means.

dont forget to add post-translational modification thought [i.e. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6779701/ ]

source: Stack Exchange

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

sort: votes


text:
How can I extract the Protein distance map data in Julia?
I calculated the distance between C-alphas using Julia,
but there is no way to extract distance map data from Julia.
I don't even know the data file form, png or csv or etc.
I don't want to extract it as an image, I just want all raw data.

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

Can you post an example of your input data, some code, any idea what your inputs and desired outputs are supposed to look like?

for someone who have same question as mine...
this helps you to extract distance map data file.
julia> using DelimitedFiles

julia> writedlm("calphadistances.out",dm.data,"")

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


text:
Bcftools 1.16 able to add F_MISSING tag?
I tried adding the F_MISSING tag using bcftools 1.16. When I run this command:
bcftools +fill-tags input.vcf.gz -- -t 'F_MISSING' | bcftools view -i 'INFO/F_MISSING<0.25' -Oz -o output.vcf.gz

I get the following error:
Error parsing "--tags F_MISSING": the tag "F_MISSING" is not supported

This command runs fine using bcftools 1.15. However, version 1.15 gives complications with other packages I use in my snakefile. Do you maybe know alternatives for how to add F_MISSING using bcftls 1.16 (if possible)?
When I type bcftools +fill-tags --version:
bcftools 1.9 using htslib 1.9
plugin at 1.9 using htslib 1.9
##SOLUTION##
Indeed the issue was that I was not installing the most recent version of Conda.
I solved it by changing the .condarc file to solelely include the following lines:
channels:
    - conda-forge
    - bioconda
    - defaults

The order is crucial as well.

Wait, what is version 1.9? Aren't you using version 1.16?

I know, I am pretty confused by that myself. I installed bcftools in a newly created conda env using conda install -c bioconda bcftools as indicated on https://anaconda.org/bioconda/bcftools and conda says v1.16

Conda might, but are you running the version installed in conda or the one on your system? Don't you need something like conda run bcftools --version?

I ran the following:
conda activate bcftools1.16  and then  conda list --show-channel-urls.
Just to check again the version of bcftools:
# Name                    Version                   Build  Channel
bcftools                  1.9                  h68d8f2e_9    bioconda

Please put all of this into the question. I don't know much about conda, but it still looks like you're running the system-wide one. What's the output of type -a bcftools and, type -a bcftools | grep -v aliased | grep -oP '/\S+' | while read c; do echo "=== $c ==="; $c --version; done?

I solved it, I added my solution to the question.

Great! But please post your solution as an answer instead so the question can be marked as slved.

Thank you for doing it for me. Very new to stackoverflow as you might have guessed.

The issue was that not installing the most recent version of Conda.
It was solved by changing the .condarc file to solely include the following lines:
channels:
    - conda-forge
    - bioconda
    - defaults

The order is crucial as well.
Converted to answer from question.


text:
Definition of genotype in demuxlet
I am reading the Online methods of the demuxlet paper.
The genotype $g$ is taken from the set $\{0, 1, 2\}$, defined as "the true genotype of the sample corresponding to $c$-th droplet at $v$-th variant".
What is the mapping between these values and the corresponding genotype?
Each variant is modeled using a binomial with $n=2,p=\frac{g}{2}$, so I would assume that 1 is heterozygous, 0 is homozygous REF and 2 homozygous ALT. Is this assumption correct?

The genotypes $\{0,1,2\}$ are mappings to the genotypes AA (homref), AB (heterozygous), BB (homalt) (from supplementary resources). This is also mapped to the number of alterations (e.g., 0 -> 0 alterations -> homref).


text:
file path of BiocManager:install()
I was trying to format my GWAS summary statistics with the MungeSumstats R package and would need to install SNPlocs.Hsapiens.dbSNP155.GRCh37 and BSgenome.Hsapiens.1000genomes.hs37d5 as reference panel.
I typed
BiocManager::install("SNPlocs.Hsapiens.dbSNP144.GRCh37") BiocManager::install("BSgenome.Hsapiens.1000genomes.hs37d5")
but failed due to Internet connection problems. So I planned to download the files from the web browser and copy them to the R library. I was wondering where I should put the files. If they were downloaded from R, where would they be saved?
I would really appreciate any answers, thanks!

First, I would try again to download via the bioconductor command. If your internet is intermittent it may just take a few tries.
If that doesn't work, I would then strongly recommend using BioCManager to install from a local source wherever you've downloaded the files. If you try to copy the files around there are a lot of ways that things can go wrong, and there are a lot of ways that it still might not work even if you put it in the right place, if you don't use the package manager as intended.

Great thanks! I have tried unzipping the .tar.gz file and copying it, and things indeed went wrong. Like mentioned, installing from local sources seemed to be the best way.


text:
Can the information regarding atomic bonds be interpreted only from atom lines in a PDB file?
The following are the first 14 atoms from the protein 1CRN:
ATOM      1  N   THR A   1      17.047  14.099   3.625  1.00 13.79           N
ATOM      2  CA  THR A   1      16.967  12.784   4.338  1.00 10.80           C
ATOM      3  C   THR A   1      15.685  12.755   5.133  1.00  9.19           C
ATOM      4  O   THR A   1      15.268  13.825   5.594  1.00  9.85           O
ATOM      5  CB  THR A   1      18.170  12.703   5.337  1.00 13.02           C
ATOM      6  OG1 THR A   1      19.334  12.829   4.463  1.00 15.06           O
ATOM      7  CG2 THR A   1      18.150  11.546   6.304  1.00 14.23           C

ATOM      8  N   THR A   2      15.115  11.555   5.265  1.00  7.81           N
ATOM      9  CA  THR A   2      13.856  11.469   6.066  1.00  8.31           C
ATOM     10  C   THR A   2      14.164  10.785   7.379  1.00  5.80           C
ATOM     11  O   THR A   2      14.993   9.862   7.443  1.00  6.94           O
ATOM     12  CB  THR A   2      12.732  10.711   5.261  1.00 10.32           C
ATOM     13  OG1 THR A   2      13.308   9.439   4.926  1.00 12.81           O
ATOM     14  CG2 THR A   2      12.484  11.442   3.895  1.00 11.90           C

These fourteen atoms are members of two THR residues.
Can the information regarding atomic bonds be interpreted only from these fourteen lines?
For example,

the first atom, N, has three bonds attached to CA.
the second atom, CA, has four bonds; three are attached to N, and one is to C.
the third atom, C, has four bonds; one is attached to CA, one to oxygen, and two are empty.
the fourth atom, O, has two bonds; one is attached to the third carbon, and one is to the fifth atom beta-carbon.
the fifth atom, beta-carbon, has four bonds; one is attached to the fourth atom, two are attached to the first branch (OG1), and one is attached to the second branch (CG2).

Or do I have to use my background knowledge of amino acids to find bonds?
Note: I am not talking about hydrogen bonds. I am talking about atomic bonds only in the atoms in the residue.

The official list of bonds in each residue in the PDB is in the Chemical Component Dictionary.
You could also use interatomic distances to infer connectivity, but that's less reliable.

You can use topology files to extract bond information.
Here is an example at CHARMM36


text:
Why does segmasker modify headers?
When parsing NCBI BLAST+ segmasker output, I realized the program modifies almost every header containing these strings: (  or  ) for no (known to me) reason.
In the majority of the strings mentioned above, the whitespace is removed:

input file:
>tr|Q9UXV4|Q9UXV4_PYRAB Pyruvate formate-lyase activating enzyme ( act-1) OS=Pyrococcus abyssi (strain GE5 / Orsay) OX=272844 GN=act-1 PE=4 SV=1,
output file:
>tr|Q9UXV4|Q9UXV4_PYRAB Pyruvate formate-lyase activating enzyme (act-1) OS=Pyrococcus abyssi (strain GE5 / Orsay) OX=272844 GN=act-1 PE=4 SV=1,
input file:
>tr|A0A5K1JWI6|A0A5K1JWI6_9APHY Cysteine protease ATG4 ) OS=Ganoderma boninense OX=34458 GN=Q523C3 PE=4 SV=1,
output file:
>tr|A0A5K1JWI6|A0A5K1JWI6_9APHY Cysteine protease ATG4) OS=Ganoderma boninense OX=34458 GN=Q523C3 PE=4 SV=1.

It isn't a significant change, but it complicates searching for headers in the output file.
However, the program is not consistent - if there is a comma (,) before the  ) string, the whitespace is not removed:

input file:
>tr|A0A3M9YA12|A0A3M9YA12_9PEZI Myeloid lymphoid or mixed-lineage leukemia 5 (Trithorax, ) OS=Verticillium nonalfalfae OX=1051616 GN=MLL5 PE=4 SV=1,
output file:
>tr|A0A3M9YA12|A0A3M9YA12_9PEZI Myeloid lymphoid or mixed-lineage leukemia 5 (Trithorax, ) OS=Verticillium nonalfalfae OX=1051616 GN=MLL5 PE=4 SV=1.

What's the reason for such a behaviour?
Is there any option to prevent it and keep headers unchanged?

edit.
The whitespace is also removed for:

  →  ,
 , → ,.

Moreover, it removes a semicolon when it is at the end of a header:
>ABCD/1-4 ABCD;ABCD; → >ABCD/1-4 ABCD;ABCD.

Weird... You could always run your input through something like sed 's/ *)/)/g; s/( */(/g' to remove the spaces yourself before using segmasker, but obviously that isn't ideal.


text:
Why am I not being able to detect all the bonds in a residue?
Let us consider the first seven atoms of the  protein 1CRN.pdb:
ATOM      1  N   THR A   1      17.047  14.099   3.625  1.00 13.79           N
ATOM      2  CA  THR A   1      16.967  12.784   4.338  1.00 10.80           C
ATOM      3  C   THR A   1      15.685  12.755   5.133  1.00  9.19           C
ATOM      4  O   THR A   1      15.268  13.825   5.594  1.00  9.85           O
ATOM      5  CB  THR A   1      18.170  12.703   5.337  1.00 13.02           C
ATOM      6  OG1 THR A   1      19.334  12.829   4.463  1.00 15.06           O
ATOM      7  CG2 THR A   1      18.150  11.546   6.304  1.00 14.23           C

This represents an amino acid called Threonine.
The general bond diagram for Threonine is as follows:

So, it has at least six bonds (I am actually not sure!) which are also evident from the seven ATOM lines.
Each atom has its radius. In general, if the distance between two atoms is less than the sum of the radii, there is a bond between these two atoms.
Based on the above information, I tried to identify the bond information, and I could detect only two:
     atom1           atom2        distance       atm1 Radi      atom2 Radi        sum radi        is bond?
         N              CA            1.50            0.65            0.70            1.35           False
         N               C            2.44            0.65            0.70            1.35           False
         N               O            2.67            0.65            0.60            1.25           False
         N              CB            2.48            0.65            0.70            1.35           False
         N             OG1            2.75            0.65            0.60            1.25           False
         N             CG2            3.86            0.65            0.60            1.25           False
        CA               N            1.50            0.70            0.65            1.35           False
        CA               C            1.51            0.70            0.70            1.40           False
        CA               O            2.36            0.70            0.60            1.30           False
        CA              CB            1.57            0.70            0.70            1.40           False
        CA             OG1            2.37            0.70            0.60            1.30           False
        CA             CG2            2.61            0.70            0.60            1.30           False
         C               N            2.44            0.70            0.65            1.35           False
         C              CA            1.51            0.70            0.70            1.40           False
         C               O            1.24            0.70            0.60            1.30            True
         C              CB            2.49            0.70            0.70            1.40           False
         C             OG1            3.71            0.70            0.60            1.30           False
         C             CG2            2.98            0.70            0.60            1.30           False
         O               N            2.67            0.60            0.65            1.25           False
         O              CA            2.36            0.60            0.70            1.30           False
         O               C            1.24            0.60            0.70            1.30            True
         O              CB            3.12            0.60            0.70            1.30           False
         O             OG1            4.34            0.60            0.60            1.20           False
         O             CG2            3.74            0.60            0.60            1.20           False
        CB               N            2.48            0.70            0.65            1.35           False
        CB              CA            1.57            0.70            0.70            1.40           False
        CB               C            2.49            0.70            0.70            1.40           False
        CB               O            3.12            0.70            0.60            1.30           False
        CB             OG1            1.46            0.70            0.60            1.30           False
        CB             CG2            1.51            0.70            0.60            1.30           False
       OG1               N            2.75            0.60            0.65            1.25           False
       OG1              CA            2.37            0.60            0.70            1.30           False
       OG1               C            3.71            0.60            0.70            1.30           False
       OG1               O            4.34            0.60            0.60            1.20           False
       OG1              CB            1.46            0.60            0.70            1.30           False
       OG1             CG2            2.54            0.60            0.60            1.20           False
       CG2               N            3.86            0.60            0.65            1.25           False
       CG2              CA            2.61            0.60            0.70            1.30           False
       CG2               C            2.98            0.60            0.70            1.30           False
       CG2               O            3.74            0.60            0.60            1.20           False
       CG2              CB            1.51            0.60            0.70            1.30           False
       CG2             OG1            2.54            0.60            0.60            1.20           False
2
Done!

How can I detect all bonds?
EDIT: As per pippo1980's request:
public class ADTResidue: Dictionary<int, ADTAtom>, IEquatable<ADTResidue>
{
    public string ProteinName { get; set; }
    public string ChainName { get; set; }
    public int ResidueNo { get; set; }
    public Point3d Coordinate { get; private set; }
    public string Code { get; set; }
    public string Full { get; set; }
    public Color ColorRGB { get; set; }

    public ADTResidue(string code, string full, Color colorRgb)
    {
        Coordinate = new Point3d();

        Code = code;
        Full = full;
        ColorRGB = colorRgb;
    }

    private readonly Dictionary<int, CovalentBond> BondDictionary;

    public List<CovalentBond> GetBonds()
    {
        return new List<CovalentBond>(BondDictionary.Values);
    }

    public ADTResidue()
    {
        BondDictionary = new Dictionary<int, CovalentBond>();
    }

    public void CalculateBonds()
    {
        foreach (ADTAtom atom1 in this.Values)
        {
            foreach (ADTAtom atom2 in this.Values)
            {
                if (atom2.Radius <= 0)
                {
                    throw new Exception("zero radius in " + atom2.Name);
                }

                if (atom1 != atom2)
                {
                    double distance = Point3d.Distance(atom1.Coordinate, atom2.Coordinate);

                    bool exists = ABond.IsBondExists(atom1.Coordinate, atom2.Coordinate, atom1.Radius, atom2.Radius);

                    if (exists)
                    {
                        int hashCode = atom1.GetHashCode(atom2);
                        if (!BondDictionary.ContainsKey(hashCode))
                        {
                            BondDictionary.Add(hashCode, new CovalentBond(atom1, atom2));
                        }
                    }
                }
            }
        }

        int count = BondDictionary.Count;

        Console.WriteLine(count);
    }

    #region equality comparison
    // IEquatable<> implementation...
    #endregion
}

public abstract class ABond: Pair<ADTAtom, ADTAtom>
{
    public ABond():base()
    {
    }

    public ABond(ADTAtom first, ADTAtom second) : base(first, second)
    {
    }

    public static bool IsBondExists(Point3d first, Point3d second, double radius1, double radius2)
    {
        double distance = Point3d.Distance(first, second);

        return (distance <= (radius1 + radius2));
    }
}

public class CovalentBond : ABond
{
    public CovalentBond():base()
    {
    }

    public CovalentBond(ADTAtom first, ADTAtom second) : base(first, second)
    {
    }
}

would you share your code with us ?

@pippo1980, its a part of a bigger framework. I am sharing a smaller part of it. Hope that helps.

unfortunately not, but could help others. Did you try running it on first T of 3NIR.pdb ?

@pippo1980, Yes. My code works now. Previously, I made mistakes with radii values.

is that above C#    ?

@pippo1980, Yes.

is it same as java ?

@pippo1980, Not same, but similar.

Your method doesn't work for two possible reasons;

The radius information may be incorrect.
Your molecule requires some minimization. This option helps to optimize distance between atoms.

You can find some information about the minimization at the following links

Wikipedia
Gromacs


text:
Help sought with GWAS and vcf files, lack phenotype labels
This question has also been asked on Biostars
Hi, I am very new to this area, and I am taking a class about bioinformatics. For an independent project assignment, I need to do a GWAS. I am using the bash terminal. I downloaded all the fastq I need,

trimmed the files,
converted them into bam/sam
then performed vcf
and then bed/bim/fam etc.

However, when I tried to perform GWAS in plink, I realized I dont have phenotype data. It supposed to have two phenotypes.
There are two groups/phenotypes of fastq files, each containing 29 samples. Let's say they are group 1 and 2.
For each group, I converted every fastq to sam then bam, then I combined 29 bam to one bam. Then I combined two bams (for the two groups) together to a vcf.gz.
Problem. There is no phenotype data in the following plink files.
Questions
What step has been wrong?
Or what I should do to incorporate the phenotype data?
I dont have to be perfect (like the QC steps), and I cannot understand too complicated codes. I am seeking a simple output such as a Manhattan plot or alternative pipeline.

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.

Please show the command lines you have run. You've asked a question about what step was wrong, but that's difficult to work out without knowing precise details about what you did.

Do you have phenotypic data associated with the samples? As noted on the biostars answer, you might need to not combine BAMs such that the relation of each sample with phenotype is maintained.


text:
Can I perform Expression Analysis on this 10x dataset?
10x Genomics: 20k Human PBMCs is the dataset.

Description of the dataset:

Inputs/Libraries
Human peripheral blood mononuclear cells (PBMCs) of a healthy male donor aged 30-35 were obtained by 10x Genomics from AllCells.

Gene expression and V(D)J libraries were generated from ~33,000 cells (18,470 cells recovered) as described in the Chromium Next GEM Single Cell 5' HT Reagent Kits v2 User Guide (CG000421) using the Chromium and sequenced on an Illumina NovaSeq 6000 to a read depth of approximately 25,000 mean reads per cell for Gene Expression, 15,000 mean reads per cell for TCR Amplified, and 25,000 mean reads per cell for BCR Amplified libraries.

Paired-end, dual indexing

Read 1: 26 cycles (16 bp barcode, 10 bp UMI)
i5 index: 10 cycles (sample index)
i7 index: 10 cycles (sample index)
Read 2: 90 cycles (transcript)

This is the RPubs Link for the Single Cell Analysis work I've performed on it.
Can someone confirm that I can do both Single-Cell RNA Seq Analysis AND Differential Gene Expression Analysis?
I know that Expression Analysis comes from having a control and a treatment group. I honestly don't know anything about bioinformatics but I'm diving into the fire to get a solid understanding once I've coded everything.
I followed this video and want to follow this other video.
Thank you guys.

@M__ I've updated the post with the name of the file I performed analysis on. Please let me know if this is helpful.

Thanks @Antonio answered below

Typically what this PMBC RNA-seq data set would be used for is for a control against individuals who would exhibit a pronounced cellular immune response. For example, following a COVID19 infection the PMBCs would be have lots of T-cell CD4+ and CD8+ and are essential for controlling the infection. Thus it is a really good approach to understand the reasons for different clinical outcomes, such chronic infection and severe infections.
I would suspect a cancer patient would show a pronounced shift in cellular immune response, but I'm not certain. Anyway, by comparing a control data set against a COVID19 patient a better understanding of the cellular immune response can be obtained and a differential gene expression analysis would be an important inroad to assess the T-cell up-regulation. DEG would show very pronounced differences between a healthy individual and a patient.
Having a hypothesis about the expected results from DEG is important, because interpreting DEG output is otherwise difficult.
With regards your question, a DEG analysis IMO would be a good analysis and is easy to perform, however there only appears to be a single sample. A DEG analysis in context is control X vs patient Y. The sample you have appears to be a control sample so you would need a second sample. Performing DEG is easy, but do stress that interpreting DEG isn't easy because everything is relative. Thus what is sought usually in PMBC comparisons healthy vs patient is evidence of T-cell CD8+/CD4+ up regulation, any indication about which possibly which HLA type is active (I dunno whether RNAseq will do this - good if it can) and information on the signalling cascade. Basically as much information on the cellular immune response that can be gleaned via RNAseq without performing laborious and difficult immunological experiments requiring specialist apparatus.
There are investigators here with more experience than I in eukaryotic RNAseq analysis, so there input would be helpful.


text:
segmasker doubles found intervals
During NCBI BLAST+ segmasker output exploration, I found out that for some sequences, the program returns the same interval twice:
>tr|A0A0D2B4V6|A0A0D2B4V6_9EURO Uncharacterized protein OS=Cladophialophora immunda OX=569365 GN=PV07_04097 PE=3 SV=1
131 - 139
131 - 139

This seems not to be a planned behaviour, but maybe someone has an explanation.
The command I ran was:
segmasker -in <input_file> -out <output_file> -locut 1.5 -hicut 1.8 -window 15

Exemplary sequences triggering such an event:
>sequence1
MDSHNSSTLIILSLPPKTFVGLDLVSFDSSPNFQGIAKIPVGIHFLYTGTDASLSIRHGRWLNLTSTHPIQVLRWNSDGETLQPVDQNDQTARNAIAAISNGRGLVDYGALQSATSELSTRESQSHDASSLAAAAAAAAADEDRAESTDWPALTAHISSPLLTRILSPDWIVSSISSAPTDTETIPGLSHLEASNALHQLPLDLLPVNLKQTWAGGDIGRARTERARDRSWYLGQLIEQAARPQGGDKTAGAREVLGELQFCFLMVLTLANYSCLEQWKRLLSVFFTCRTALDEVEGYFVQVIRVLRLQIKHVDDVEGGLFELRDESGSAWLRTLWARFRSVVDDADKGERVSLKKEVDALEKLFEGKYGWQSERDILRRGMLELEDGERVEVTMHGVDEDEETGEYAPVIVET
>sequence2
VDKLFKKNAMGNNIFHEIAIQGSLMLLQTIRDNVDGQMDTYLSEINDQNESCITIVADRHRGLLAKDLIEIFVDLGADINETDNWGNTALHYSVFHRDHVLAEWLCQQPGINLNAINNDEHTPLGLAILWNIQSIKTLLESKGARFNNIEFSDSETSDDDDDDDDDGHDENDE

What command are you running? What options?

Right, I should have include it. segmasker -in <input_file> -out <output_file> -locut 1.5 -hicut 1.8 -window 15.

Could there be something silly like having duplicate sequences in your input? Which are then combined in the output?

Nah, I've just run segmasker on the given example and got.

I'm using segmasker locally. Installed last week.


text:
Keep first 20 nucleotides of sequences using R
I have an Excel sheet from the CRISPR library where I have sequences of gRNAs (with 30 nucleotides) in a column and I only need to keep the first 20 nucleotides for those gRNAs and delete the rest nucleotides. The column has 1000 rows that have 30 nucleotide gRNA sequences. How can I delete extra nucleotides using R

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 will need to "read" your data. There are packages that can read Excel files directly but I would save the data as csv first and then use read.csv(). The following would create a new column with your desired output:
data_frame$new_col <- substring(data_frame$col_with_seqs, 1, 20)

Alternatively you can just use Excel's LEFT() since your data is already in Excel:
=LEFT(A1, 20)


text:
Matches, mismatches and indels
I have a bam-file of reads and reference genome. How can I located matches, mismatches and indels in three different files with Phred score in each data of the files?

Try one of these:
GATK:
https://gatk.broadinstitute.org/hc/en-us/articles/360036194592-Getting-started-with-GATK4
freebayes:
https://bioinformaticsworkbook.org/dataAnalysis/VariantCalling/freebayes-dnaseq-workflow.html#gsc.tab=0


text:
How to extend the x axis in Dimplot Seurat
I have performed a Seurat PCA via Dimplot.
How do I extend the x axis?
As you can see in my figure the double x axes overlap.

Could you please add some context around your question, so we can understand why you want to do this? What are your trying to produce? How did you generate this plot?

Looking at the source code for Seurat's DimPlot() it's based on ggplot2 graphics (https://github.com/satijalab/seurat/blob/HEAD/R/visualization.R#L794), so you can manipulate the dimplot object and save it using ggplot2 commands. Using ggsave() you can specify the figure dimensions; if you increase the width of the final figure it should 'fix' your label overlap issue, i.e.
#BiocManager::install("Seurat")
library(Seurat)
#> Attaching SeuratObject
library(ggplot2) # for plotting
library(cowplot) # for the theme used by Seurat

# example data
data("pbmc_small")
p <- DimPlot(object = pbmc_small, split.by = "ident")

# example dimplot
p

# increase font size to make the labels overlap
p +
  theme_cowplot(font_size = 28) +
  theme(strip.background = element_blank())

# save the 'overlapping label' figure and specify a relatively large width to 'space out' the three facets
ggsave("dimplot_test.png", width = 40, height = 10, units = "cm")

You could also decrease the font size and keep the original dimensions to fix the problem, or alter the breaks on the x axis, e.g.
# reduce the size of the x axis text labels
p +
  theme_cowplot() +
  theme(strip.background = element_blank(),
        axis.text.x = element_text(size = 8))

# change the number of breaks
p +
  theme_cowplot() +
  theme(strip.background = element_blank()) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 2))

# Specify the breaks you want to use
p +
  theme_cowplot() +
  theme(strip.background = element_blank()) +
  scale_x_continuous(breaks = c(-20, 0, 20, 40))

Created on 2022-12-01 with reprex v2.0.2

As a direct answer to your question, restrict the output to producing a single column:
Dimplot(..., ncol=1)

But this doesn't seem like a useful answer. Please provide more context to your question if you want more help than this.


text:
Could someone help me understand if all sample within this are control group?
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE117872
I'm teaching myself sequence analysis.
This is the count data:

and then this is, what I believe to be, the coldata:

All I want to know is if I can perform Differential Gene Expression analysis on this using just the drug_status column since it seems like the second dataset already has results to it.
What would be the reference if listings are either Sensitive, Resistant, or Holiday?
Thank you.


text:
What are the conditions to perform Gene Expression Analysis on data you've performed RNA Sequence Analysis on?
I want to perform Differential Gene Expression Analysis. I recently completed RNA Sequence Analysis using the Seraut Tutorial up to this point: Finding Differentially Expressed Features (Cluster Biomarkers).
My project can be found at this Rpubs Link which contains the link to the dataset on the 10x genomics website.
My question is... is the following table derived from performing the steps found in the first link I provided?

If so, once I derive that table ( ID, Name, L2FC, p-value, etc.), does that mean I can perform Differential Gene Expression Analysis?
I am a complete newbie and am teaching myself all of this. Thanks, y'all.
UPDATE:
Yes, it did provide me that table. With this table, am I able to perform gene expression?

The table pictured here is the result of a differential gene expression (DGE) test. In the table, each row is a gene (labeled with the gene ID and name) and each column (besides the first two) are the results of DGE tests. Each cluster has two sub-columns: L2FC (log2 fold-change) which tells you the average relative expression changes between cells of the cluster vs. others; and p-value, which is how significant the change is.
You can make this table yourself by performing DGE tests for all clusters. You've already done this yourself in your Rpubs script! To make it easier though, so you don't have to manually test each cluster, use the FindAllMarkers function. It will compare each cluster to all other clusters and report the L2FCs and p-values.
However, I recommend using the p_val_adj (adjusted p-value) instead of the p-value. This is due to a problem called multiple testing: you are performing many, many tests (each gene in each cluster is a separate statistical test), and we expect 5% of differentially expressed genes (DEGs) to be DE by random chance. Over something like 100,000 tests, that means you'd expect 5,000 of those to be false positives (which is huge). So p-value is adjusted to account for multiple tests.
Once you filter the DEG table (p-adj<0.05), you can do other downstream analyses to describe what has changed between the conditions. You already have a good start in your Rpub script with the violin plots. Cluster 6 is clearly distinct in terms of those genes' expression levels. You could try to figure out what's special with that cluster. Good luck!

With my violin plots, does that mean that those two genes , c("MS4A1", "CD79A"), are found within multiple clusters?

@Antonio Each violin plot shows the expression of that gene in all cells (the black dots). As you see, most values are 0 (which is very common for scRNA-seq), along with some random low-level noise. The important bits to look at is anything above the noise level (in this case, ~0.5). Any violin with color is usually somewhat significant. For the MS4A1 gene, you can see the only clear violin is in cluster 6, while CD79A is in 6 and 20. If you look at your UMAP, cluster 6 is off to the side by itself. This may be a gene that is cell-type specific or only expressed in certain conditions.


text:
How to solve Nextflow error: "Trace file already exists"?
When trying to run epi2me-labs/wf-artic, I get the following error:
❯ nextflow run epi2me-labs/wf-artic \
 -r v0.3.18 --fastq  ~/Downloads/barcode95.fastq.gz \
 --scheme_version Midnight-ONT/V3

N E X T F L O W  ~  version 22.10.1
Launching `https://github.com/epi2me-labs/wf-artic` [ecstatic_marconi] DSL2 - revision: a60a1e1e73 [v0.3.18]
Trace file already exists: /Users/corneliusromer/Downloads/BJGMC_ONT_Run71/output/execution/trace.txt

It seems that nextflow doesn't like to rerun after it failed. Is there a convenient way to force nextflow to rerun, in a similar way to snakemake --force-rerun? Or do I need to manually delete the trace file every time? That works, but it's tedious.

Rotation of report files (trace, execution, timeline etc) was removed recently (since version 22.10.0 in a762ed5). To overwrite the existing trace file you will need to add the following to your Nextflow config:
trace.overwrite = true

See: https://github.com/nextflow-io/nextflow/issues/3317


text:
Command line based interaction with the UCSC genome browser, automated data upload and track generation
Is there a way to submit data to ucsc genome browser for result/track display, on the fly, as a local bedgraph file with example contents like below:
browser position chr19:49302001-49304701
browser hide all
browser pack refGene encodeRegions
browser full altGraph
#   300 base wide bar graph, autoScale is on by default == graphing
#   limits will dynamically change to always show full range of data
#   in viewing window, priority = 20 positions this as the second graph
#   Note, zero-relative, half-open coordinate system in use for bedGraph format
track type=bedGraph name="BedGraph Format" description="BedGraph format" visibility=full color=200,100,0 altColor=0,100,200 priority=20
chr19 49302000 49302300 -1.0
chr19 49302300 49302600 -0.75
chr19 49302600 49302900 -0.50
chr19 49302900 49303200 -0.25
chr19 49303200 49303500 0.0
chr19 49303500 49303800 0.25
chr19 49303800 49304100 0.50
chr19 49304100 49304400 0.75
chr19 49304400 49304700 1.00

Normally, you can input the data in the data box like below:

and submit - and display results after another [GO] button that's pushed manually on the following web page to produce the output below:

I am hoping there is a command line option for this where it takes the data input above as a local file name and generates the output that normally one can generate with a few manual clicks (including pasting of the data text into the window or upload of the data file) on the ucsc web page.
Thanks in advance


text:
Coexpression network analysis
I have a question related with the most popular frameworks for network construction: WGCNA for weighted networks and igraph package for unweighted networks.
What are the advantages or disadvantages of applying these approaches and when it is recommended to use one over another.

Hi, this is a very open-ended discussion question that doesn't work well for the Stack Exchange format. Do you have a particular application in mind that could be linked to your question to make it more targeted?

I suspect the OP has been recommended an igraph recipe or concept for an unweighted graph theory calculation of co-expression  and they believe this represents a "plug and play" package like WGCNA (albeit thats a guess). The answer is they are not really the same thing.

igraph is a serious general network theory framework for any given data science application. It permits full directionality and describes vertices (nodes), edges (branches)  and does cool stuff like De Bruijn graphs and Laplacian graphs. Understanding graph theory would be helpful and will be a steep learning curve. It's a serious data science package.
WGCNA is a specific ready made RNA expression package which appears to be network theory and is undirectional, i.e. no directionality which is fine for RNA expression. It will have far less overall power but if you've a load of RNA-seq data and want 'plug and play' type stuff WGCNA is the only way. It's a correlation network, based on r2 (Pearson's correlation), to assess co-expression. I would assume the weighting is geared towards identifying co-expression. There is an error in Fig. 2 legend BTW Langfelder & Horvath (2008): the legend states r^2 what they mean is r2. I guess they messed up on their LaTeX.
If this is a one of coexpression analysis WGCNA. If this is a big project, graph theory is central to that effort and you're looking towards developing a data science tool set then its igraph. There are other graph theory packages BTW.

In case my answer wasn't clear igraph will do weighted edges via is_weight. The information is here, basically it will do anything.
If your question is "what does weighting contribute to coexpression analysis", that is a good question. I know graph theory, I dunno the biology of gene expression versus the maths. What WGCNA is doing is statistically is simple if weighting is excluded ... my reasoning is it must be pivotal to its success.

To address the OPs comments.
Statistical power In the pure sense, I was incorrect to use this term because technically 'lack of power' is the inability to reject a hypothesis when it is false. For example, parametric statistics are more powerful than non-parametric statistics for this reason.
What I meant was two things:

In graph theory you can do anything within the "network", what is called a "graph", e.g. find the centre of the graph, minimal distance across the "graph", cyclisation, anything you can think of.

1a. WGCNA is fixed it only does what is described in Langfelder & Horvath (2008). Nothing can be done outside that without rewriting the underlying code and it would be better to use igraph instead if that was the objective in my opinion.
Thus igraph is more "powerful" because it will do a lot more stuff.

My concern of WGCNA is its over-reliance on r2. This is a parametric statistic and may not always be applicable but from my reading of the paper (albeit briefly) it was entirely dependent on that. To centre everything on r2, personally  I would at a minimum use r in addition (Spearman's correlation). This is extremely complex data and there is a risk of over-simplification. r  and r2 can certainly give different results.

The answer to the question however it is not what I think but everyone else thinks and Langfelder & Horvath (2008) has almost 10000 citations. I am not aware that the calculation has been overturned, e.g. DSeq2 has been overturned and demonstrated to be inappropriate. So I would assume it must be capturing the overall picture, so the positives outweigh potential negatives.

thanks @M__ for your insights into both packages, it has been very useful and  it has been explained in a very concise way the objectives,  do you know why WGCNA present less statistical power.

@MaríaJosé I've updated updated the post.


text:
Struggling to install scvi-tools on Mac M1
I was wondering if somebody here has experience with installing scvi-tools on a Mac with M1 chip. I tried the following command:
conda install scvi-tools -c conda-forge

But got the following error:
Collecting package metadata (current_repodata.json): failed
CondaSSLError: OpenSSL appears to be unavailable on this machine. OpenSSL is required to
download and install packages.
Exception: HTTPSConnectionPool(host='conda.anaconda.org', port=443): Max retries exceeded with url: /conda-forge/osx-64/current_repodata.json (Caused by SSLError("Can't connect to HTTPS URL because the SSL module is not available."))

Then I tried to install OpenSSL with following command:
conda install -c anaconda openssl

and I got the exact same error as the conda install scvi-tools -c conda-forge command.
I do not understand why collection of metadata is failing for openssl. Any inputs would be highly helpful.
Thanks!

There might be system incompatibilities, thats what my error log suggests. Thus I was unable to easily install as well. However this might be a Python version  issue (could be Python 3.6, I used 3.11 [not the best idea]). However the error log should have flagged that. scvi-tools is a bioinformatics tool, but please try and edit the question with the biological problem. I'd guess this is RNA-seq PCA type analysis

Wait, how did you end up with a non-working conda install on that M1 mac in the first place?  Looks like it's not complaining about a missing package dependency; it's complaining about not having openssl to do HTTPS downloads within conda.

Did you try:
pip install scvi-tools

?

Would need venv (virtual env) with pip, but I agree I think this is the best approach in this particular scenario

https://scvi-tools.org/


text:
How to make bar bigger in UpsetR plot
I have tried to use UpsetR to visualize the input file which can be found here.

library("UpSetR")
orthogroups_df<- read.table("orthogroups.GeneCount.tsv",  header=T, stringsAsFactors = F)

#All species
selected_species <- colnames(orthogroups_df)[2:(ncol(orthogroups_df) -1)]
selected_species
 [1] "Atha" "Cann" "NQLD" "Natt" "Ngla" "Nlab" "Nsyl" "Ntab" "Ntom" "Slyc" "Stub" "Vvin"

head(orthogroups_df)
  Orthogroup Atha Cann NQLD Natt Ngla Nlab Nsyl Ntab Ntom Slyc Stub Vvin Total
1  OG0000000    0    0  965    0    0    3    0    0    0    0    0    0   968
2  OG0000001    0    1    3    0    0  448    0    0    0    0    0    0   452
3  OG0000002    0    1  313    0    0  120    1    0    1    0    0    0   436
4  OG0000003    0   93   15   21   46   16   33   63   36   25   39   26   413
5  OG0000004    1   42    2   34  109    6    8  154   11    9    4    0   380
6  OG0000005    0    2   61    1   34   44   91   70   43   20    1    0   367

ncol(orthogroups_df)
[1] 14

orthogroups_df[orthogroups_df > 0] <- 1
# we only show intersections of interest
intersections=list(list(selected_species),
                   list("Nlab", "NQLD", "Ngla"),
                   list("Nlab", "NQLD", "Nsyl"),
                   list("Nlab", "Ngla", "Nsyl"),
                   list("NQLD", "Nsyl", "Ngla"))
upset(orthogroups_df,
      text.scale = c(1.4, 1.4, 1.4, 1.4, 1.4, 1.4),
      sets=rev(selected_species),
      nsets = ncol(orthogroups_df),
      #keep.order=T,
      #mb.ratio=c(0.5,0.5),
      #order.by='degree',
      #order.by='freq',
      intersections = intersections,
      sets.x.label="Total nbr of orthogroups",
      mainbar.y.label = "Nbr of orthogroups",
      #empty.intersections = "on",
      group.by = "intersections")

How is possible to make the bar bigger of the double digits values or is a way to rotate the plot?

I'm not quite sure what you mean.  You can alter the % of height that the bottom part of the graph takes up, that will make all the vertical bar parts proportionately taller, but I don't think you should because you have so many groups to display in the bottom part.
I think what you'd want to do is to alter the y-axis so that it stops being totally proportional, to compress the size of the tallest bar, so you can better see the difference in bar sizes between the other comparisons?
I'm not sure you can do that.  Making the y axis logged with "scale.intersections" might help.
You'd probably be best off making a second graph that omits that "common to all" intersection.  Once that is out of the picture, the difference between the others will be more visible.


text:
Maximum likelihood estimation for tree construction
Could someone direct me to a resource (textbook chapter, lecture notes, online video, etc.) that demonstrates, in a mathematically or computationally rigorous way, how to use maximum likelihood estimation to generate a phylogenetic tree, given a set of aligned sequences?
I've found several resources which describe the procedure for maximum likelihood estimation, such as:

giving the formula and the parameters that need to be determined (which I already know);
how to use it for aligned sequences (I already know this),
I am seeking a worked-out example.

Alternatively, if someone is willing to work out a small example, how would you use maximum likelihood estimation to generate a tree for the following three sequences?
organism_1 ATCG
organism_2 ATCC
organism_3 ATTT

@DataScienceNovice ... here's a rough explanation. Just to note a tree with three taxa doesn't really work as a tree - you need 4 taxa.

If we have tree 1 ((A,B),(C,D)) the probability of the tree is expressed as the summed log likelihood for each position 1, 2, 3 .... in the alignment. The probability is derived from a mutation matrix, thus certain mutations e.g. transitions A<->G; C<->T are more probable (high likelihood), e.g. transversions A<->C; A<->T; G<->C; G<->T are less probable in accordance with the matrix. This hypothetical matrix (most matrices) is scoring transitions as much more probable than transversions.

The 'total probability' expressed  log likelihood for each column, i.e. each alignment position is added together leveraging an interesting property of logs (see final sentence of this paragraph). That becomes the 'score' or 'likelihood' for that tree. Note in the example you provided you have 4 alignment positions, so 4 probabilities of the log-likelihood will be added. Note, the probabilities are actually multiplied but log-likelihoods added is the same thing.

We will now rearrange the tree lets try tree 2 ((C,D),(B,G)), again the hypothetical matrix scores transition/transversions, thus essentially if this tree has higher transitions than the tree in point 1. that will provide a higher likelihood. The tree that maximises the transitions as sister groups will maximise the likelihood. Equally maximising transversions as the most divergent groups will maximise the likelihood of that tree.

Re-arrange the tree again tree 3 ... ((C,B),(D,A)) ... which tree has the highest transitions vs transversions? ...

In this example Tree 2 has the highest likelihood. Thus that is the maximum likelihood. Simple?

In the example you have the outcome will depend on the matrix. Alignment position 3 comprises a transition, whilst position 4 comprises a transition and a transversion. Generally transversions are 1:2 to 1:4 less likely even though there are twice as many. So the example you gave the 'tree' (you can't really have a tree with 3 taxa) would be ((organism2,organism3),organism1). This is against the total number of mutations because organism2 and organism3 has two mutations whilst organism2 and organism1 only have one, but the mutation types are different (also see point 3 below). In other words the occurrence of a transversion between organism1 against organism2,organism3, mean organism2 and organism3 are more likely sister groups even though they comprise 2 mutations (transitions) against the single mutation between organism1 and organism2.
A good maximum likelihood program for phylogenetics is RAxML.

The reality is it's not that simple because the mutation matrix is a function of the tree, which is a function of the mutation matrix and in technical tree building the matrix is also calculated alongside the total probability of the tree via maximum likelihood. Just to be clear it can be calculated as a separate maximum likelihood calculation and then fixed for the tree calculation, but it's not a good idea because of non-linear interactions.

In addition, mutation rate (heterogeneity) is also calculated via maximum likelihood.

The final thing to note is back-mutations are also calculated and that refers to past posts on this site, which I will link if I find them.

Complexity If the parameters are calculated via maximum likelihood simultaneous to the tree the parameter space becomes EXTREMELY large and the calculation is searching the tree landscape against the parameter landscape which becomes a multi-dimensional calculation.
Reality In reality much of a maximum likelihood calculation becomes about the search strategy because its impossible to traverse all of the parameter space so a heuristic estimate needs to be performed and the calculation can get stuck on lower plateaus of probability and fail to identify the tree maximum likelihood. This is the nuts and bolts of real world tree building.
Note However if you fix the matrix and the rate, which a lot of people do, then points 1 to 5 describe the maximum likelihood tree.


text:
Downloading genes from NCBI in fasta format
I'm pretty new in bioinformatics.
I need to download FASTA sequences of several genes. The list of genes I have assembled consists of 140 genes, so I'd rather do this through via code than download each gene manually from the NCBI database.
The genes belong to 2 different organisms (bacteria) I have chosen, most of them are found in both organisms and some are specific to one or the other bacterium.
I've tried running this code which I found here: https://www.ncbi.nlm.nih.gov/books/NBK179288/
esearch -db nucleotide -query "aminoacyl-tRNA\ hydrolase" | efetch -format fasta

The gene of interest being aminoacyl-tRNA hydrolase here. But I don't know how to filter it for the specific organism, Neisseria meningitidis. Tried running this, without luck:
esearch -db nucleotide -query "aminoacyl-tRNA\ hydrolase" | efilter -organism [Neisseria meningitidis] | efetch -format fasta

Appreciate any help :)

Hi how does the list of 140 gene you are interested in look like? do you have gene accessions, or Ids or just gene symbols?, ie symbol: BRCA1 ; id :
672 ; accsession NC_000017.11 , are all the same gene.

what output do you get from your commands? nothing, an error, the wrong gene?

Yes, please [edit] your question and add a few lines of your input file. Presumably, you have a list of genes, so please show us the list so we can understand the input format.

@user12256545 Right now I have the list on a word-document with locus name and gene-ID (that would be "pth" and 61281026 for aminoacyl-tRNA hydrolase). All genes seem to have the same accession nr. so that won't be useful unfortunately.

If I put the species name in quotes, I get expected (by me) output:
% esearch -db nucleotide -query "aminoacyl-tRNA\ hydrolase" | efilter -organism "Neisseria meningitidis" | efetch -format fasta
# outputs:
>NZ_JWNU01000010.1 Neisseria meningitidis strain Nm440716 Scaffold10, whole genome shotgun sequence
GACCGCCTATTCGCTTGCCGCCGGCGGCCCCATCATGCAGGCAGGATTACACGCCTTCACACTCGTACCC
# [snip, tons of output]

Note that filtering down to only N. meningitidis still leaves all of the different sequenced accessions of this species, such that you will still have a lot of different annotated copies of this gene. You will need to filter more specifically to reduce that.

This was helpful, thank you! What would you suggest I use as identifier to filter more specifically? I tried substituting "aminoacyl-tRNA\ hydrolase" with gene-ID (61281026) but that doesn't seem to work.

@Vjon I would suggest trying an NCBI taxid in place of the species name.


text:
Breseq error (code 137). Any ideas?
The code I ran was here, so nothing fancy.
breseq -l 110 -o RifR_align -r Big_burk_assembly.fasta RifRNano_nanopore.fastq.gz

This was the output.
NOW PROCESSING Read alignment to reference genome
[system] bowtie2-build -q "RifR_align/data/reference.fasta" "RifR_align/02_reference_alignment/reference"
    [system] bowtie2 -t --no-unal -p 1 -L 31 --ma 1 --mp 3 --np 0 --rdg 2,3 --rfg 2,3 --ignore-   quals --local -i S,1,0.25 --score-min L,1,0.9 -k 2000 --reorder -x "RifR_align/02_reference_alignment/reference" -U "RifR_align/01_sequence_conversion/RifRNano_nanopore.converted.fastq" -S "RifR_align/02_reference_alignment/RifRNano_nanopore.stage1.sam" --un "RifR_align/02_reference_alignment/RifRNano_nanopore.stage1.unmatched.fastq"
    Time loading reference: 00:00:00
    Time loading forward index: 00:00:00
    Time loading mirror index: 00:00:00
    Killed
    Use of uninitialized value $l in scalar chomp at /home/lts/miniconda3/envs/assemblers/bin/bowtie2 line 588, <BT> line 4626.
Use of uninitialized value $l in substitution (s///) at /home/lts/miniconda3/envs/assemblers/bin/bowtie2 line 589, <BT> line 4626.
Use of uninitialized value $l in print at /home/lts/miniconda3/envs/assemblers/bin/bowtie2 line 593, <BT> line 4626.
(ERR): bowtie2-align exited with value 137
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!> FATAL ERROR <!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Error running command:
[system] bowtie2 -t --no-unal -p 1 -L 31 --ma 1 --mp 3 --np 0 --rdg 2,3 --rfg 2,3 --ignore-quals --local -i S,1,0.25 --score-min L,1,0.9 -k 2000 --reorder -x "RifR_align/02_reference_alignment/reference" -U "RifR_align/01_sequence_conversion/RifRNano_nanopore.converted.fastq" -S "RifR_align/02_reference_alignment/RifRNano_nanopore.stage1.sam" --un "RifR_align/02_reference_alignment/RifRNano_nanopore.stage1.unmatched.fastq"
Result code: 35072
FILE: libbreseq/common.h   LINE: 1411
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!> STACK TRACE <!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Backtrace with 6 stack frames.
breseq(+0x44cf4) [0x56102c3c3cf4]
breseq(+0x4a4f1) [0x56102c3c94f1]
breseq(+0x6d142) [0x56102c3ec142]
breseq(+0x36cc4) [0x56102c3b5cc4]
/lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf3) [0x7f71ba61e0b3]
breseq(+0x40559) [0x56102c3bf559]
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

The programs seemed to have no trouble aligning the reference genome, however the Bowtie2 package is throwing error code 137 and terminating the program.

That first Killed message you get suggests the process was killed by the system because it was taking too much RAM. It could be that all the other errors are just because the first was killed. Does it run through if you use a smaller dataset?

@LiamTSullivan please do consider the comment and the answer of MaximillanPress via upvote/accept.

[Not quite an answer but too long for a comment]
I am guessing that your assembly Big_burk_assembly.fasta is a rather large assembly. Assembly indexing can be very memory intensive for complex assemblies. If it is eating enough memory to kill the bowtie2 indexing process (as noted in comment by @terdon), I'd consider filtering it, for example if there are lots of very short contigs that are unlikely to be meaningful, consider removing those.
Just looking at the way you've named the assembly file, it looks like it might be a somewhat ad hoc assembly, so maybe there is a meaningful way to filter it. We might be able to suggest more disciplined ways to do this if you tell us more about it.
Alternately, just get more memory as suggested in comments by @terdon. This will be hard if you are running on a personal computer, but on a cluster you can usually specify this in a job parameter file. Here (section "submit scripts") is how to set it on SLURM, for example.


text:
How to make unrooted tree for Likelihood mapping result by using IQ2-tree?
I am a biologist, and I do not fully understand the tree topology of the experimental species.
I used four-taxon set (4 sequences) to identify the Four-cluster Likelihood-Mapping by using
iqtree2 -s ***.fa -lmap 2000 -n 0

I have obtained the svg file - triangle (attached - 1), however there were no data mentioned about T1, T2 and T3 unrooted tree based on my input 4  sequences  (Haur_2, Shad_6, Equa_oga2 and apa). I used single orthologous copy genes as the supermatrix for this data.
Therefore, I do not understand which topology of tree is reliable depends on the sequence data. I knew that my data set showed similar result topology from Iq2 tree output.
There is no species sequence data available for the species that I am looking for a reliable topology of those species.
If anyone knows how to understand the assigned sequence on the unrooted tree based on Four-cluster Likelihood-Mapping by using iqtree2? Could you kindly answer this question? I want to know the orange unrooted tree (image 2) for my analyzed data. Thank you very much.
output ***.iqtree
Quartet support of regions a1, a2, a3 (mainly for clustered analysis):
     #quartets    a1 (% a1)        a2 (% a2)        a3 (% a3)    name
-----------------------------------------------------------------------------
1   2000     648     (32.40 ) 712     (35.60 ) 640     (32.00 )  Haur_2
2   2000     648     (32.40 ) 712     (35.60 ) 640     (32.00 )  Shad_6
3   2000     648     (32.40 ) 712     (35.60 ) 640     (32.00 )  Equa_oga2
4   2000     648     (32.40 ) 712     (35.60 ) 640     (32.00 )  apa
-----------------------------------------------------------------------------
2000     648     (32.40 ) 712     (35.60 ) 640     (32.00 )

You have provided two different trees. The second was the Arthropoda and labeled "image 2".
The first example the explanation is that it is an exact polytomy. I know this because the values are equally split between the 3 possible alternatives for a 4-taxon tree. Thus there is no bifurcated solution to the data, any bifurcating combination is equally valid therefore it is simply expressed as all 4-taxa sharing the same node (hypothetical common ancestor). In other words there are 3 equally valid solutions.
Thus, the tree is 4-taxa tree with 1-node. Thats the solution. A normal 4-taxon tree will have 2-nodes, this is a star-tree, i.e. one node, technically a polytomy.
The example below concerns Arthropoda and demonstrates the top tree ((Chelicerta, Myriapoda),(Hexapoda, misc.)) is the correct tree. That 90.4 value demonstrates that. What this says is that neither myriapods or Chelicerta share a direct common ancestor with hexapods, moreover "misc" (whatever that is) shares a direct ancestor with the hexapods. "Image 2" is a classical bifurcating tree and represents only one solution.

So to specifically answer the question. The tree for arthropoda is:

and the other tree is,

I should have removed the scale bar in this example, because I didn't have the data, but not when you have genetic quantification of the distances.

Thank you very much for your answer. I have searched about polytomy. It is right that I used the data from a single ancestral lineage. But then, I do not understand the recent paper Nicolás Mongiardino Koch, et al. eLife (2022) showed a likelihood-mapping analysis in IQ-TREE to visualize the phylogenetic signal for alternative resolutions of the quartet involving these three lineages (Arbacioida, Salenioida, and Stomopneustoida + Camarodonta) and their sister clade (Irregularia). My data was also sister clade. Why my analyzed result won't be like Koch, et al. eLife (2022) Fig2. B?

@RioKashimoto this is a separate question. I would need to go through their paper in detail and would take a while. Their analysis is not necessarily correct. The paper is here https://elifesciences.org/articles/72460 ... I have been through the paper briefly and there could be areas of concern.


text:
CRISPR genome wide screen analysis using MAGeCK-- Loss of reads/sgRNA using MAGeCK count function?
Sorry if this question seems basic. I am trying to run the MAGeCK pipeline to analyze CRISPR knockout screen data produced by someone in my lab around 5 years ago. I was given the data as bam files and also have the sequencing sample statistics from when the samples were run including the number of total reads, aligned reads, and unaligned reads. I converted the bam files back to fastq files using
samtools view -h -F 2048 filename.bam > tmp.bam

bedtools bamtofastq -i tmp.bam -fq filename.fastq

After converting the bam files back to fastq files, I ran the MAGeCK count function. I saw that none of the total reads matched between the MAGeCK count summary output file and the original sequencing statistics.
More problematically, a few samples had significantly fewer total reads (like 50,000) shown on the MAGeCK count summary output files relative to the original sequencing statistics (~44 million). As a result, many of the sgRNAs are not represented at all in the sample (there are ~61,000 zero counts relative to the total 77441 sgRNAs) which is impacting my analysis.
Can anyone help me understand why I might be losing sgRNAs using MAGeCK? Is it a problem with my source files? I would appreciate any help!

It's hard to say without looking more closely at the data, but you may be able to find a few guides that are drastically different in your counting, and see how prevalent they are within the original BAM file (before filtering). Perhaps the -F 2048 being used to filter out chimeric reads is overfiltering, and causing the discrepancy? Are unaligned reads included in the BAM - if not, does this explain the discrepancy?
Since the total reads doesn't match the original, I'd look into this. Are the extracted FASTQs consistent read length? Maybe some trimming was performed that removed the guides?
Do you know what reference the BAM files were mapped to? Perhaps this reference did not include the guide sequences, and were mapped to the genome instead? You may be able to extract this information from the BAM like this:
 samtools view -H filename.bam > filename.header

Ideally you'd have the original FASTQs to work with, but it sounds like those are unavailable...


text:
DNA, molecular descriptors
Are there any programs (preferably command line tools) for calculating molecular descriptors of DNA? I am looking for something like Chemopy or RDKit but specifically for DNA.
Thanks in advance!

I would like to perform classification after all, so the more descriptors are the better for me.

Please [edit] your question and add some details. Explain exactly what "descriptors" you are thinking of, tell us what the input to this program would be and also tell us what operating system you are using.

Hi @BalintBiro, thanks for your feedback and answering your own question. Could you kindly mark this as "accepted", it helps with site stats.

Just found the right solution.
iFeatureOmega does exactly what I was looking for. Sequence based and has CLI.
https://ifeatureomega.erc.monash.edu/iFeature2/index.html
Additionally it supports protein, DNA/RNA and also small molecule inputs.

Thanks for your answer on iFeatureOmega, cool. Could you please mark the answer as "accepted"?


text:
Can I run dotPlotly from Rstudio?
I'd like to run dotPlotly from RStudio but I can't find any description/help for this. Is it at all possible?
Another side to this question: does anyone know any tool for dot plots for genome comparison apart from pafr (doesn't provide colour legend) and d-genies (we have problems with installation)? Is there a way to make dot plots for a minimap2 output with ggplot2?
To specify:

I need help in running dotPlotly (preferably from Rstudio)
I'd like to build dot plots for genome comparisons showing % of identity

Thanks

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

For the issue of just making a dotplot in ggplot2 from a PAF file, it is pretty trivial to do a simple one, though it is necessarily not the most beautiful.
I grabbed 2 random E. coli genomes and mapped them:
% minimap2 ~/Downloads/GCF_000157115.2_Escherichia_sp_3_2_53FAA_V2_genomic.fna.gz ~/Downloads/GCF_000308975.1_ASM30897v2_genomic.fna.gz > 2_ecolis.paf

Then I opened Rstudio and ran the following (which is only slightly modified from what I was first shown by an old colleague, Zev Kronenberg):
ecoli_paf = read.table("~/2_ecolis.paf", sep="\t", fill=T)
ggplot() +
geom_segment(
  data=ecoli_paf, mapping=aes(
    x=ecoli_paf$V3, y=ecoli_paf$V8, xend=ecoli_paf$V4, yend=ecoli_paf$V9, color=ecoli_paf$V6
  )
)

I am far from a ggplot2 expert as you can tell, and there are a few issues here (reverse complement is shown opposite to what you'd expect, linear interpolation of mappings from mapped endpoints), but it makes a somewhat informative dotplot:

I am confident that with minimal fiddling you can make this look a lot better. I will try to revisit and beautify.

Thanks a lot, this is useful! My problem is that I'd also need some reflection of similarity, like the % of identity. Specific tools provide it. Is there a way to show it in this simple ggplot approach?

Thank you @SKay.

@SKay yes there should be, you would have to use another layer like perhaps line thickness or transparency, or mapping colors to that instead and tracking scaffold id differently. PAF output from minimap2 captures similarity so everything is there. If i get a moment to improve my answer I will try a few approaches and see what seems best.

Thank you very much!


text:
Visualise multiple disappearing pdb files in Pymol
I have a sequence of pdb files which I want to load in pymol window so that it gives an impression of a video. I am using cmd.load function in a loop and deleting the loaded file in the next step. But for some reason the window pauses until all iteration steps have been completed and giving nothing else in the end. I know that the files are loaded because without the “delete” line I can find all pdbs overlapped on one another. Any suggestions? Thanks and regards
def trial(pathname, num_files):
    for i in range(num_files):
        filename = (pathname + str(i + 1)  +'.pdb')
        print(filename)
        # print(i)
        cmd.load(filename, (i+1))
        cmd.hide()

https://pymolwiki.org/index.php/Frame is the way —you are after moviemaking which is rich in tutorials, so I am just leaving the link here. Hopefully a tutorial about working with Frames with answer your Q.


text:
How do I generate a variant list using Illumina reads from a Salmonella genome?
I am planning on performing phylogenetic analysis of Salmonella specimens using WGS data from PulseNet and GenomeTrakr, for the purpose of public health surveillance & to provide context for foodborne outbreaks.
I have been using the Nextstrain workflow to generate and visualize trees for SARS-CoV-2 specimens, and I am planning on adapting the bacterial workflow outlined here for my bacterial isolates: https://docs.nextstrain.org/en/latest/tutorials/creating-a-bacterial-pathogen-workflow.html
What approach would you use to process my reads into VCF files for phylogenetic analysis? I've been referencing information from this similar post from 5 years ago: How do I generate a variant list (i.e. VCF file) using Illumina reads from a human genome?
How much of this information is still current/relevant and applicable to Salmonella isolates?

I would express caution about cross-over because for both SARS-CoV-2 and the human genome you can obtain a good assembly from a reference genome. The SARS-CoV-2 genome is so small and read depth high enough this is not needed  however. SARS-CoV-2 there are critical issues in assemblies for certain regions which have not been addressed, but for CDS there's no real difficulties here.
A reference genome assembly is  not good practice for what are often environmental bacteria (bacteria per se) and de novo assembly is strongly recommended, even if it is to complement a reference genome assembly. Bacterial evolution is very different from either human or virus evolution, it's full of indels, missing genes, genes 'flying' around on plasmids and potential non-homologous recombination, risk of rearrangements etc ...
The challenge you face is automating the assemblies and then polishing those assemblies and there has been enormous progress recently, "closing the gaps" manually is a thing of the past - ancient history. There should be no under-estimation as the volume of recent work that has been dedicated particularly towards long read data. Even pre-assembly, pre-polishing the progress shouldn't be underestimated, albeit I suspect you're familiar with that. The old ideas of QUAST are now historic. The raised ante is noticeable.
Back to the original question consistently using a reference genome risks something going wrong that could cause an outlier artefact in the VCF and in my opinion its no longer cool for bacteria.
The rationale is clear, mis-assembly  would cause an artefact in the tree and therefore in the epidemiological inference of the outbreak.
It's worth doing specific background reading because albeit I've deliberately not mentioned a given approach except for QUAST, the rate of progress has been significant. The world has changed.


text:
Questions about ChIP-seq reads and MACS
I had a few questions raised by the below diagram in relation to how ChIP-seq reads look like (source). I would be very grateful for your insight on them!
From the below diagram, it seems that, in a ChIP-seq experiment, we get information about the binding site in relation to both the Watson and Crick strands. That is, from our experiment, we get reads complimentary to both strands.
(Q1) I suppose this is because proteins bind to both strands of DNA, so that in our ChIP-seq experiment, we ultimately sequence reads from both the Watson and the Crick strands? Is that true?
Also, my impression from the MACS paper was that "ChIP-Seq tags represent only the ends of the ChIP fragments, instead of precise protein-DNA binding site" (Background).
(Q2) Therefore, shouldn't it be that we get reads mapping to either side of the binding site on BOTH strands? That is, we get reads more towards the 3' end of the Watson strand, and reads more towards the 5' end of the Watson strand, and similarly for the Crick strand.
(Q3) The diagram below shows we only get reads more towards the 5' end of the strands. Why is that?
(Q4) Finally, I know that MACs shifts reads by some bases towards the 3' end for each strand, in order to interpret the precise binding site location. It seems to assume what I am asking about in (Q3). However, if the answer to my (Q2) is "yes, we do get 4 such clusters of reads", then why doesn't MACs account for that?


text:
How can I align two novel sequences to a reference genome and build a phylogenetic tree?
I have two new amino acid gene sequences that I  want to align to a reference genome. Is it a good idea to have a combined phylogenetic tree? what will this type phylogenetic analysis called?

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.

What is the actual question you are trying to answer here? Are you trying to identify the gene coding for these proteins? Is the genome annotated? Are you working on a species that has introns? Alternative splicing?

Edit: just to point out a tree needs 4 taxa (sequences), so rather than one reference genome, 3 reference genomes. Thus for Blast you might need to produce a database (using 3 reference genomes) before Blasting. There are other Blast approaches to obtaining a large meaningful sample size ... thats a separate question

It called a combined alignment, or concatenated genes.
The theory is simply that if two genes evolve at noticeably different rates they can't be combined into a single alignment and make a tree from them. The theory goes back a long but the modern test is in IQTREE.
You take your two alignments and test them for equal rates via IQTREE
iqtree2 -s example.phy -p example.nex --symtest-only

Testing for a combined alignment is also described in my past answer here: The previous question was can you combined an rRNA gene with a protein gene and the answer is ...

Okay, you can't concatenate genes into a singular matrix particularly if these are rRNA genes versus protein genes, except in rare situations where the tests say 'thats okay' (below). This is simply because they evolve at very different rates rRNA genes have loci which evolve very slowly and are highly conserved.

The first part of the question was mapping the gene to a reference genome. Anything will do that ... local Blast works and using that is described in my answer to this question here:
If you loads of fasta seqs you'll need to make a database...
makeblastdb -in example700-1000.fa -parse_seqids -blastdb_version 5  -out laura700-1000.db -title "Laura 700 - 1000" -dbtype nucl

Then do this ...
blastn -query laura.fas -db laura700-1000.db -evalue 1e-6 -num_threads 4 -out blasthits.txt

However if its just reference vs. query then you can just use the above command.

The alignment is best done via a specialist alignment package Muscle Super5 and that is described in my answer here. This is done after your obtained the homologous regions in your reference genomes.

HMM Super5 is a hidden Markov model (HMM) where gap penalties and the substitution matrix are the parameters at the centre of the model. That is the answer, thats why you can't change it - it does it itself and thats its purpose. Thus one algorithm (super5) is HMM for the gap-penalty and the other one I suspect isn't.

Can you really get meaningful results using just two proteins, each aligned to its own reference genome, as input? Surely any tree built from such limited data will be unreliable, right? Why would you even build a tree if all you have are two proteins? Also, blast really isn't a good tool for this unless we know the OP is working on species with no introns, and in any case, they would need tblastn if they start from proteins.

@terdon a tree needs 4-taxa ... so 3 reference genomes and the sequence of interest. Low reliability ... no ancient history, most tree methods moved on, its not affected by the sample size - absolutely not. Blast, despite being slow, is the work behind most published phylogenies.  .. introns .. it doesn't matter most quality tree building methods (IQTREE for example), ignore indels (introns are excluded) - why is a separate question -- thats a long explanation -- but based on their mutation matrix (indels are absent).

Edit comment above ... introns are excluded as phylogenetic information if they result in indels, things like Beast, IQTREE, RAxML all ignore indels. Thus even if they is a single indel for 10000 sequences that entire column is excluded for that particular nucleotide/amino acid.

OK, but how would you even build a tree using just two alignments, of two sequences each? As you said you need at least 4 taxa, and if (big if) I understand the question correctly, the OP has at best 2 species (no idea if they're the same taxon or not) and at worse just two genes from the same species.

@terdon I edited the posted previously ...

Brilliant, thanks!

I would strongly recommend to find an outgroup or at least a fourth sequence. Building a phylogenetic tree with 3 sequences is doable but is tricky and error prone.
If you still want to do a 3-sequence analysis, you can align the two homologous proteins to the reference genome with exonerate, spaln or miniprot, extract the reference subsequence in the alignment, do a proper multiple-sequence alignment (MSA) with muscle or MAFFT, build a tree and put the root at the middle point of the longest branch (this assumes a molecular clock). After that, you can do typical phylogenetic analysis.
Whether the above works depends on the divergence of your sequences. If they are closely related and one sequence is a clear outgroup to the other two, the procedure could give you something, though having more sequences is still recommended.
EDIT in response to the comments: midpoint rooting assumes an approximate molecular clock across the whole tree. Severe violation of the assumption may lead to a wrong root. This approach is not very reliable due to the clock assumption.

Definitely agreed on 4 sequences. IMO an unrooted tree obviates the need to assign an outgroup. BTW mid-point rooting doesn't assume a global molecular clock, most tree building produce phylograms - non-clocks (Beast is an exception, options depending). UPGMA method assumes a clock.

@M__ I am repeating myself, so I don't think this discussion goes anywhere. I have deleted my comments and summarized them in the answer. Feel free to disagree but I am confident in my knowledge as well.

Thanks @user172818 I've done likewise.


text:
Nextflow with Minimap2 not working
I'm trying to create a Nextflow script that maps assemblies to reference gene sequences. This is what I have so far:

#!/usr/bin/env nexflow
nextflow.enable.dsl = 1

params.assmblreads = '/Users/rimo/Desktop/Assemblies/Assemblies_Flye/flye_assembly.fasta'
params.assmblref = '/Users/rimo/Desktop/Assemblies/gene_targets/GENE.fsa'

refgenome_file = file(params.assmblref)
fasta_files = Channel.fromPath(params.assmblreads)
                     .map { file -> tuple(file.baseName, file) }

process minimapping {

    input:
    file genome from refgenome_file
    set val(readID), file(readFile) from fasta_files

    output:
    set val(readID), file("${readID}.paf") into mapped_reads

    script:
    """
    ./minimap2 $readFile $genome > ${readID}.paf
    """

}

I'm running this on a macOS so I have to include the "./" before a command or it will give the error: command not found"
But everytime I run the script with the command: ./nextflow run trials.nf it results in this error message:
(base) rimo@Ritas-MBP minimap2 % ./nextflow run ../attempt_out.nf
N E X T F L O W  ~  version 22.10.0
Launching `../attempt_out.nf` [nasty_baekeland] DSL1 - revision: 98d1deaebd
executor >  local (1)
[ae/47e3c0] process > minimapping (1) [  0%] 0 of 1
Error executing process > 'minimapping (1)'

Caused by:
executor >  local (1)
[ae/47e3c0] process > minimapping (1) [100%] 1 of 1, failed: 1 ✘
Error executing process > 'minimapping (1)'

Caused by:
  Process `minimapping (1)` terminated with an error exit status (1)

Command executed:

  ./minimap2 flye_assembly.fasta GENE.fsa > flye_assembly.paf

Command exit status:
  1

Command output:
  (empty)

Command error:
  .command.sh: line 2: ./minimap2: No such file or directory

Work dir:
  /Users/rimo/Desktop/AWS/minimap2/work/ae/47e3c0708f33b1949315fc83066b95

Tip: you can try to figure out what's wrong by changing to the process work dir and showing the script file named `.command.sh`

This is my first time ever writing a Nextflow script so I'm very new and confused to this new way of scripting. If someone could help me out that would be great!

Command error:
  .command.sh: line 2: ./minimap2: No such file or directory

Work dir:
  /Users/rimo/Desktop/AWS/minimap2/work/ae/47e3c0708f33b1949315fc83066b95

The problem is that your command is expecting the minimap2 binary inside your working directory (i.e. the 'work dir', as defined above), but it could not be found. To avoid using relative (or absolute) paths, and to help keep your workflow portable, instead move your minimap2 binary into a directory in your $PATH. Then, you should be able to call minimap2 without a relative (or absolute) path to it. Here's one way using the new DSL2, which comes enabled by default with Nextflow 22.10.0:
params.assmblreads = '/Users/rimo/Desktop/Assemblies/Assemblies_Flye/flye_assembly.fasta'
params.assmblref = '/Users/rimo/Desktop/Assemblies/gene_targets/GENE.fsa'

process minimap2 {

    input:
    tuple val(readID), path(readFile)
    path genome

    output:
    tuple val(readID), path("${readID}.paf")

    """
    minimap2 "${readFile}" "${genome}" > "${readID}.paf"
    """
}

workflow {

    refgenome_file = file( params.assmblref )
    fasta_files = Channel.fromPath( params.assmblreads )
                         .map { file -> tuple( file.baseName, file ) }

    minimap2( fasta_files, refgenome_file )
}

To make your pipeline more portable, and to avoid manually installing and configuring your own tools, another way is to use Conda environments. Nextflow has built-in support for Conda and gives you access to all of the tools, including Minimap2, available in Bioconda. Install miniconda to get started, then add conda.enabled = true to your nextflow.config and the conda directive (e.g. conda 'minimap2=2.24') to your process definition for example. In the same way, you could also use Docker (or Singularity) containers with the container directive. Note that Docker containers are also built by the Bioconda build system and these are uploaded to quay.io as part of the Biocontainers project (e.g. biocontainers/minimap2).


text:
converting coordinates to sequences using bedtools getfasta - [segmentation fault] errors
I use bedtools to convert genomic coordinates (in the form of bed files, chrX x1:x2) to genomic sequence content in the form of fasta files ('ACGT...G').
The bedtools format for extracting sequence (.fa file) from coordinates (.bed file) is:
bedtools getfasta -fi  input_sequences.fa -bed  genomic_coordinates.bed -fo output_sequences.fa

Very frequently I get a "Segmentation fault (cored dumped)" during the run.
Warning: the index file is older than the FASTA file.
Segmentation fault (core dumped)

The strange thing is that the error will go away, if for example I add additional data (chromosomes) to the input fasta file (input_sequences.fa), OR remove some chromosomes?? I am having trouble finding a logical path to explain this behavior so I can work around it.
There is no changing of fasta files while bedtools is working - in fact there is nothing in the primary directory other than the two files, the fasta file input_sequence.fa and the bed file genomic_coordinates.bed. There is 256GB of RAM memory on the machine, but that's not an issue, especially since the error goes away sometimes when additional chromosomes added to the fasta file.
Anyone experienced issues with bedtools getfasta erratic Segmentation faults? Known workarounds? I am using the current version of bedtools v2.30.0
What are some of the alternatives to "bedtools getfasta" for converting genomic coordinates to sequences?
Thanks...

This sounds like C thing. I don't do eukaryotes (mostly). I am aware of the strong reputation of bedtools, but thats exactly what C does if the data, code and RAM don't mix

Thanks @M__  you are right, the error is likely caused by bedtools accessing (inside it's C/C++ code) un-allocated memory space. It's possible an array get's defined under certain conditions where the expected/allocated size ends up not being correct - often it's a programming oversight - what I am trying to understand is the condition that may be leading up to it, so I can avoid the condition: something related to makeup of the input fasta file as related to the coordinates in the bed file.

Have you tried taking a sample of your bed? Say, for example, splitting it into 3 parts and running the code again. Try to narrow down what is causing that error. Is the bed a big file? Maybe, you don't have enough RAM?

The error about old index files is insightful to me. Bedtools presumably assumes that the input files are static / unchanging. If bedtools is automatically indexing the files before extracting sequences (when index files don't already exist), and the input files keep changing, that would be a good explanation for why problems are happening.
Do these input fasta files have bedtools and/or samtools index files in the same directory (usually fileName.fa.fai). Are you re-indexing the fasta files after changing them (e.g. via samtools faidx fileName.fa)?
I don't mean changing during the running of bedtools, I mean changing between different bedtools runs. If you create a file input_sequences.fa, then the first time it is subset, it will create an index file input_sequences.fa.fai. The next time bedtools is run to subset an input_sequences.fa within the same directory, it will check to see if an index file already exists; if one exists, it won't be regenerated. If the file input_sequences.fa has changed in the interim, then the index file will be incorrect, and the results of subsetting will be unpredictable.
This is why Bedtools reports the error: it can cause strange things to happen with output (intermittent segmentation faults, for example).
In fact, bedtools is explicit about this index generation in its output:
$ ls
sampleTagsNew.bed  sampleTagsNew.fa

$ bedtools getfasta -fi sampleTagsNew.fa -bed sampleTagsNew.bed
index file sampleTagsNew.fa.fai not found, generating...
>NewSampleTagX01:0-20
GTTGTCAAGATGCTACCGTT
>NewSampleTagX02:1-20
TTGTCAAGATGCTACCGTT

$ ls
sampleTagsNew.bed  sampleTagsNew.fa  sampleTagsNew.fa.fai

$ perl -i -pe 's/NewSampleTagX01/NewSampleTagX01_updated/' sampleTagsNew.fa

$ bedtools getfasta -fi sampleTagsNew.fa -bed sampleTagsNew.bed
Warning: the index file is older than the FASTA file.
>NewSampleTagX01:0-20
updatedGTTGTCAAGATG
>NewSampleTagX02:1-20
TagX02GTTGTCAAGATG

$ rm sampleTagsNew.fa.fai

$ bedtools getfasta -fi sampleTagsNew.fa -bed sampleTagsNew.bed
index file sampleTagsNew.fa.fai not found, generating...
WARNING. chromosome (NewSampleTagX01) was not found in the FASTA file. Skipping.
>NewSampleTagX02:1-20
TTGTCAAGATGCTACCGTT

In the example above, I carry out the following operations:

List the files in the directory
Run bedtools getfasta to subset a fasta file. Bedtools reports that an index file did not exist, so it is generating a new one.
List the files in the directory. Note that there is a new file: sampleTagsNew.fa.fai, which is the fasta file index.
Modify the sample tags fasta file [to change the name of a sequence]
Run bedtools getfasta again with an identical command line. It produces the same warning you have seen, and produces unexpected output, i.e. incorrect results.
Remove the fasta index file
Run bedtools getfasta again with an identical command line. This time the output is correct, stating that the first sequence cannot be found in the input fasta file.

Thank you @gringer... update is SUPER helpful.  Indeed the fasta file contents do change between subsequent bedtools calls, while the input fasta file name stays the same. There IS an *.fai file which I deleted and reran the bedtools call that was consistently causing segmentation fault. This time I got no segmentation fault!!! As you suspected index file *.fai seems to be the issue, but not all outdated *.fai files cause problems. I  just tested that by mixing and matching  all kinds of outdated *.fai files, no segmentation fault.  For now I am happy to manually delete *.fai's. Thanks again!!


text:
What does it mean when a gene and its transcript have opposite orientations in a GFF3 file?
I was working with a given GFF3 file, and I observed that some transcripts have orientation opposite to their transcripts. Here is a snippet:
chr1    .       gene    2189548 2194772 .       +       .       ID=c5eaf8c1-2bb2-566f-e154-6a9be797c11f;Name=FAAP20;gene_biotype=protein_coding;gene_id=34a96673-bfab-ebbf-b2b7-f9cb73163a4a;gene_name=FAAP20
chr1    .       transcript      2189548 2194772 .       -       .       ID=f47feeac-74ce-1543-53b6-afd9bbb8ab19;Name=rna-NM_182533.4;gene_biotype=protein_coding;gene_id=34a96673-bfab-ebbf-b2b7-f9cb73163a4a;gene_name=FAAP20;Parent=c5eaf8c1-2bb2-566f-e154-6a9be797c11f;transcript_biotype=protein_coding;transcript_name=rna-NM_182533.4;protein_id=NP_872339.3

How do I interpret this? I checked UCSC Browser and in each case it seems that the browser agrees with the transcript in each case, which seems like it should be canonical.
Moreover, I observed in this file that gene orientation is always +. Whereas a highly curated GFF3 of the human genome represents both directions:
# my file
$ awk '$3 == "gene" {print}' filtered_potential_genes.gff | cut -f7 | sort | uniq -c
   1600 +
# highly curated annotations
$ awk '$3 == "gene" {print}' /MANE.GRCh38.v0.95.select_refseq_genomic.gff | cut -f7 | sort|  uniq -c
   9393 +
   9191 -

According to the person who provided me the file, the existence of bidirectional transcription was indeed the motivation for this. However, it's obviously confusing to have it still be a definitive + rather than . or something to signify unstrandedness.
It seems likely that this is due to an artifact in my file. Are there other interpretations? Is the standard simply that gene orientation is meaningless and transcript orientation is important?

Looks like a poor annotation to me. I could understand a few scattered ncRNAs being arranged in the opposite direction from the gene, but all genes being oriented identically doesn't make sense.