text:
How can I bypass 'out of memory' problem using a HISAT2 for human genome indexing?
I am trying to perform a HISAT2 indexing using a GCF_000001405.39_GRCh38.p13_genomic.gtf and GCF_000001405.39_GRCh38.p13_genomic.fna files downloaded from https://www.ncbi.nlm.nih.gov/assembly/GCF_000001405.39 using the procedure described in  Pertea et al. (https://www.ncbi.nlm.nih.gov/pubmed/27560171), e.g. first I create extract_splice_sites.py and extract_exons.py and then my code looks like this
#!/bin/bash
#SBATCH --time=02:30:00
#SBATCH --account=def-myusername
#SBATCH --mem=160000M
#SBATCH --mail-user=mymail@gmail.com
#SBATCH --mail-type=All

module load nixpkgs/16.09  intel/2018.3
module load hisat2/2.1.0

hisat2-build --ss imam.ss --exon imam.exon GCF_000001405.39_GRCh38.p13_genomic.fna imam002

I am trying to receive a set of imam002*.ht2 files but every time the procedure ends with a Ran out of memory message. From here https://github.com/griffithlab/rnaseq_tutorial/wiki/Indexing I learned that I must allocate 160Gb of RAM if I want to index the entire genome. Even when I do that and submit my script via sbatch on a HPC I still got this
slurmstepd: error: Detected 1 oom-kill event(s) in step 30827017.batch cgroup.

Just in case here's the window that I have when I don't sbatch my script, I still have an error
Inactive Modules:
  1) openmpi/2.1.1

The following have been reloaded with a version change:
  1) gcccore/.5.4.0 => gcccore/.7.3.0     2) icc/.2016.4.258 => icc/.2018.3.222     3) ifort/.2016.4.258 => ifort/.2018.3.222     4) imkl/11.3.4.258 => imkl/2018.3.222     5) intel/2016.4 => intel/2018.3

Settings:
  Output files: "imam002.*.ht2"
  Line rate: 7 (line is 128 bytes)
  Lines per side: 1 (side is 128 bytes)
  Offset rate: 4 (one in 16)
  FTable chars: 10
  Strings: unpacked
  Local offset rate: 3 (one in 8)
  Local fTable chars: 6
  Local sequence length: 57344
  Local sequence overlap between two consecutive indexes: 1024
  Endianness: little
  Actual local endianness: little
  Sanity checking: disabled
  Assertions: disabled
  Random seed: 0
  Sizeofs: void*:8, int:4, long:8, size_t:8
Input files DNA, FASTA:
  GCF_000001405.39_GRCh38.p13_genomic.fna
Reading reference sizes
  Time reading reference sizes: 00:00:23
Calculating joined length
Writing header
Reserving space for joined string
Joining reference sequences
  Time to join reference sequences: 00:00:13
  Time to read SNPs and splice sites: 00:00:15
  Ran out of memory; automatically trying more memory-economical parameters.
./hisat2_index.sh: line 24:   943 Killed                  hisat2-build --ss imam.ss --exon imam.exon GCF_000001405.39_GRCh38.p13_genomic.fna imam002

When I use the same script for another genome (not human) - I don't have this problem.
Is it really a problem of RAM and, if it is, how is usually human genome gets indexed by HISAT2? Is it even possible?

Perfectly possible, only need more RAM. Set it higher than 160G then, SLURM does not like it when processes reach the maximum allocated RAM. You can also download indices from the hisat2 website: https://ccb.jhu.edu/software/hisat2/manual.shtml

You haven't allocated the requisite 160GB of RAM required. Or rather, you allocated 160GB of memory, which is what you were told to do, but they most likely meant 160GiB.
Computers use binary, not decimal, so a gigabyte is 1024 megabytes. This is the same annoying mix-up people run into when their hard drives don't have the space the marketing department says it does. By requesting 160000M in your job script, you actually requested a few gigabytes less than what your source says you need.
$$
160000 \text{ megabytes} \times \frac{1 \text{ gigabyte}}{1024 \text{ megabytes}} = 156.25 \text{ gigabytes}
$$
$$
160 \text{ gigabytes} \times \frac{1024 \text{ megabytes}}{1 \text{ gigabyte}} = 163840 \text{ megabytes}
$$
If "160GB" is what you need, then 163840M is the correct value for the job script.
I would honestly suggest cutting the Gordian knot and not dealing with this "gigabytes vs. gibibytes" ordeal by simply using the G suffix.1
#SBATCH --mem=160G


text:
Primer-Blast Custom Database [error]
I'm using Primer-Blast to compare a PCR template against a custom database of reference sequences, using Primer-Blast's Custom database functionality.
However whenever I try to upload a FASTA file as the custom database and then run Primer-Blast, I get the following error:
Exception error: Fasta Reader: sequence id ends with 80 valid nucleotide characters. Was the sequence accidentally placed in the definition line?
No matter how I try to format my FASTA file, I get this error.
Does anyone know how to properly format a FASTA file to have it accepted as a custom database of reference sequences to compare against the PCR template? Or is this a bug in the software?
Here's an example FASTA file, to upload both for the template and the custom database.
>test1
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAACGAACTTTAA
AATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAACTAATTACTGTCGTTGACAGG
ACACGAGTAACTCGTCTATCTTCTGCAGGCTGCTTACGGTTTCGTCCGTGTTGCAGCCGATCATCAGCACATCTAGGTTT
CGTCCGGGTGTGACCGAAAGGTAAGATGGAGAGCCTTGTCCCTGGTTTCAACGAGAAAACACACGTCCAACTCAGTTTGC

Could you provide a sample of your FASTA file, pls?

@mrhd just added

@mrhd It's almost like it's not recognizing the newline escape character at the end of the sequence ID.

That might actually be the cause of the error. Did you create the FASTA on Windows or Linux/Mac?

@mrhd The actual FASTAs I want to use were downloaded from GISAID. However I am on a Mac.

Ok so it looks like the workaround for this is to copy-paste the template into the text box, and then upload your custom database. If you upload both the template and the custom database, you will get the error.

I just copied your FASTA sequence into a file (on Mac OS) and uploaded it to Primer BLAST, the search did not give an error.
If the error is caused by the line endings, maybe you can fix your downloaded files with a tool like dos2unix.
Update:
On a whim, I tried adding a description to the FASTA header:
>test1 description
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAACGAACTTTAA
AATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAACTAATTACTGTCGTTGACAGG
ACACGAGTAACTCGTCTATCTTCTGCAGGCTGCTTACGGTTTCGTCCGTGTTGCAGCCGATCATCAGCACATCTAGGTTT
CGTCCGGGTGTGACCGAAAGGTAAGATGGAGAGCCTTGTCCCTGGTTTCAACGAGAAAACACACGTCCAACTCAGTTTG

This time, there was no error.
Does that work for you as well?

You uploaded it both as the PCR template and as the custom database? Because it works fine when you just upload it as the template, and then set the database to something like Refseq representative genomes. It's only when you set the database to custom that it errors out.

Sorry, you're right. When uploading the sequence as database, I get the same error. -.-

The error does not appear when pasting the sequence into the text box under "Custom", though. Strange...

Ok, for some crazy reason, I'm able to get it to work when I copy-paste the template but upload the custom database. But it doesn't work when you upload both the template and the custom database, even though it does work if you upload the template but don't use the custom database option.

Very odd indeed..

I did it by copy and paste. Here's the top line of the output:
    Primer pair 1
Sequence (5'->3')   Template strand Length  Start   Stop    Tm  GC% Self complementarity    Self 3' complementarity
Forward primer  TGTCGTTGACAGGACACGAG    Plus    20  148 167 59.97   55.00   5.00    2.00
Reverse primer  TTACCTTTCGGTCACACCCG    Minus   20  264 245 59.97   55.00   3.00    2.00
Product length  117
Primer pair 2
Sequence (5'->3')   Template strand Length  Start   Stop    Tm  GC% Self complementarity    Self 3' complementarity
Forward primer  ATCTGTGTGGCTGTCACTCG    Plus    20  82  101 60.04   55.00   5.00    3.00
Reverse primer  ACTCGTGTCCTGTCAACGAC    Minus   20  168 149 59.97   55.00   4.00    3.00
Product length  87
Primer pair 3
Sequence (5'->3')   Template strand Length  Start   Stop    Tm  GC% Self complementarity    Self 3' complementarity
Forward primer  TACTGTCGTTGACAGGACACG   Plus    21  145 165 60.00   52.38   7.00    3.00
Reverse primer  ACACCCGGACGAAACCTAGA    Minus   20  251 232 60.54   55.00   4.00    2.00
Product length  107
Primer pair 4
Sequence (5'->3')   Template strand Length  Start   Stop    Tm  GC% Self complementarity    Self 3' complementarity
Forward primer  GCTTACGGTTTCGTCCGTGT    Plus    20  192 211 60.94   55.00   4.00    0.00
Reverse primer  GTTGAAACCAGGGACAAGGC    Minus   20  292 273 59.33   55.00   4.00    2.00
Product length  101

A command line solution would work, based on the bug reported by the OP. However there's nothing in NCBIs e-utilities and Biopython doesn't do Primer-Blast, but will do Blast. You could do it Perl via the LWP Agent, e.g. here.

source: Stack Exchange

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

sort: votes


text:
Methods to predict one connected structure from two known separate structures
Currently, I am facing such a question: how to predict one connected protein structures from two separate known structures?
It should be different from the protein-protein docking method, which is non-covalent interactions. My question is that the two proteins should "contently" link together.
Any suggestions are appreciated.

MODELLER can do this as well.

There are two approaches here.
Physiological
If domain A and B of your protein interact with each other physiologically, then you need to protein-protein docking. This is not too straightforward as protein are flexible.
Several online servers exists such as Haddock and include ClusPro, which use more traditional approaches. While the online server DEMO (I-Tasser group) does a structural homology based search and a final MD based repacking to better minimise the structures (no protein-protein docking).
Many of these work by magic as far as the user is concerned. But understanding what is going on is really important, even if it is overly laborious. When you don't know where a protein goes, it's called blind or global docking, which is computational expensive and so imprecise that the effort required to prove it is correct for peer-review is too steep for it to be useful. If you are doing local docking it is okay, which may be the case here depending on the linker. Specifically, you need to generate more 'decoys' (potential poses) the more unsure you are: tens of thousands or more. Then you have find clusters and see how each scores following more precise refinement. As an example of a tediously and manually set up run, Rosetta has several docking protocols, which explain the steps required —of note are the many judgement calls, such as in the local refinement, where one has to judge whether setting -dock_pert 3 8 (3 Å and 8° perturbation) is too strict etc.
Pretty picture
The protein-protein docking approach will result in a globular protein, which may be incorrect and majorly is not great for pictures.
If you want a pretty picture of two or more domains most likely you would want them placed like beads on a string. If you have a linker that is not solved the recommended distance to space them is $(3.5 N_{missing})^{0.5}$ Å —unfortunately I do not know where this lore comes from. If it helps I wrote a script that does this in automatic but manual placement works best.
Linkers
In the situation you have a good confidence about how the two parts fit together (or have placed them like beads on a string), you may want to add linkers. If you absolutely don't care about correct bond topology here is a blog post on how to join bits (shoddily) with PyMOL and here is one to do it properly with Rosetta or Pyrosetta.
If one wanted to do a better job, one tool to use is Rosetta remodel, which has a bit of a learning curve, but is very powerful for protein design. This requires the PDB to be on a single chain (pymol alter command, followed by create will do the trick: more here) and actually the amino acids to not have gaps —PDB number being the same as the pose numbering. Rosetta gives better results if the poses are energy minimised, but for this it does not matter too much. Then a "blueprint file" formatted as discussed in the documentation. Briefly, in a way that residues left alone are in the form 1 M . (where residue 1 is a methionine), while the residues before and after the insertion are like 100 D L PIKAA D (meaning D100 is changed to a Asp in loop SS; not actual numbers or residue), while between entries are added 0 X L PIKAA G etc. To model loops, -generic_aa G is a good option as this tells the default rough fitting (centroid mode) to use glycine and not valine as generic residue.
The top pose will actually not have a superb score, so a 3 (regular) or more cycles (15 thorough) with Rosetta Relax application either of the whole protein or with a movemap specific to the linker, may be done.
There actually is an option in Rosetta remodel to have a domain insertion, i.e. it would swing around like a medieval flail. However, it will not repack and therefore properly dock so this option is a bad idea for this specific application.

thanks for your suggestion. I will learn that. Besides protein-protein docking, I want to physically connect two protein together. For example, I have known protein A structure, and protein B structure. Then I would like to connect the C terminal of protein A and N terminal of protein B. Then it will generate one final protein. If I get that connected protein, the next coming question is how to optimize the big structure to make it low-energy conformation because the connected structure may include clash residues or bad conformations?
Thank you.

I have expanded the answer, but it is very advisable finding out where the second protein goes first before adding the linkers. If the docking is flexible interacting sidechains will be repacked to be energy minimised anyway. The length of a straight amino acid (say βsheet) is about 3.5Å, so if your linker were 4 residues long, you'd have a maximum reach of 14Å, which is a lot, hence why local docking first is a better option.

Thanks so much for the explanation, which is very helpful. I will have a look at Rosetta remodel related functions.

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


text:
how could i find uniprot ac in blast?
I have made a search in p-blast.
I would like to search for the protein results in uniprot. how can I do it directly?
I can go to GENBANK but i would like to go to UNIPROT.
Thank you!!!

Ah ... yes thats tricky. The only thing I can suggest is using Biopython for Genbank and Uniprot because Biopython facilitates both.

thank you a lot for your answer!!

UniProt also provides a BLAST service (https://www.uniprot.org/blast), or you could use the Expasy BLAST too, at https://web.expasy.org/blast.
If you have performed your similarity search at NCBI or elsewhere, you can map the identifiers from that result to UniProt by using the UniProt IDmapping service at https://www.uniprot.org/uploadlists. This can also be used programmatically, cf https://www.uniprot.org/help/api_idmapping

This is all very useful information... just want to add that NCBI Protein BLAST also has "UniProtKB/Swiss-Prot" as a Database option.


text:
Matching files by chromosome during wdl scatter gather
I'm trying to implement a polygenic score pipeline in WDL and I'm quite new to this pipeline management.
The first step of the pipeline takes a file of GWAS summary statistics and splits it by chromosome.
The corresponding wdl task:
task split {
    input {
        File gwas
        String output_prefix
    }

    command {
        ./splitGwas -i ${gwas} -o ${output_prefix}
    }

    output {
        Array [File] gwas_by_chr = glob("${output_prefix}_*.assoc")
    }
}

The next part of the process is to compute posterior effects of SNPs by chromosome. The inputs to this task are going to be the files split by chromosome in the previous step and an LD matrix for the chromosome.
The corresponding wdl call would be something like this:
Array [File] ld_matrices = read_lines(file_of_ld_matrices_by_chr)

scatter(pair in zip(split.gwas_by_chr, ld_matrices)) {
    call sbayes.run {
        input:
            GCTB = gctb_executable_path,
            gwas = pair.left,
            ld_matrix = pair.right,
            output_prefix = out
    }
}

Since I'm grabbing the split files using glob in the previous step, I do not particularly know which element in the array corresponds to which chromosome. So I would be making a mistake if I zip the summary stats for one chromosome with the LD matrix corresponding to a different chromosome.
The programmatic way would be to iterate the scatter over chromosome name but then I would lose the ability to specify the dependency between the two tasks.
Is there a better way to do this?
Cross posted on biostars (https://www.biostars.org/p/434582/) so pardon me if that is discouraged.

Two questions: 1. Would it be possible for file_of_ld_matrices_by_chr to be in JSON format, instead of just a list a list of file paths? 2. What does the name of the file produced by splitGwas look like (ie. what's in between the prefix and the .assoc?

Yes, I can put the LD matrices in JSON
The file names produced by split GWAS have the names ${prefix}_chromosome.assoc so if the user supplied prefix string is T2D as in type 2 diabetes, the outputs are T2D_1.assoc, T2D_2.assoc ... T2D_22.assoc

You could use the sub function to obtain the chromosome from the globbed file names and then use that to retrieve the correct LD matrices from a Map (see read_map or read_json).
eg.:
Map[String, String] ld_matrices = read_json(file_of_ld_matrices_by_chr)

scatter(gwas_for_chr in split.gwas_by_chr) {
    String chr = sub(sub(gwas_for_chr, prefix, ""), "\.assoc$", "")

    call sbayes.run {
        input:
            GCTB = gctb_executable_path,
            gwas = gwas_for_chr,
            ld_matrix = ld_matrices[chr],
            output_prefix = out
    }
}

The file_of_ld_matrices_by_chr would have to look like this:
{
  "chr1": "/path/to/file1",
  "chr2": "/path/to/file2"
}

This works perfectly.


text:
Can you give me a detailed explanation about PrimeBlast?
I am trying to write a research paper on PRIMER-BLAST. I have read something about it on
https://www.ncbi.nlm.nih.gov/pubmed/22708584/
but still, it is not enough. So can you please tell me more about it or suggest some resources.

I'm not sure if it's just the language barrier and you actually meant something else, but...if you are unaware of basic tools for finding scientific papers (such as Google Scholar, Web of Science), and of good practices (such as following up citations of interest), then perhaps you shouldn't be writing a research paper just yet?

I would start at the Primer3 website here.
Once you understand all the variables on the website via the manual here and the associated publications Koressaar et al (2007) and Untergasser et al (2012) you will understand the gist of the algorithm used in primer design. This gives you optimised primers for a single strand of DNA/RNA
Primer-Blast looks in addition at maximising sequence conservation within an alignment of different sequences. The alignment is created using a Blast (I think its Mega-Blast. Thus Primer-Blast looks at one more parameter accounting for mutation heterogeneity which might be encountered in the PCR/RT-PCR and the number of primers it recommends is therefore reduced over say PRimer3.
Prove it for yourself take a given sequence and run it through Primer3 and then Primer-Blast and assess the differences and similarities.
References
Untergasser A, Cutcutache I, Koressaar T, Ye J, Faircloth BC, Remm M, Rozen SG (2012) Primer3 - new capabilities and interfaces. Nucleic Acids Research 40(15):e115
Koressaar T, Remm M (2007) Enhancements and modifications of primer design program Primer3 Bioinformatics 23(10):1289-91


text:
After completing the evaluation, how to predict PTM lysine or any other modification with unlabeled data for making a web server?
I have created a model with SVM classifier and evaluated it with 5-fold cross-validation. The dataset was sequence data containing lysine modified sites. Now I want to build a web server where anyone can input any sequence in FASTA format. Apparently, I have no label information at this point. So how can I will be able to extract features or use any feature selection method which requires label information? How can I extract or create the test set to predict? Suppose I have used f-score feature selection method in my model and selected 50 features from the given dataset. But for web server, having no label information, how can I select 50 features from the sequence that are provided by an user and then test it with my previously trained model?
As being a novice in the world of bioinformatics, I am seeking the answer which might be improper for others. But it will be very much helpful for me.
I am giving some link of some papers and its corresponding web servers.

iPTM-mLys: identifying multiple lysine PTM sites and their different type (https://academic.oup.com/bioinformatics/article/32/20/3116/2196525)
Web server link : http://www.jci-bioinfo.cn/iPTM-mLys
iHyd-PseCp: Identify hydroxyproline and hydroxylysine in proteins by incorporating sequence-coupled effects into general PseAAC
(https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5190098/)
Web server link : http://www.jci-bioinfo.cn/iHyd-PseCp

There are too many questions really and not enough background, but I presume the gist of the question is this, you've a trained model and you want to know how to turn it into a website. So really you are looking for basic information on how to reproduce your model via a website. Is that correct?

Yes, it is @Michael. Basically I am using sequence-coupling effect to extract features and ANOVA F-test to select features where label information is necessary. Therefore, to reproduce my model in a website, I have to extract the query protein sequence given by an user and apply F-test again on it and later use it for detecting the class. Is this the process or there is other rules for it?
I have given two links that contain some background I hope.

Interesting application. I would simply suggest using Flask, which is described in O'Reilly's Programming PyTorch, chapter 8. PyTorch in Production. Obviously this is Deep Learning rather than ML but the basic gist will hold. You should also think about building a docker for distribution (same chapter). We're in lockdown so I don't have this book with me.

Thanks @Michael for suggesting me a path! Hope this will work. I look forward to reading it.

Interesting application. I would suggest using Python's Flask system, which is described in a few pages in O'Reilly's Programming PyTorch, chapter 8. PyTorch in Production. Django is not really recommended (don't say that on Stackoverflow however). O'Reilly do an entire book on Flask, but thats to make the site pretty.
Obviously, the explanation here is Deep Learning rather than ML but the basic gist will hold. You should also think about building a docker for distribution which is easy (same chapter). We're in lockdown so I don't have this book with me.
The difficult part will be web architecture because you shouldn't have the CPU/GPU on the same server as the quiry interface, otherwise you are prone to a denial of service attack or just sytem overload. To do this you need to consider  Representational state transfer (RESTful) and that is not easy at all (at least for me). Flask facilitates RESTful, but thats like systems admin stuff. I guess before getting into RESTful I'd see how popular your website is, but if its the same server as your institute website you'd need to be really careful because an external user could bring your institutional website down (you wouldn't be popular).

That would be very helpful. As I have to go through the book you have mentioned, I am not quite sure it is related to feature extraction or feature selection procedures. I just wanted to know that without any class labels, how can I classify a protein sequence. Such as an F-test requires class labels for statistical analysis. A query protein which is completely new and not included in my trained model's dataset, is to be classified as either positive or negative.


text:
Can't install snakemake v5.15.0 from conda
I'm trying to install the latest version of snakemake from conda into a conda env, following the instructions here. I.e.
conda create -n snakemake -c conda-forge -c bioconda snakemake
However, this seems to only give me the option to install snakemake v5.3.0, while the latest no-arch version on conda is v5.15.0.
Any help greatly appreciated.
Solving environment: done

## Package Plan ##

  environment location: /home/ubuntu/anaconda3/envs/snakemake

  added / updated specs:
    - snakemake

The following packages will be downloaded:

    package                    |            build
    ---------------------------|-----------------
    aioeasywebdav-2.4.0        |        py36_1000          18 KB  conda-forge
    aiohttp-3.6.2              |   py36h7b6447c_0         544 KB
    appdirs-1.4.3              |   py36h28b3542_0          15 KB
    async-timeout-3.0.1        |           py36_0          12 KB
    attrs-19.3.0               |             py_0          39 KB
    bcrypt-3.1.7               |   py36h7b6447c_0          40 KB
    blas-2.14                  |         openblas          10 KB  conda-forge
    boto3-1.9.66               |           py36_0         105 KB
    botocore-1.12.189          |             py_0         3.4 MB
    datrie-0.8                 |   py36h7b6447c_0         147 KB
    docutils-0.16              |           py36_0         668 KB
    dropbox-9.4.0              |           py36_0         756 KB
    filechunkio-1.6            |           py36_0           7 KB  bioconda
    ftputil-3.2                |           py36_0          85 KB  bioconda
    google-api-core-1.16.0     |           py36_1          86 KB
    googleapis-common-protos-1.51.0|           py36_2          72 KB
    idna_ssl-1.1.0             |           py36_0           7 KB
    importlib_metadata-1.5.0   |           py36_0          48 KB
    jsonschema-3.2.0           |           py36_0          95 KB
    liblapacke-3.8.0           |      14_openblas          10 KB  conda-forge
    multidict-4.7.3            |   py36h7b6447c_0          68 KB
    networkx-2.4               |             py_0         1.2 MB
    numpy-1.17.0               |   py36h99e49ec_0          24 KB  r
    numpy-base-1.17.0          |   py36h2f8d375_0         5.2 MB  r
    pandas-1.0.3               |   py36h0573a6f_0         8.6 MB
    protobuf-3.11.4            |   py36he6710b0_0         635 KB
    psutil-5.7.0               |   py36h7b6447c_0         315 KB
    pygraphviz-1.3.1           |           py36_0         205 KB  bioconda
    pynacl-1.3.0               |   py36h7b6447c_0         1.1 MB
    pyrsistent-0.16.0          |   py36h7b6447c_0          93 KB
    pysftp-0.2.9               |           py36_0          31 KB  bioconda
    ratelimiter-1.2.0          |        py36_1000          12 KB  conda-forge
    rsa-3.1.4                  |           py36_0          87 KB  bioconda
    s3transfer-0.1.13          |           py36_0          79 KB
    snakemake-5.3.0            |           py36_1           4 KB  bioconda
    snakemake-minimal-5.3.0    |           py36_1         280 KB  bioconda
    sqlite-3.31.1              |       h62c20be_1         2.0 MB
    typing_extensions-3.7.4.1  |           py36_0          40 KB
    wrapt-1.12.1               |   py36h7b6447c_1          49 KB
    yarl-1.4.2                 |   py36h7b6447c_0         132 KB
    zipp-2.2.0                 |             py_0          12 KB
    ------------------------------------------------------------
                                           Total:        26.2 MB

The following NEW packages will be INSTALLED:

  _libgcc_mutex      pkgs/main/linux-64::_libgcc_mutex-0.1-main
  aioeasywebdav      conda-forge/linux-64::aioeasywebdav-2.4.0-py36_1000
  aiohttp            pkgs/main/linux-64::aiohttp-3.6.2-py36h7b6447c_0
  appdirs            pkgs/main/linux-64::appdirs-1.4.3-py36h28b3542_0
  asn1crypto         pkgs/main/linux-64::asn1crypto-1.3.0-py36_0
  async-timeout      pkgs/main/linux-64::async-timeout-3.0.1-py36_0
  attrs              pkgs/main/noarch::attrs-19.3.0-py_0
  bcrypt             pkgs/main/linux-64::bcrypt-3.1.7-py36h7b6447c_0
  blas               conda-forge/linux-64::blas-2.14-openblas
  boto3              pkgs/main/linux-64::boto3-1.9.66-py36_0
  botocore           pkgs/main/noarch::botocore-1.12.189-py_0
  ca-certificates    pkgs/main/linux-64::ca-certificates-2020.1.1-0
  cachetools         pkgs/main/noarch::cachetools-3.1.1-py_0
  cairo              pkgs/main/linux-64::cairo-1.14.12-h8948797_3
  certifi            pkgs/main/linux-64::certifi-2020.4.5.1-py36_0
  cffi               pkgs/main/linux-64::cffi-1.14.0-py36h2e261b9_0
  chardet            pkgs/main/linux-64::chardet-3.0.4-py36_1003
  configargparse     pkgs/main/noarch::configargparse-1.1-py_0
  cryptography       pkgs/main/linux-64::cryptography-2.8-py36h1ba5d50_0
  datrie             pkgs/main/linux-64::datrie-0.8-py36h7b6447c_0
  decorator          pkgs/main/noarch::decorator-4.4.2-py_0
  docutils           pkgs/main/linux-64::docutils-0.16-py36_0
  dropbox            pkgs/main/linux-64::dropbox-9.4.0-py36_0
  expat              pkgs/main/linux-64::expat-2.2.6-he6710b0_0
  filechunkio        bioconda/linux-64::filechunkio-1.6-py36_0
  fontconfig         pkgs/main/linux-64::fontconfig-2.13.0-h9420a91_0
  freetype           pkgs/main/linux-64::freetype-2.9.1-h8a8886c_1
  fribidi            pkgs/main/linux-64::fribidi-1.0.5-h7b6447c_0
  ftputil            bioconda/linux-64::ftputil-3.2-py36_0
  gitdb              pkgs/main/noarch::gitdb-4.0.2-py_0
  gitpython          pkgs/main/noarch::gitpython-3.1.1-py_1
  glib               pkgs/main/linux-64::glib-2.63.1-h5a9c865_0
  google-api-core    pkgs/main/linux-64::google-api-core-1.16.0-py36_1
  google-auth        pkgs/main/noarch::google-auth-1.13.1-py_0
  google-cloud-core  pkgs/main/noarch::google-cloud-core-1.3.0-py_0
  google-cloud-stor~ pkgs/main/noarch::google-cloud-storage-1.27.0-py_0
  google-resumable-~ pkgs/main/noarch::google-resumable-media-0.5.0-py_1
  googleapis-common~ pkgs/main/linux-64::googleapis-common-protos-1.51.0-py36_2
  graphite2          pkgs/main/linux-64::graphite2-1.3.13-h23475e2_0
  graphviz           pkgs/main/linux-64::graphviz-2.40.1-h21bd128_2
  harfbuzz           pkgs/main/linux-64::harfbuzz-1.8.8-hffaf4a1_0
  icu                pkgs/main/linux-64::icu-58.2-h9c2bf20_1
  idna               pkgs/main/noarch::idna-2.9-py_1
  idna_ssl           pkgs/main/linux-64::idna_ssl-1.1.0-py36_0
  importlib_metadata pkgs/main/linux-64::importlib_metadata-1.5.0-py36_0
  jinja2             pkgs/main/noarch::jinja2-2.11.1-py_0
  jmespath           pkgs/main/noarch::jmespath-0.9.4-py_0
  jpeg               pkgs/main/linux-64::jpeg-9b-h024ee3a_2
  jsonschema         pkgs/main/linux-64::jsonschema-3.2.0-py36_0
  ld_impl_linux-64   pkgs/main/linux-64::ld_impl_linux-64-2.33.1-h53a641e_7
  libblas            conda-forge/linux-64::libblas-3.8.0-14_openblas
  libcblas           conda-forge/linux-64::libcblas-3.8.0-14_openblas
  libedit            pkgs/main/linux-64::libedit-3.1.20181209-hc058e9b_0
  libffi             pkgs/main/linux-64::libffi-3.2.1-hd88cf55_4
  libgcc-ng          pkgs/main/linux-64::libgcc-ng-9.1.0-hdf63c60_0
  libgfortran-ng     pkgs/main/linux-64::libgfortran-ng-7.3.0-hdf63c60_0
  liblapack          conda-forge/linux-64::liblapack-3.8.0-14_openblas
  liblapacke         conda-forge/linux-64::liblapacke-3.8.0-14_openblas
  libopenblas        conda-forge/linux-64::libopenblas-0.3.7-h5ec1e0e_6
  libpng             pkgs/main/linux-64::libpng-1.6.37-hbc83047_0
  libprotobuf        pkgs/main/linux-64::libprotobuf-3.11.4-hd408876_0
  libsodium          pkgs/main/linux-64::libsodium-1.0.16-h1bed415_0
  libstdcxx-ng       pkgs/main/linux-64::libstdcxx-ng-9.1.0-hdf63c60_0
  libtiff            pkgs/main/linux-64::libtiff-4.1.0-h2733197_0
  libuuid            pkgs/main/linux-64::libuuid-1.0.3-h1bed415_2
  libxcb             pkgs/main/linux-64::libxcb-1.13-h1bed415_1
  libxml2            pkgs/main/linux-64::libxml2-2.9.9-hea5a465_1
  markupsafe         pkgs/main/linux-64::markupsafe-1.1.1-py36h7b6447c_0
  multidict          pkgs/main/linux-64::multidict-4.7.3-py36h7b6447c_0
  ncurses            pkgs/main/linux-64::ncurses-6.2-he6710b0_0
  networkx           pkgs/main/noarch::networkx-2.4-py_0
  numpy              r/linux-64::numpy-1.17.0-py36h99e49ec_0
  numpy-base         r/linux-64::numpy-base-1.17.0-py36h2f8d375_0
  openssl            pkgs/main/linux-64::openssl-1.1.1g-h7b6447c_0
  pandas             pkgs/main/linux-64::pandas-1.0.3-py36h0573a6f_0
  pango              pkgs/main/linux-64::pango-1.42.4-h049681c_0
  paramiko           pkgs/main/noarch::paramiko-2.7.1-py_0
  pcre               pkgs/main/linux-64::pcre-8.43-he6710b0_0
  pip                pkgs/main/linux-64::pip-20.0.2-py36_1
  pixman             pkgs/main/linux-64::pixman-0.38.0-h7b6447c_0
  prettytable        conda-forge/noarch::prettytable-0.7.2-py_3
  protobuf           pkgs/main/linux-64::protobuf-3.11.4-py36he6710b0_0
  psutil             pkgs/main/linux-64::psutil-5.7.0-py36h7b6447c_0
  pyasn1             pkgs/main/noarch::pyasn1-0.4.8-py_0
  pyasn1-modules     pkgs/main/noarch::pyasn1-modules-0.2.7-py_0
  pycparser          pkgs/main/noarch::pycparser-2.20-py_0
  pygraphviz         bioconda/linux-64::pygraphviz-1.3.1-py36_0
  pynacl             pkgs/main/linux-64::pynacl-1.3.0-py36h7b6447c_0
  pyopenssl          pkgs/main/linux-64::pyopenssl-19.1.0-py36_0
  pyrsistent         pkgs/main/linux-64::pyrsistent-0.16.0-py36h7b6447c_0
  pysftp             bioconda/linux-64::pysftp-0.2.9-py36_0
  pysocks            pkgs/main/linux-64::pysocks-1.7.1-py36_0
  python             pkgs/main/linux-64::python-3.6.10-hcf32534_1
  python-dateutil    pkgs/main/noarch::python-dateutil-2.8.1-py_0
  python-irodsclient conda-forge/noarch::python-irodsclient-0.8.2-py_0
  pytz               pkgs/main/noarch::pytz-2019.3-py_0
  pyyaml             pkgs/main/linux-64::pyyaml-5.3.1-py36h7b6447c_0
  ratelimiter        conda-forge/linux-64::ratelimiter-1.2.0-py36_1000
  readline           pkgs/main/linux-64::readline-8.0-h7b6447c_0
  requests           pkgs/main/linux-64::requests-2.23.0-py36_0
  rsa                bioconda/linux-64::rsa-3.1.4-py36_0
  s3transfer         pkgs/main/linux-64::s3transfer-0.1.13-py36_0
  setuptools         pkgs/main/linux-64::setuptools-46.1.3-py36_0
  six                pkgs/main/linux-64::six-1.14.0-py36_0
  smmap              pkgs/main/noarch::smmap-3.0.2-py_0
  snakemake          bioconda/linux-64::snakemake-5.3.0-py36_1
  snakemake-minimal  bioconda/linux-64::snakemake-minimal-5.3.0-py36_1
  sqlite             pkgs/main/linux-64::sqlite-3.31.1-h62c20be_1
  tk                 pkgs/main/linux-64::tk-8.6.8-hbc83047_0
  typing_extensions  pkgs/main/linux-64::typing_extensions-3.7.4.1-py36_0
  urllib3            pkgs/main/linux-64::urllib3-1.25.8-py36_0
  wheel              pkgs/main/linux-64::wheel-0.34.2-py36_0
  wrapt              pkgs/main/linux-64::wrapt-1.12.1-py36h7b6447c_1
  xmlrunner          conda-forge/noarch::xmlrunner-1.7.7-py_0
  xz                 pkgs/main/linux-64::xz-5.2.5-h7b6447c_0
  yaml               pkgs/main/linux-64::yaml-0.1.7-had09818_2
  yarl               pkgs/main/linux-64::yarl-1.4.2-py36h7b6447c_0
  zipp               pkgs/main/noarch::zipp-2.2.0-py_0
  zlib               pkgs/main/linux-64::zlib-1.2.11-h7b6447c_3
  zstd               pkgs/main/linux-64::zstd-1.3.7-h0b5b093_0

Proceed ([y]/n)?

Does it help if you include the specific version? conda create -n snakemake -c conda-forge -c bioconda snakemake=5.15

@mrhd - yes, it does, thanks for the suggestion

This is an example of the conda solver behaving incorrectly. In this case, it's choosing to pick a version that minimizes the number of dependencies required. Version 5.3.0 has 25 fewer dependencies, so if you let conda choose the version that's what it will pick. There's nothing you can really do about this other than specifying the version of snakemake that you want.
Update: It's very likely that mamba will correctly do what you want by default. Though I haven't used it, it's supposed to be a drop in replacement for conda here.

That works! Thanks Devon (I'm sure that's one of the things I tried, but perhaps when my channels were still in the wrong order...)

From next release on, installing via mamba is the default recommendation for Snakemake.

Actually, the issue does not seem to be the number of dependencies, but rather a lack of information in the stripped down current_repodata.json. If that is disabled and conda uses the full repodata, results are correct again, see https://github.com/conda/conda/issues/9905

I don't know if it means anything but...

It seems you are using anaconda instead of miniconda as required by the snakemake docs. Maybe reinstall conda using miniconda3.
conda is installing python 3.6 instead of the more recent 3.8... odd

anaconda is just miniconda with a bunch of additional packages installed.


text:
How can I interpret multiple alignment results?
I would like to ask how I can interpret objectively multiple alignment results. I have used JALVIEW  and TCOFFEE.

I would like to ask if consistency score plays any role in the interpretation?
I would like to ask for a region of a protein to be considered conserved, all the amino acids should have similar properties? So even if only one protein has a different amino acid (among the proteins that are included in the multiple alignment), then there is no conservation? for example, in this picture,
I am curious why there is no dot beneath some columns with almost identical amino acids (expect one sequence for example).

Good question but we will need to know what the purpose of conservation analysis is. Make primers? State "region of conservation"? Design a vaccine?

the ultimate goal is to find some information about the structure and the function of an unknown protein sequence.

Thanks I'll respond in a few days, my apologies for the delay in advance.

no worries!thank you so much for your help!i am grateful!


text:
"megahit: command not found" after correct installation with miniconda
I want to use MEGAHIT (git: https://github.com/voutcn/megahit) to assemble metagenome shotgun sequencing reads (illumina).
I install the package with the command conda install -c bioconda megahit as described in on the github page which yields:
Solving environment: done

## Package Plan ##

  environment location: /home/datura-workstation/miniconda3

  added / updated specs:
    - megahit

The following NEW packages will be INSTALLED:

  megahit            bioconda/linux-64::megahit-1.2.9-h8b12597_0

Proceed ([y]/n)?

Where I click yes and the packages is installed. However, after installation I am not able to use the command megahit as this results in the error: megahit: command not found.
I use Ubuntu 18.04.3 LTS and have installed miniconda3.
I reproduced the installation of megahit on a windows computer with a virtual machine using the same Ubuntu version. On my virtual Ubuntu I CAN use the megahit command.
Could anyone please help me solve this problem? I searched the internet but cannot find a good solution on the issue (also not on the github issue page of this package).

Why not just install Anaconda3?

Did you install conda without any custom modifications? Check the bin folder of the miniconda root folder for presence of the executable. If it is there then add the bin to PATH.

/home/datura-workstation/miniconda3/bin/megahit ? it should be there

Did you run make install after building?

@JoseFernandoLopezFernandez Anaconda is a package manager that replaces compilation. Therefore this is not necessary. This must be a PATH or environment problem

@ATpoint Ah, gotcha. Thanks

In minconda3/bin there is no executable file named megahit. (I checked the folder on my test virtual machine and there the file is indeed there.) That seems to be the problem. I can find the megahit on my workstation in the directory /home/datura-workstation/miniconda3/pkgs/megahit-1.2.9-h8b12597_0/bin. Why is it there? Can I just copy these files to the miniconda/bin?  All files of megahit are placed in this directory while they should be located in the miniconda3/bin directory.

EDIT: I re-installed megahit and now there the megahit files are present in the miniconda3/bin. However, I still receive the error.

EDIT2: after activating conda activate base, megahit became functional. Thanks all for the help!

type conda info --envs.  Look at your list of environments.  Say it's base that has megahit in it, do conda activate base, then just do megahit.

even in the base environment, the command megahit is not found.

I re-installed megahit and now when I activate conda activate base, megahit does work. Thank you!!


text:
Model the effect of treatment on mutation burden
I am struggling a bit to model the following problem. Basically, I would like to model tumor mutation burden (mutations per megabase, a continuous) as a function of treatment (categorical). The following is an example of how data looks like
patient-id therapy-1 therapy-2 therapy-3 cancer-type TMB
PID-1      1         NA        1         CT-1        1.3
PID-2      0         0         1         CT-1        0.14
PID-3      NA        0         1         CT-2        5
PID-4      1.        1         NA        CT-3        2.3
PID-5      1         1         NA        CT-3        10
...
...

1 - treated; 0 - not treated; NA - No information
I was trying multiple linear regression, but in most cases one patient got multiple therapies. This is where I am confused on how to interpret the results of linear regression. So far, I haven't included the cancer-type information into the model. However, IMO, including cancer-type into model will yield specific treatment effects in specific cancer type? Or is it over-expectation of what we can achieve from linear model?
Also, is there a better way to address the problem to get reasonable associations between treatment and TMB in this case?
Thank you for the help. Stay safe!

Could you post your model is it something like TMB ~ cancer-type*patient-id ?

Depends how you consider "linear regression", if you mean Pearsons, Spearman's analysis then general linear modelling (GLM) is a far better approach if you have continuous data. I don't really understand TMB, but that is a critical output then GLM will work. Logistic regression is an alternative approach data type depending. You have to perform a Q-Q plot and the residual against against any transformations you make (you might need to transform TMB). The lower the residual the better, the straight the Q-Q plot the better.
To explain how GLM works long hair and height correlate in a straight regression analysis. The implication is if you want to be tall - cut your hair! Of once you account for gender and send it through a GLM then gender and hairlength will correlate far stronger than hair length and height.
Mac OSX has a very easy to use pointy clicky, but highly professional, GLM package on its AppStore called Wizzard (bit expensive). It doesn't do general linear mixed modelling however. Everything does GLM of course, its just this package is really easy to use and you can figure out GLM by trial and error.


text:
How to make Venn diagram for the Macs peak calling output of two data sets?
I have two outputs of macs2 peak calling for two of my data sets. I wanted to plot the Venn diagram to see how many peaks are shared. I mainly work on Unix/Linux. Do you know any way that I can have those diagrams plotted?
Much Appreciated

Check out the tool Intervene. It's exactly designed for this. Github link.

Dear @Cat if you describe Intervene for those of us who are unfamiliar with it?

Venn diagrams can be helpful, but for peak calling they can be somewhat misleading.
This topic is addressed in the DiffBind documentation, in the "Comparison of occupancy and affinity based analyses" section.
Suffice to say, asking "which peaks are shared?" is not the same as asking "is there evidence for differential peaks, and where are they?".
For this reason, performing a differential binding analysis might be a better option.
If you're just looking for Venn diagram plotting software, though, the VennDiagram (note, I have personally contributed to this tool in the past) and UpSet R packages are useful tools.


text:
Gene information from chromosome location
I have a file that contains species name, chromosome, start and end locations like the following format
hg19,chr8,77778733,77779026. I would like to know some information about the regions 77778733-77779026. How can I quickly find out if it lies in an intron or exon of a gene? Or it is in the upstream or downstream regions of a gene? Is there a tool or web service from where I can get this information?

you could probably use bedtools intersect or bedops bedmap

How can I quickly find out if it lies in an intron or exon of a gene?

I have an answer on Biostars that might help: https://www.biostars.org/p/124515/#124522

Or it is in the upstream or downstream regions of a gene?

One option is to build a table of annotated bins and then map your regions-of-interest to those annotations.
These annotations would consist of whether the bin is inside a gene, or, after picking the nearest gene by strand, locate the bin to an upstream region of the gene TSS or to an downstream region of the gene TES.
To start, one could download stranded gene annotations:
$ GENCODE_URL=ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.basic.annotation.gff3.gz
$ GENCODE_GENES=${PWD}/../data/gencode.v28.genes.bed
$ wget -qO- "${GENCODE_URL}" | gunzip -c | awk '($3=="gene")' | gff2bed - | awk -v FS="\t" -v OFS="\t" '{ n=split($10,a,";"); p=a[4]; n=split(p,a,"="); $4=a[2]; print $0; }' | cut -f1-6 > ${GENCODE_GENES}

Chop up nuclear chromosomes of the genome of interest into 1kb bins (or whatever gradation of bin size you prefer):
$ ASSEMBLY=hg38
$ GENOME=${PWD}/../data/${ASSEMBLY}.bed
$ CHOP_SIZE=1000
$ CHOPPED_GENOME=${PWD}/../data/chopped.${ASSEMBLY}.bed
$ fetchChromSizes ${ASSEMBLY} | awk -v FS="\t" -v OFS="\t" '{ print $1, "0", $2 }' | grep -ve "\.*_\.*" | grep -v "chrM" | sort-bed - > ${GENOME}
$ bedops --chop ${CHOP_SIZE} ${GENOME} > ${CHOPPED_GENOME}

Then assign the midpoints of bins inside a gene, or to the nearest of any particular set of intergenic categories.
Here's a script I wrote which uses the output of closest-features to do that. This could be easily modified for any set of distance categories:
#!/usr/bin/env python

import sys
import numpy as np

'''
Input should be bin midpoints run through closest-elements, against some set of BED6+ (stranded) gene annotations
'''

def assign_to_category(dist):
  if dist < -50000:
    return '-50kb_to_-Inf'
  elif dist < -20000 and dist >= -50000:
    return '-20kb_to_-50kb'
  elif dist < -10000 and dist >= -20000:
    return '-10kb_to_-20kb'
  elif dist < -5000 and dist >= -10000:
    return '-5kb_to_-10kb'
  elif dist < -1000 and dist >= -5000:
    return '-1kb_to_-5kb'
  elif dist < 0 and dist >= -1000:
    return 'TSS_to_-1kb'
  elif dist == 0:
    return 'Gene'
  elif dist > 0 and dist <= 1000:
    return 'TES_to_+1kb'
  elif dist > 1000 and dist <= 5000:
    return '+1kb_to_+5kb'
  elif dist > 5000 and dist <= 10000:
    return '+5kb_to_+10kb'
  elif dist > 10000 and dist <= 20000:
    return '+10kb_to_+20kb'
  elif dist > 20000 and dist <= 50000:
    return '+20kb_to_+50kb'
  elif dist > 50000:
    return '+50kb_to_+Inf'

for in_line in sys.stdin:
  (bin_str, upstream_gene_str, upstream_gene_dist, downstream_gene_str, downstream_gene_dist) = in_line.rstrip().split('|')
  try:
    upstream_gene_dist = int(upstream_gene_dist)
  except ValueError:
    upstream_gene_dist = np.nan
  try:
    downstream_gene_dist = int(downstream_gene_dist)
  except ValueError:
    downstream_gene_dist = np.nan
  upstream_gene_elems = upstream_gene_str.split('\t')
  downstream_gene_elems = downstream_gene_str.split('\t')

  '''
  We have six possible scenarios, depending on the minimum distance and strandedness of closest elements
  '''

  '''
  1) If either distance value is NaN, the bin is assigned to the non-NaN distance element
  '''
  if upstream_gene_dist is np.nan and downstream_gene_dist is np.nan:
    raise SystemError("Bad NaN input on:\n[{}]\n".format(in_line.rstrip()))
  elif upstream_gene_dist is np.nan:
    gene_id = downstream_gene_elems[3]
    if downstream_gene_elems[5] == '+':
      # bin is closer to downstream TSS
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(downstream_gene_dist))), gene_id))
    else:
      # bin is closer to downstream TES
      sys.stdout.write('{}\t{}\t\n'.format(bin_str, assign_to_category(abs(downstream_gene_dist)), gene_id))
    continue
  elif downstream_gene_dist is np.nan:
    gene_id = upstream_gene_elems[3]
    if upstream_gene_elems[5] == '+':
      # bin is closer to upstream TES
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(upstream_gene_dist)), gene_id))
    else:
      # bin is closer to upstream TSS
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(upstream_gene_dist))), gene_id))
    continue

  '''
  2) If either distance value is zero, the bin is assigned to "Gene" to state overlap with a gene
  '''
  if upstream_gene_dist == 0:
    gene_id = upstream_gene_elems[3]
    sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(0), gene_id))
    continue

  if downstream_gene_dist == 0:
    gene_id = downstream_gene_elems[3]
    sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(0), gene_id))
    continue

  '''
  3) If the upstream strand is positive, if the downstream strand is positive, then:

   -- the upstream distance is a signed (-) distance to the upstream gene's TES
   -- the downstream distance is a signed (+) distance to the downstream gene's TSS

  '''
  if upstream_gene_elems[5] == '+' and downstream_gene_elems[5] == '+':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TES
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(upstream_gene_dist)), gene_id))
    else:
      # bin is closer to downstream TSS
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(downstream_gene_dist))), gene_id))
    continue

  '''
  4) If the upstream strand is negative, if the downstream strand is positive, then:

   -- the upstream distance is a signed (-) distance to the upstream gene's TSS
   -- the downstream distance is a signed (+) distance to the downstream gene's TSS

  '''
  if upstream_gene_elems[5] == '-' and downstream_gene_elems[5] == '+':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TSS
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(upstream_gene_dist))), gene_id))
    else:
      # bin is closer to downstream TSS
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(downstream_gene_dist))), gene_id))
    continue

  '''
  5) If the upstream strand is positive, if the downstream strand is negative, then:

   -- the upstream distance is a signed (-) distance to the upstream gene's TES
   -- the downstream distance is a signed (+) distance to the downstream gene's TES

  '''
  if upstream_gene_elems[5] == '+' and downstream_gene_elems[5] == '-':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TES
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(upstream_gene_dist)), gene_id))
    else:
      # bin is closer to downstream TES
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(downstream_gene_dist)), gene_id))
    continue

  '''
  6) If the upstream strand is negative, if the downstream strand is negative, then:

   -- the upstream distance is a signed (-) distance to the upstream gene's TES
   -- the downstream distance is a signed (+) distance to the downstream gene's TES

  '''
  if upstream_gene_elems[5] == '-' and downstream_gene_elems[5] == '-':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TSS
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(upstream_gene_dist))), gene_id))
    else:
      # bin is closer to downstream TES
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(downstream_gene_dist)), gene_id))
    continue

  '''
  Sanity check: If we get here, something went wrong
  '''
  raise SystemError("Bad logic on:\n[{}]\n".format(in_line.rstrip()))

To use it, e.g.:
$ MP_LEFT=500
$ MP_RIGHT=-499
$ MP_LEFT_UNDO=-500
$ MP_RIGHT_UNDO=499
$ CATEGORIZED_BINS=${PWD}/../data/categorized_bins.${ASSEMBLY}.bed
$ ./assign_chops.py < <(bedops --range ${MP_LEFT}:${MP_RIGHT} --everything ${CHOPPED_GENOME} | closest-features --dist - ${GENCODE_GENES}) | sort-bed - | bedops --range ${MP_LEFT_UNDO}:${MP_RIGHT_UNDO} --everything - > ${CATEGORIZED_BINS}

The output from this file annotates bins with the closest feature to a gene. It looks something like this:
$ more categorized_bins.hg38.bed
...
chr1  32000 33000 TES_to_+1kb RP11-34P13.3
chr1  33000 34000 TES_to_+1kb FAM138A
chr1  34000 35000 Gene  FAM138A
chr1  35000 36000 Gene  FAM138A
chr1  36000 37000 Gene  FAM138A
chr1  37000 38000 TSS_to_-1kb FAM138A
chr1  38000 39000 -1kb_to_-5kb  FAM138A
...

For example, the region chr1:33000-34000 is closest to the 1kb window downstream of the TES of FAM138A.
Once you have this file, you can quickly map it to your sorted regions-of-interest file:
$ bedmap --echo --echo-map --delim '\t' regions-of-interest.bed categorized_bins.hg38.bed > answer.bed

The file answer.bed will contain your regions of interest and the annotated bins which overlap them.

How can I quickly find out if it lies in an intron or exon of a gene?
  Or it is in the upstream or downstream regions of a gene? Is there a
  tool or web service from where I can get this information?

Since you emphasise speed and suggest a web service, I'm guessing you just want to check a few individual loci. In that case, just use a genome browser like those from UCSC or Ensembl. If you look up those coordinates in UCSC you can see the sequence falls within the last exon of the ZFHX4 protein coding gene and other information.

There is already some very useful answers out there but in my experience Homer does a very efficient job at annotating bed files. This gives you a lot of details about nearest transcription start sites or exonic regions and more.
Be sure to convert your data to bed format first


text:
Pindel has an extremely large output
I am wondering if I am doing something wrong here.
I ran Pindel (alongside two other structural variant [SV] callers) on four WGS samples. The two other callers report approximately 50000 SVs and have a decent running time.
However, Pindel takes long to run and its output files are very large (VCF files of 9 GB, more than a million detected SVs). I never used Pindel before, but I am wondering if I am doing something wrong here.
My command to call Pindel (from Python to cmd) is:
pindel -f {} -i {} -o {} -T {}.format(reference, config_file, output_path, threads)

PS: I use Pindel version 0.2.5b9
I am mainly confused how Pindel reports so many more calls in comparison to the other two callers. I appreciate any suggestions or remarks!

Some other SV callers have higher thresholds for SV sizes (I think Delly does not call anything small). If you would filter Pindel variants by the same treshold, would you get something comparable?

I see. I could try to play around with the settings a little indeed. Thanks for the info!

You could just check histograms of SV sizes, that would be probably the quickest way to figure if the difference has something to do with SV sizes (I am not sure, it was just a guess)


text:
How to run Cox proportional Hazard model for each group in R?
My population is divided into 10 groups. I want to run the cox-proportional hazard model separately for each of these 10 groups so that in the end i can compare which groups has highest hazard ratio depending uopn the covariates.
How can i do this in R? Earlier i was using this formula.
Model_1 <- coxph(serve_object ~ GROUP + age+ Smoking+ Diabetes, data = DF
summary(Model_1)

But the output is very confusing. it gives hazard ratio of 9 out of 10 groups. I think it dealt group as covariate.

The fact that you supplied GROUP as a covariate, you are getting 9 coefficients as one of the ten levels is used as the reference level.
If you would like to fit one model per group (and in this scenario you would not be correcting for GROUP), you should first split your data frame and fit your models on the resulting list's elements:
lapply(split(DF, DF$GROUP),
       function(x) coxph(serve_object ~ age + Smoking + Diabetes, data = DF))

Example with the mtcars data set when grouped by number of cylinders, cyl having values of 4, 6 or 8:
> lapply(split(mtcars, mtcars$cyl),
+        function(x) lm(disp ~ hp, data = x))
$`4`

Call:
lm(formula = disp ~ hp, data = x)

Coefficients:
(Intercept)           hp
    59.0369       0.5579

$`6`

Call:
lm(formula = disp ~ hp, data = x)

Coefficients:
(Intercept)           hp
   290.9178      -0.8799

$`8`

Call:
lm(formula = disp ~ hp, data = x)

Coefficients:
(Intercept)           hp
   320.2083       0.1572

Can you help me in interpreting result.
.$Group3
Call:
coxph(formula = so_ukb ~ age_at_recruitment.210220.0 + Smoking +
Hypertention, data = UKB_FILTRD)
                            coef exp(coef) se(coef)      z      p

age_at_recruitment.210220.0  0.06961   1.07208  0.00141  49.35 <2e-16
Smoking.L                   -0.21220   0.80880  0.01335 -15.90 <2e-16
Hypertention.L              -0.44221   0.64262  0.01361 -32.49 <2e-16
Likelihood ratio test=5142  on 3 df, p=< 2.2e-16
n= 448287, number of events= 11377
(3270 observations deleted due to missingness)

I am just confused that how to get one cumulative hazard ratio by adjusting model for these covariates. It is giving HR for age, smoking, hypertention separatively.. I want to see just cumulative.

Thank you for above lapply code, it worked so well.

Let's take the coefficient you get for age: The exp() of this coefficient is the change in in HR when age changes by one unit and while all the other variables are accounted for. In that sense, it takes all the variables into account so can be said to be "cumulative", although I am not sure such terminology exists for linear models. You might want to use interaction terms if your covariates depend on each other.

Can you give an idea of how we can export the summary(coxph) into csv?

or may be first convert it into data frame and then export df intom csv.

You can give a try with the broom package. Its tidy() function would convert output of lm() and many other related functions to data frames. Would probably work with the coxph()'s output.


text:
remove sequences from fasta file matching a string in the header
I have a file with 16S sequences. some headers contain species information. For my purposes I would like to exclude a number of species from the file, therefore I would like to do a pattern matching within a 16S header (that contains other characters beyond the species name) to then remove that header and the following sequence.
So, I have:
>NCBI_ID1;high_rank_taxonomy;...;species_name1
ATCCTT
>NCBI_ID2;high_rank_taxonomy;...;species_name2
ATTTTCCC
>NCBI_ID3;high_rank_taxonomy;...;species_name3
TTTCCCCAAAA
>NCBI_ID4;high_rank_taxonomy;...;species_name4
CCCTTTAAA

I would like to exclude, say, species 1,2 and 4 to retain:
>NCBI_ID3;high_rank_taxonomy;...;species_name3
TTTCCCCAAAA

I found this solution in an adjacent question:
awk 'BEGIN{while((getline<"ids.txt")>0)l[">"$1]=1}/^>/{f=!l[$1]}f' seq.fa

where I have my species names in the file ids.txt and sequences in file seq.fa
but that would work if the whole of my header was just the species name (contents of each line in ids.txt), which is not the case.
I thought I could be smart and use *species_name1* instead of species_name1 as contents of ids.txt but that apparently doesn't work.
Could someone help me with modifying the current solution to accommodate the existence of additinal strings in the headers?
Many thanks!
A

my awk is really bad.. would doing it in R or python help?

Say yes to that, @SW's Python is really good. The other way in is a Perl one-liner.

Or, you know, Heng Li's bioawk would be a good middle ground.

@StupidWolf sure thing, I am very new to bioinformatics, so I am happy to operate anything that does the job and learn it in the process : ) R would be preferable at the moment

Hi @Irrussional, see whether the R solution posted below works? I think it should be ok

yes, on it currently! thanks

you may use simple R commands, using the 'seqinr' package
#load package
library(seqinr)
#load file containing sequences
data<-read.fasta("filename.fa")
#create a vector containing species names: these are the last part of the string
vec.names<-unlist(lapply(strsplit(names(data), ";"), function(x)x[length(x)]))
#find names to keep: indices which are not in the species to remove
species.to.remove<-c("species_name1", "species_name2", "species_name4")
vec.tokeep<-which(! vec.names %in%  species.to.remove)
#write the final output
write.fasta(sequences=data[vec.tokeep], names=names(data)[vec.tokeep], file.out="output.fa")


text:
Protein PDB/CIF text file reader that can fold a same section into one line
Is there a software that can recognise the protein PDB/CIF files, and fold multiple lines of a same section into one line, so that it is easier to handle? I need to understand the format in CIF so as to create my own CIF file.
For example, when Notepad++ is used to read a .py file, it can fold the same section into one line.
So for a PDB/CIF file, e.g. all the atom lines for a same chain can be folded into one line, and I can expand it if needed.
Like this picture, assuming "This PC" is the whole protein, while the other objects are the different chains.

You could try opening the PDBXML version of the file in an XML viewer that lets you collapse sections.

Thank you, the PDBXML file is 10 times the size of a CIF file. The 5TX1 has 100 MB. Can you recommend a software that is good for PDBXML?

Programs that work with PDB and mmCIF files may present the structure as a tree, for example:

(this is a Find dialog from ccp4MG)
but I'd be surprised if a text editor had such a feature. Normally, you don't edit PDB files in a text editor.

Thank you, I need to understand the format in CIF so as to create my own CIF file

@lanselibai do you want to create an mmCIF file in a text editor? why?

I am trying to do homology modelling for a virus capsid based on an existing CIF (5TX1, as the template) in the PDB database. The CIF template contains the assembly information, which I want to understand. But the CIF file is too long, which is very easy to get lost.

@lanselibai to understand assemblies you may start with Introduction to Biological Assemblies and the PDB Archive at RCSB's PDB101.

Sure, I will read it! Thank you!


text:
Can I use BioEdit to quantify the sequence similarity?
I am using BioEdit. It can perform an align of multiple sequences. Can this sequence similarity be outputted as a value, e.g. 70%?

No Bioedit is an alignment editor only. It is not in anyway an analysis program.
Dump the format in fasta or mega (its the option is available) and import it into Mega X here. You can quantify the diversity as a pairwise matrix, population analysis (if you have a clear a priori that you can subdivide your sequences into different groups) or phylogeny ... which provides diversity as mutations per alignment position. Its pointy and clicky and easy to use.


text:
Generate an appropriate output from DNaSP6 to Arlequin program
I am performing a haplotype and nucleotide diversity analysis of my Sanger sequences (from Watermelon mosaic virus coat protein region) using the Arlequin program. My aim is to see the probability that two randomly sampled alleles are different, and the expansion and genetic variation of the virus respectively.
Some of my sequences have gaps (-) because the PCR reaction could not cover all my region of interest. After making the alignment in Mega and export in a Nexus file to DnaSp6 I saw that there are many options to export to the Arlequin program after DnaSP6. I do not know what is the best option to generate the output from DNAsp6 to Arlequin because I must choose one option for the sites with the alignment gaps option and another option for the invariable sites option.
These are the options:
Sites with Alignment Gaps options:

Not considered: These sites are ignored (complete deletion).
Considered: Gaps are considered just like another nucleotide variant (fifth state).
Only gaps are considered: Only gaps information is considered to build haplotypes.

Invariable Sites options:

Removed: Invariable (monomorphic) sites will not be included in the output file.
Included: Invariable (monomorphic) sites will be included in the output file.

From my point of view I think that it is more certain for my analysis (above explained) to use the Not consider gaps option (because they can buffer the information of the other sequences) and include invariable sites option because I want to analyse these parameters of my dataset.
Is it sensible? I would like to know what is the best option for my dataset.

Yes, your question makes alot of sense and its a good question because hits at the heart of differences between population genetic analysis and phylogenetic analysis. However, when you speak about gaps you are referring to 'indels' (insertion/deletion).
Indels The conservative approach is to simply delete sites with gaps. Arlequin uses a lot of maximum likelihood solutions and gaps would automatically be deleted in an Arlequin max. likelihood analysis, because it is the output per alignment per matrix. If the matrix doesn't have a 'gap'.
Indels and pop gen Thats tricky, it is better to remove them. HOWEVER the indel may contain underlying population signal, if its a maximum likelihood analysis (pop gen does this too) I suspect they will be ignored. If its pairwise, such Fst including indels might give you a better answer because there is more polymorphisms to assess population structure. So if you're doing a pairwise analysis, I'd do it both ways and then report the methodology correctly if there are differences in results and simply explain the theory.
Invariable sites This is different from an indel. An invariable site is a site which carries no phylogenetic signal, it carries information but a phylogenetic analysis could draw any conclusions because all the nucleotides/amino acids are identical at a give alignment position.
The output you obtain in allele frequency will be different and here's the rub ...
Phylogenetics analysis Arlequin is really pop gen, but it will contain some phylogenetics programs. If you see anything with 'Jukes Cantor correction', 'Kimura 2 (or 3) parameter model', 'HGW model', '"REV" model' (eek I can't remember its name any its a reverible 6-parameter model of nucleotides),  you are dealing with phylogenetics. Keeping invariable sites is normal here because the rate matrix is calculation "per alignment site".
Population genetics analysis Fst, Fis, linkage disequilibrium statistics, REMOVE THEM. Why? Because you want SNPs (single nucleotide polymorphisms). The invariable sites will get in the way of the statistics. For example is the locus recombines and you want to assess a linkage disequilibrium statistic (recombination), invariable sites will recombine, but the algorithm can't see this because the nucleotides are all the same. The only thing the algorithm can see and analyse are SNPs, which 'representative' of all the surrounding invariant sites. Thus removing invariants is essential. Arelquin might do this automatically (but I don't know the preparatory algorithms), but if you submit multiple loci Arlquin will want to know which are different loci. It will then recode the loci automatically if a single locus comprises more than one SNP. We are going into the underlying theory of linkage disequilibrium here and I'm not sure if say Fst/Fis is particularly bothered with this or is just concerned with the populations of genetic diversity.
Summary So gaps ... if they are genuine indels and not just 'N', ie indeterminant nucleotide calls, remove them for phylogenetics. If its pop gen do the calculation both ways and assess the output, if there is no difference report the 'gap/indel removal' result.
Invariant sites .. for pop gen you remove them but Arlequin might do this automatically to prepare a data set for analysis (I don't know). For phylogenetics you never remove invariant sites.
Goodluck

Wow, thanks you Michael for your answer. It looks great but I have some doubts. Perhaps my lack of knowledge in this field does not allow to me unterstand well your answer. I am analyzing the Fst fixation index too. I though that conservative regions may introduce information about my sequences and if all my sequences consider this regions in the FST analysis the results will be sense, aren't them? Thank you for your answer.

Hi @AdrianPL you are welcome. For Fst I would remove invariants (no change, if there is one SNP+ keep it) because it give stronger population structure, but indels/gaps are more complicated (because they might contain information and honest I would just leave them in). Arlequin was particularly written in the days of microsatellites where every position was an allele. The Fst I know is the classic pairwise calculation, there are maximum likelihood versions, but I don't know how they work. Fixation, if you mean homozygosity, also consider Fis (which is what is designed for-separate question).

Thanks you Michael. I removed the indels of my dataset and I think that your explanation is very sense because the Fst values gave to me a stronger population structure than those produced by Rstudio or Dnasp6. I prefer make all analysis with the same program, Arlequin because It gives to me Standard error and permutation values. Is there any article that remove the invariant sites to calculate the Fst value? Thanks you!!

In NGS era there are lots, anything on NGS Plasmodium falciparum population genetics, but it would standard practice and Dnasp6 isn't regarded in the same calibre as Arlequin, which is an 'old school standard'. R however has a lot of pop gen approaches, so its difficult to generalise here. You could ask this as a separate question if you thought it was critical, but from my point of view Fst is a measure of population migration, not genetic diversity, so I would see the removal of invariants as a positive step from an emperical mathematical perspective.

Thanks you a lot Michael!


text:
Is there a single cell RNAseq equivalent of GTEx or TCGA?
Or do I need to find individual studies and obtain data the long way!

Does the question read as "Is there a database of Human scRNA seq data?"? I think https://bgee.org/ contains some scRNA seq data, probably other databases too. Did you try to simply google it?

I think the most well-known scRNA dataset that is similar to GTEx in the Human Cell Atlas (https://www.humancellatlas.org). It is pretty much still ongoing but there are some preliminary data.


text:
Map gene IDs to protein structure if exist
I want to map some gene ids to their proteins structures if they exist.
The genes are rows of this CSV. There are other information that may help also such as Product description but that may lead to anbiguities. I am looking for a programmatic way, biopython or pdb API?

From EMBL SIFTS you can download a table that goes from Uniprot ID to PDB and chain and numbers that can give you any offset. Do note that offset correction is a problem for about 10% —I use python's pymol2 module to fix these (yup, silly but best option). Also if the organism is common you might wanna check out SwissModel's premade threaded models of close homologues(20% id min)

Not sure if you want the structure exactly or a handle of some sort, but the best way i know to get ALL your 5000 identifiers to map to uniprot/pdb is uniprot's mapping service.
You can either copy your ids straight from excel or parse with python ex.:
import pandas as pd
import os

file_name = 'aap7847_Zhang_SM_Table-S5.xlsx'
sheetname = 'Table S5'
df = pd.read_excel(file_name, sheet_name=sheetname, skiprows=1)
ids = list(df['Gene_ID'])

name_of_file_to_save_to = "geneIds.txt"
with open(name_of_file_to_save_to, "w") as newf:
    for id in ids:
        print("{}\n".format(id), file=newf)
    newf.close()

You can then save the Uniprot's response and either add it as a column to your file or map to other parts of Elixir ecosystem with SIFTS, for example, as @Matteo Ferla says.

In the case that you maybe need a single id's information, you could hit their API too, although this wouldn't work for bulk or at least didn't work for me. You can map the uniprot accession you get from this to PDB or whatever else you need!
import requests
def get_uniprot_accession(geneid, format='list'):
    resp = requests.get(
        "https://www.uniprot.org/uniprot/?query={}&sort=score&format={}".format(geneid, format)).text
    if format == 'list':
        resp = resp.strip()
    return resp


text:
How to convert Illumina IDs to Gene Names?
I am new to the area of bioinformatics, so apologies if this is too obvious of a query.
I need to analyze the RNA Seq data from GSE98455.
The dataset is of the following format:
|---------------------|------------------|
|    Some Id          |      Counts      |
|---------------------|------------------|
|    13101.t00001     |       392        |
|---------------------|------------------|
|    13101.t00002     |        20        |
|---------------------|------------------|

I would like to extract the counts for certain genes only. So how do I match the ID column with gene names? Is there a reference table or something I need to look at?
Side Note: Usually for microarray data, I refer to the GPL files which provide a reference between accession ID and Gene Names. I was not able to do that for this RNA Seq Data set.
Update May 2, 2020: I reached out to the authors earlier today. I looked at the paper associated with this dataset, it seems the annotation belongs to Oryza sativa japonica reference genome v 6. I can't find anything related to it on the web. If I get any response from the authors, I will update it here.

you have to ask the authors. I checked https://plants.ensembl.org/Oryza_sativa/Info/Index?db=core and the gffs don't have any transcripts that match it

This is not illumina IDs. you can read the data processing, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM2596386, they aligned it to the ensembl genome and i have no idea what annotation is used


text:
vst() from DESeq2 vs voom() from limma
I have used both to transform my leukaemia RNA-Seq data for subsequent hierarchical clustering. The result is quite different. Some subtypes of leukaemia only form a cluster (or at least sit closer) in one of the transformed data set. While I understand both can be used to transform counts into gene expression matrix that can be use for clustering, I want to know more about the difference. From my really poor statistic background I assume voom() is very similar to rlog(CPM) and can generate a microarray-like data set. On the other hand, the name of vst() is quite self-explanatory and is described to be better at handling varying size factors than rlog() does.
What are the actual difference and is there one that you think is better at handling jobs like clustering analysis?
I am hoping this to be a more educational question than technical one. If it violates any rules please delete it.


text:
What does it mean if aligned genomes are not of the same length?
I am aligning 4 genomes of Psuedomonas chlororaphis using progressive mauve software. Command used:
progressiveMauve --output=outputfile --disable-backbone chl189.fas chl17415.fas chl50083.fas chlTAMOak81.fas

At the end, the resulting aligned genomes are not of the same length. What could be the potential problems?

This is a global alignment. Which is why im surprised that they aren't all the same length after alignment?
Also, while analysing the subsequent genomes, I noticed that some region which were present in only 2/4 genomes were misaligned in the end result. How do I proceed troubleshooting this problem?
I am a student and new to this field, could you please tell me how to proceed with manually checking to see whats wrong?

I’ve seen cases like this. For me, there were plasmids at end that didn’t get aligned.

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

Hi @SethCommichaux please do remember to provide supporting information for your answer, i.e. not simply to confirm the OP's observation but offer a solution.


text:
how to evaluate neural networks
I have made some architectures of neural networks in order to classify proteins in three categories. I have calculated true positive, true negative, false negative, false positive, Mathews, sensitivity, specificity, accuracy for each category. My question is what is the right way to evaluate the architecture. Except for the c-index, could I calculate a cumulative sensitivity, accuracy etc?or should I evaluate separately for each category? Thank you!

Wow, I would just go for accuarcy as the singular most important measure, but be careful of sampling bias.
If the classification training had an unbiased sample* and you have an accuracy of >90%, in general ANN that wouldn't be great, but in bioinformatcs thats really good. If your training comprises sampling bias that gets tricky and you either need to accomodate that in your accuarcy ROC etc ... or else retrain to eliminate the sample bias.
You can account for all of them of course by wrapping everything into a single analysis.
*, by this I mean that each of the 3 classifications comprises 33% in each training set. If this is badly skewed, this requires consideration and gets complicated. If it is precisely 1/3 each then 'accurary' is the key stat.

thank you!how can i calculate accuracy for all the three categories? is it correct to calculate the mean ? (sorry in advance for this question)

Thanks, firstly this looks very cool work, again... its cool.  Did you do 3 independent training sets, each with a single classification OR a single training with 3 classifications?

a single training with 3 classifications.thank you!

I get it, sorry. Yes and no. Yes: for your own internal purpose yes and a line or two in the report/paper. You don't want major skews between the classifications (these are different skews from the sampling skews). No: for the overall accuracy score, its one overall number (obviously one which isn't hiding a skew). Essentially the answer is its just one number, but you need to be confident its representative. I hope thats ok.

If you get 90% accuracy with a 33% sample for training of each classification, then its taken as read that the misclassification has reasonable distribution. Leading bioinformatics publications have present been as low as 80%.

so to sum up, instead of  making a 3*3 matrix to measure accuracy, I could calculate the mean (because now I have constructed matrices for each category). and if this is higher than 90% (but even less) then my model is good. am I correct? Thank you so much!!!

Basically yes in my opinion 90% would be fantastic and you are sorted, look at past papers (Bioinformatics journal, to verify the accuarcy against their work ...  Your confusion matrix (the formal name for a 3*3 matrix) helps diagnose your training, if your accuracy is low, scikit-learn has an automatic function for this. Personally I think 85%+ is good, but just my opinion.

thank you!in order to select one architecture than another how much should be the difference in accuracy? for example from 0.856 to 0.860 is it considered imortant difference?from 0.85 to 0.86?or it needs to be even bigger lets say from 0.85 to 0.90 or more?thanks!!

Either 0.856 or 0.86, sounds like what you really have, I think is fine and shows archecture convergence. Please refer to other comparable papers however for accuracy comparision. Obviously 90% is better. Things get very complicated indeed increasing the accuracy.


text:
Downloading all COI sequences from BOLD fails
I have metabarcoding sequence data (COI) from bulk animal samples (including arthropoda, nematoda, annelida, mollusca) and I want to BLAST all of these sequences. I used following command to do this: blastn -remote -db nt -query COI_all.fasta -num_alignments 2 -out COI_blasted.txt. However this results in errors similar to this post: https://www.biostars.org/p/359971/ .
These errors probably appear due to the number of sequences in my file (around 700) and
the remote connection is thus interrupted.
I found that a solution would be to use blastn with a local database and since the samples are so diverse, I would like to download ALL animal COI sequences from BOLD (or gen bank). It would not be a problem if non-animal (e.g. plant) sequences would also be included.
I think the BOLD database would be great to BLAST my sequences to. However, I'm currently struggling to find a good way to download all animal COI sequences from BOLD.
When entering COI-5P as search term on http://v4.boldsystems.org/index.php/Public_SearchTerms I receive error: Your search terms resulted in too many matching terms. Please try again with more specific search criteria.. I could likely download the sequences from all the phyla etc seperately and merge them, but I'd rather just download 1 file.
I also tried to use the API by running: wget http://v4.boldsystems.org/index.php/API_Public/sequence?marker=COI-5P. A download starts but around 3.7 MB download, it is stuck and the file I receive only contains ~5000 sequences.
UPDATE: I've contacted BOLD about the stalling behavior and this is their reply: "This issue is because of the large API request that retrieves millions of records, which our system does not handle. Please break up the search by smaller groups, such as classes."
Does anyone have a solution to download all COI sequences from BOLD in one file?
I could also download COI sequences from gen bank using the ftp://ftp.ncbi.nlm.nih.gov/blast/db/ URL, but I'm not sure which exact files I need. For 16S, 18S,.. it is obvious, but not for COI. Any suggestions?
Thanks for the help.

I couldn't help but looking a bit deeper into this. There is actually an R package that uses the BOLD API.
However, they also say that there will be timeouts when downloading large sets of data. Luckily, they suggest a workaround (downloading sequences separately for each taxonomic Class).
So you should be able to download the marker sequences using the suggested steps:
# install packages (if necessary)
install.packages("bold")
install.packages("taxize")
install.packages("seqinr")
# load packages
library(tidyverse)
library(bold)    # API interface to BOLD
library(taxize)  # for NCBI taxonomy lookup
library(seqinr)  # for FASTA output

# get class-level taxa within "Animalia" from NCBI taxonomy
taxa <- downstream("Animalia", db = "ncbi", downto = "class")
# check if taxa present in BOLD
checks <- bold_tax_name(taxa$Animalia$childtaxa_name)
taxa_bold <- checks[!is.na(checks$taxon),]$taxon

# download sequences from BOLD for each class-level taxon
sequences <- map(taxa_bold[1:3], bold_seq, marker = 'COI-5P') %>%
  flatten() %>%
  bind_rows()

# write sequences to file
write.fasta(
  sequences = as.list(sequences$sequence),
  names = as.list(sequences$id),
  nbchar = 80,
  file.out = 'coi5p.fasta')

Note that I selected the first tree taxa (taxa_bold[1:3]) for testing, so that it runs quicker. If you want to get all taxa, just remove the [1:3].
One last importand caveat: In the help for bold_seq it says that

Notes from BOLD on the marker param:  "All markers for a specimen matching the search stringwill be returned. ie. A record with COI-5P and ITS will return sequence data for both markers evenif only COI-5P was specified."You will likely end up with data with markers that you did not request - just be sure to filter thoseout as needed.

So you might want to double-check that the downloaded sequences are actually COI-5P. But if you are using them as a BLAST database, you might also be happy with COI-5P and ITS markers being present in the database (it affects the e-values, though).

I tried it and when downloading all sequences I receive following warnings: Warning messages: 1: In .f(.x[[i]], ...) :   the request timed out, see 'If a request times out' returning partial output 2: In .f(.x[[i]], ...) :   the request timed out, see 'If a request times out' returning partial output 3: In .f(.x[[i]], ...) :   the request timed out, see 'If a request times out' returning partial output
it downloaded 1012676 sequences, so quiet alot but I don't think that all are downloaded. Again, these warnings seem to be related to some kind of download limitation..

If there are timeouts one must assume that not all sequences have been downloaded. You could also try going a to a deeper taxonomic rank than class (order?) in case that some classes have so many sequences that it leads to a timeout.

I think I will use the esearch as final tool, as it is more to the point than the BOLD downloads. Thanks for your extensive help regarding the issue, for which I awarded the bounty.

Thank you @Robvh! If I had to choose for myself I would also go with GenBank as a source.

@mrhd hey so I took interest in this script. The bind_rows() command returns the following error.
Error in bind_rows(): ! Argument 1 must be a data frame or a named atomic vector. Run rlang::last_trace() to see where the error occurred.
Can you lend a helping hand after all these years?

A straightforward way of getting the sequences from GenBank is from their website using this URL:
https://www.ncbi.nlm.nih.gov/nuccore?term=COI[Gene Name]
There you can click on "Send to" -> select "Gene Features" and "FASTA Nucleotide" -> Click "Create File".
There are >3 mio hits, so the download could take a while.

Thank for your suggestion, this is certainly a possibility for the Gen Bank data. For now, I see Gen Bank data as a good alternative, but am still looking for a solution to the BOLD database.

I tried to download the sequences from BOLD with wget using the API link you provided above, and got the same stalling behaviour that you described. Maybe they apply some sort of limit per download. I'm not really familliar with BOLD.

Thanks for replicating my problem, then I know my computer is not the problem. I will contact them and ask for any limitations. I will post back their reply.

downloading all COI sequences with the method you provide in this answer also does not work. Around ~150MB download the downloading stops with error "failed-network error". I assume this is also some download limitation as we do not experience any network errors?

check the update of the original post to see the answer of BOLD regarding the stalling behavior.

I'm now downloading all COI gen bank sequences withesearch from the E-Utilities package in the Ubuntu terminal using command: esearch -db nucleotide -query "COI[GENE]" | efetch -format fasta > COI_all_genbank.fasta. This should do the trick.

That's a quite resourceful idea, @Robvh! Didn't know about that command line tool.


text:
Filtering VEP annotation file
I am filtering a VEP annotated vcf, trying to maintain just those variants classified as deleterious by SIFTand as damaging (probably or possibly included) by PolyPhen.
I am using:
filter_vep -i "$input" -o "$output" -filter "SIFT = deleterious and PolyPhen match damaging"

and I also tried:
filter_vep -i "$input" -o "$output" -filter "SIFT match deleterious and PolyPhen match damaging"

The output contains just variants with some value in these two fields, but some of the variants are classified as tolerated by SIFT and as benign by PolyPhen, so do you know why the filter is not working?
As an example: this two variants shouldn't exist in my filtered VCF, because the first one is classified as benign by PolyPhen and the second is classified as benign by both fields:
    CHROM  POS  REF  ALT  FILTER  INFO
    1       1934441 .       G       T       .       PASS    deleterious(0.01)|benign(0.138)
    1       207085100       .       G       A       .       PASS   tolerated(0.16)|benign(0.037)

As an example of a variant that has pass properly those filters:
    CHROM  POS  REF  ALT  FILTER  INFO
    2       3938449 .       T       A       .       PASS    deleterious(0.01)|probably damaging(0.67)

Those are 3 variants among those that I have in my filtered VCF, but I should be getting just those variants which are classified as deleterious and damaging by both SIFT and PolyPhen, as the third variant I've shown.

I think what you're trying to do could be achieved by adding parenthesis to the conditions and by specifying the format of the input file by using the flag --format vcf, as specified in the documentation.
The final command should look like this:
filter_vep -i "$input" --format vcf --filter "(SIFT = deleterious) and (PolyPhen match damaging)" -o "$output"

Here is another way to achieve your goal using the fuc vcf_vep command I wrote:
$ fuc vcf_vep in.vcf 'SIFT.str.contains("deleterious") and PolyPhen.str.contains("damaging")' > out.vcf

For getting help:
$ fuc vcf_vep -h
usage: fuc vcf_vep [-h] [--opposite] [--as_zero] vcf expr

This command will filter a VCF file annotated by Ensemble VEP. It essentially wraps the `pandas.DataFrame.query` method. For details on query expression, please visit the method's documentation page (https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.query.html#pandas-dataframe-query).

examples:
  $ fuc vcf_vep in.vcf 'SYMBOL == "TP53"' > out.vcf
  $ fuc vcf_vep in.vcf 'SYMBOL != "TP53"' > out.vcf
  $ fuc vcf_vep in.vcf 'SYMBOL == "TP53"' --opposite > out.vcf
  $ fuc vcf_vep in.vcf 'Consequence in ["splice_donor_variant", "stop_gained"]' > out.vcf
  $ fuc vcf_vep in.vcf '(SYMBOL == "TP53") and (Consequence.str.contains("stop_gained"))' > out.vcf
  $ fuc vcf_vep in.vcf 'gnomAD_AF < 0.001' > out.vcf
  $ fuc vcf_vep in.vcf 'gnomAD_AF < 0.001' --as_zero > out.vcf

positional arguments:
  vcf         Ensemble VEP-annotated VCF file
  expr        query expression to evaluate

optional arguments:
  -h, --help  show this help message and exit
  --opposite  use this flag to return records that don’t meet the said criteria
  --as_zero   use this flag to treat missing values as zero instead of NaN

You can also do above with Python:
from fuc import pyvcf, pyvep
vf = pyvcf.VcfFrame.from_file('in.vcf')
expr = 'SIFT.str.contains("deleterious") and PolyPhen.str.contains("damaging")'
filtered_vf = pyvep.filter_query(vf, expr)
filtered_vf.to_file('out.vcf')

If you are working under a Linux environment, then you may easily use a piped grep search. For example:
grep deleterious file.vcf |grep damaging >filtered_file.vcf

The grep output will not be a properly formed VCF.


text:
MSA (protein) with biopython or something else?
I am very new to bioinformatics (and python in general), but I would like to use python to more efficiently analyse enzymes both in terms of structure and functio, using Jupiter notebook. I would like to ask what is the best program/source code for multiple sequence alignments (amino acids) to identify conserved binding sites etc. I see that biopython has a few ways of creating alignments, but I have the impression that it is more focused on nucleotide sequences?

It works for both —but for global alignments you'd want system calls to something better (muscle say) than the inbuilt pairwise method. I don't know about MSA viewers for Jupyter, but I'd like to say that nglview is a really nice module to view your 3D protein structure in the notebook —it is despotic to set up at times, but totally worth the struggle.

Thanks, yea I got nglview to work after some time! Still working on the MSA in python, most sites seem to focus on nucleotide alignments.

I found a nice blog post about plotting MSAs in Jupyter notebooks, which I thought was relevant to this.

In my personal experience, MUSCLE is the easiest program to use in conjunction with Biopython. Biopython features a command line wrapper for this program, which makes it very easy to use. Make sure to download the appropriate MUSCLE program from drive5 and save it somewhere. E.g., if you are using Jupyter in Linux:
!wget https://www.drive5.com/muscle/downloads3.8.31/muscle3.8.31_i86linux64.tar.gz
!tar -xzvf muscle3.8.31_i86linux64.tar.gz
!cp muscle3.8.31_i86linux64 /usr/local/bin
!chmod 755 /usr/local/bin/muscle3.8.31_i86linux64

Then, you can run MUSCLE like so:
from Bio.Align.Applications import MuscleCommandline
def runMUSCLE(infile, outfile):
    muscle_exe = r"/usr/local/bin/muscle3.8.31_i86linux64" #Here is where we installed MUSCLE
    muscle_cline = MuscleCommandline(muscle_exe,
                                     input=infile,
                                     out=outfile,
                                     clwstrict=True #Output in clustal format (more visually pleasing), otherwise the output is in FASTA. Whichever you need.
                                    )
    muscle_cline()

You don't need to specify that your sequence is amino acidic, however keep in mind that the input file must be in FASTA format.


text:
Is there a built-in parser for the textual feature table that Entrez.efetch returns?
Entrez.efetch helpfully gives me the feature table of a nuccore entry.
If I understand correctly, efetch can only return the feature table in a textual format (e.g., as mentioned here).
Is there any built-in Python/Biopython function that parses this textual format of a feature table?
Here is a code snippet that uses efetch and writes the feature table to a file:
from Bio import Entrez
Entrez.email = "orenmn@gmail.com"
with Entrez.efetch(db='nuccore', id='U00096.3', rettype='ft') as handle:
    feature_table = handle.read()

with open('feature_table.txt', 'w') as f:
    f.write(feature_table)

Here are the first 20 lines of feature_table.txt:
>Feature gb|U00096.3|
190 255 gene
            gene    thrL
            gene_syn    ECK0001
            locus_tag   b0001
            db_xref ASAP:ABE-0000006
            db_xref ECOCYC:EG11277
            db_xref EcoGene:EG11277
190 255 CDS
            product thr operon leader peptide
            transl_table    11
            protein_id  gb|AAC73112.1|
            transcript_id   gnl|b0001|mrna.b0001
            db_xref UniProtKB/Swiss-Prot:P0AD86
337 2799    gene
            gene    thrA
            gene_syn    ECK0002
            gene_syn    Hs
            gene_syn    thrA1
            gene_syn    thrA2


text:
Why this ballgown error (first column of pData does not match the names of the folders) comes out
I am trying to create a ballgown object, and I would like to set a phenodata information.
For that I do
phenodata <- read.csv('l002.csv')

and my l002.csv file contains the following info
id,state
RNA-CCC_8N_S5_L002,control
RNA-CCC_8T1_S6_L002,tube
RNA-CCC_12N_S7_L002,control
RNA-CCC_12T_S8_L002,tube
RNA-CCC_13N_S105_L002,control
RNA-CCC_13T_S106_L002,tube
RNA-CCC_23N_S107_L002,control
RNA-CCC_23T_S108_L002,tube

A ballgown folder contains folders that have exactly the same names, e.g. RNA-CCC_8N_S5_L002, RNA-CCC_8T1_S6_L002, etc...
Then, I use this piece of code for bg object creation
bg = ballgown(dataDir = "ballgown", samplePattern = "RNA-CCC_", pData = phenodata )

but every time I run it, I have the same error every time:
Error in ballgown(dataDir = "ballgown", samplePattern = "RNA-CCC_", pData = phenodata) :
  first column of pData does not match the names of the folders containing the ballgown data.

What can be a problem since first column of pData matches exactly the names of the folders containing the ballgown data..
If I get rid of phenodata ballgown object creation works, but I need this phenodata for further manipulations.

Well, it looks like if you're using Windows for managing files - that's not a good idea - it sorts names of your ballgown folders in alphabetical order while you need it differently.
That is why the one should use Linux - when you're ls -ltr your ballgown folder in Linux it gives you the "right" sequence of folder names to be put inside a ballgown object using pData.
This is what helped me to get rid of this error:
In Linux enter your ballgown folder and do ls -ltr | awk '{print$9}' and you will have the right sequence that you have to put inside you .csv file using Notepad++, indicating what's what. Below's my example compared to the above mentioned initial example:
id,state
RNA-CCC_12N_S7_L002,control
RNA-CCC_12T_S8_L002,tube
RNA-CCC_13N_S105_L002,control
RNA-CCC_13T_S106_L002,tube
RNA-CCC_23N_S107_L002,control
RNA-CCC_23T_S108_L002,tube
RNA-CCC_8N_S5_L002,control
RNA-CCC_8T1_S6_L002,tube

I don't have my error anymore and after
phenodata <- read.csv('my_phenodata_file.csv')
bg = ballgown(dataDir = "ballgown", samplePattern = "RNA-CCC_", pData = phenodata)

my results are nicely wrapped
Sat May 02 22:11:57 2020
Sat May 02 22:11:57 2020: Reading linking tables
Sat May 02 22:11:59 2020: Reading intron data files
Sat May 02 22:12:04 2020: Merging intron data
Sat May 02 22:12:06 2020: Reading exon data files
Sat May 02 22:12:14 2020: Merging exon data
Sat May 02 22:12:16 2020: Reading transcript data files
Sat May 02 22:12:20 2020: Merging transcript data
Wrapping up the results
Sat May 02 22:12:22 2020

The problem is explained here as well
https://www.biostars.org/p/228389/
Hope this will help anyone

Thank you very much @'Dmitrii Trubetskoy' it certainly helps.