text:
How to obtain the standard gene name
I have a list of genes and proteins extracted from literature. How can I convert them into standard gene names.

It would be helpful to have some example of the gene/protein names that you extracted so solutions can more easily be checked for validity.

Please define "standard gene names", do you mean gene symbols?

Please [edit] your question and: i) give us a few examples of the various name formats you have, ii) show us the output you would want from those names (do you want gene name? Gene symbol? From what database? GI accession? RefSeq? Ensembl? The human name from HUGO?). What do you consider a "standard" gene name? And in what species?

Gene name
Using an Ensembl ID, Uniprot ID or NCBI ID it is much easier to look up your gene. So I would strongly reccomend using one of those as a key. Personally, I'd use Uniprot, but Ensembl is probably better if your genes are human.
One reason for this is also that different resources will use different names, especially for humans.
There are a few considerations, so I will first say how this could be done properly...
Custom Uniprot matcher with weights
Some names are used by different protein, so which one is which would depend on context, co-occurrence with another gene.
Your dataset might be extremely noisy. Therefore, you likely have assigned a confidence score to your matches, yes? It would be nice to give a score to the name based on whether it is the first name of the protein or the fourth synonym in the list.
Therefore to do a conversion properly you really ought to make your own conversion system based upon context and Uniprot data. Uniprot contains the most data on synonyms and notes and binding partners etc. The dataset/API is in XML and Biopython does not have a parser, but I have some notes on the matter and can share some code if needed as I have done something similar (e.g. JSON for human gene/protein names to Uniprot IDs).
Lazy match
The NCBI entrez API is rather flexible and good at presenting to you matches for a strange names. So I would suggest using that.

If you have a list of human gene names from various sources that may include old or alternative names, then I'd suggest using the HUGO database (www.genenames.org). From there you can get a table of approved names and symbols as well as alternative symbols* . With this you can search for your gene names in the table and assign the approved name.
I wouldn't use  Ensembl ID as some genes have multiple IDs that indicate the haplotype of the source genome, which is probably not useful for most purposes.
*Names are the full description of the protein, e.g. "B-Raf proto-oncogene, serine/threonine kinase" whereas the symbol is the short version you usually see, e.g "BRAF"

bioDBnet might have what you're looking for if the gene name is in standard format
https://biodbnet-abcc.ncifcrf.gov/db/db2db.php
For proteins Uniprot can be used to look up a protein and find it's information
https://www.uniprot.org/


text:
snakemake get environment variable for rule name and log path
When a rule fails in snakemake workflow, it returns a descriptive message in log file like following:
Error in rule step1:
    jobid: 99
    output: results/foo/step1/sample5/outdir/output1.tsv, results/foo/step1/sample5/outdir/output2.tsv
    log: logs/step1/sample5_out.log (check log file(s) for error message)
    shell:
    ...
        (exited with non-zero exit code)

Since I run workflow over many samples, I have defined a postrun.sh bash script in jobscript.sh directive that runs after snakemake {exec_job} section and notifies me via userbot in slack channel when a step failed for one or more samples.
jobscript.sh
        ########################### START SNAKEMAKE CMD ############################
    {exec_job}
        ############################ END SNAKEMAKE CMD #############################
    exitstat=$?
    export exitstat

    echo "END at $(date)"
. "${{HOME}}"/.smk_confs/postrun.sh

postrun.sh
## execute after running each of snakemake job
echo 'Running postrun.sh'
ERRMSG="MYJOB ID: $PBS_JOBID failed at $(pwd) on $(hostname) for ${USER} with exit status: ${exitstat}."

## ping to slack
ping_slack -m "$ERRMSG"
sleep 5

I am looking for a snakemake exported environment variable that I can use in postrun.sh to know which rule ID failed for one or more samples. For now, I can think of using job scheduler variable matching job name, e.g, PBS_JOBNAME, but better if I have a variable from snakemake API for rule ID or preferably, path to log file.
Thanks

You could add an exit trap to the beginning of the shell scripts inside your Snakefile, which is executed after each script is done and sends a message when there was an error:
 shell:
    """
    trap 'if [ $? -neq 0 ]; then ping_slack -m "Error in rule xxx with input {input}"; exit 1' EXIT
    ...

I believe this needs to be added to all of workflow rules. I am hoping to set a trap similar to existing onerror handler of snakemake but prefer to execute it at completion of each rule rather an entire workflow.

yes, that would be needed in each rule.

source: Stack Exchange

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

sort: votes


text:
how to generate an scRNA counts file
I have barcodes.tsv, genes.tsv and matrix.mtx file from my RNA seq data and I want to generate a scRNA counts files using R. When I tried using the readMM function, it gave me the following error:
Error: cannot allocate vector of size 185.1 Gb.

This happens even when I run it as a job. How can I generate a combined csv file using R?
Thank you.

Hi and welcome to the site! It is very hard to help without a concrete example to work on. Could you show us a few lines of a couple of your input files and then the output you require from those example lines? That way we can understand what you need better and we can also make sure any solutions we come up with actually work for your data.

Please [edit] your question to add extra information. Comments are hard to read, easy to miss and can be deleted without warning.

Also posted on biostars: https://www.biostars.org/p/395275/

Most likely you tried other things with these files and you are not reading a sparse matrix anymore. a good place to start is usually the documentation associated with the software. For cell ranger v3.0 you can find this here. https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/matrices.  If you scroll down you will see a section on how to generate a csv file. Restart with the zipped output of the filtered_feature_bc matrix. If you still get error update your post with code + error

I use the DropletUtils package from Bioconductor. Getting counts begins like this:
library(DropletUtils)
filePath <- "path_to_data/raw_feature_bc_matrix/"
singleCellExperiment <- read10xCounts(filePath, col.names = TRUE)

You probably need to load the matrices part by part in order to avoid loading a large size vector. One way is to use features like nlines or skip of the scan() function. Also, this is suggested code by 10xgenomics:
library(Matrix)
matrix_dir = "/opt/sample345/outs/filtered_feature_bc_matrix/"
barcode.path <- paste0(matrix_dir, "barcodes.tsv.gz")
features.path <- paste0(matrix_dir, "features.tsv.gz")
matrix.path <- paste0(matrix_dir, "matrix.mtx.gz")
mat <- readMM(file = matrix.path)
feature.names = read.delim(features.path,
                           header = FALSE,
                           stringsAsFactors = FALSE)
barcode.names = read.delim(barcode.path,
                           header = FALSE,
                           stringsAsFactors = FALSE)
colnames(mat) = barcode.names$V1
rownames(mat) = feature.names$V1

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


text:
samtools mpileup skipping read
I run the command:
samtools mpileup -O -s -q20-B -Q20 -f hg19.fa -r chr1:569929-569931 myFile.bam

and get:
chr1    569929  G   7   ...,,., EEEEEEE NTTVTTU 53,48,42,60,30,29,27
chr1    569930  C   6   ...,.,  EEAEEE  NTTTTU  54,49,43,31,30,28
chr1    569931  G   6   ...,.,  EEEEEE  NTTTTU  55,50,44,32,31,29

I am interested in position 569930.
I inspect the bam with:
samtools view myFile.bam chr1:569929-569931

and I see:
NS500355:NS500355:HHYN5BGXB:2:21311:2303:9255   99  chr1    569869  0   76M =   569900  107 CTCCATAACGCTCCTCATACTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGA    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE<EEEEEEEAEEEEEEE    NM:i:2  MD:Z:5C3C66 MC:Z:76M    AS:i:66 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9321,76M,1;
NS500355:NS500355:HHYN5BGXB:4:22603:25857:2799  99  chr1    569877  45  76M =   569896  95  CGCTCCTCATACTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCAC    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE    NM:i:1  MD:Z:1C74   MC:Z:76M    AS:i:74 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9329,76M,1;
NS500355:NS500355:HHYN5BGXB:1:12105:11868:12011 99  chr1    569882  51  76M =   569900  94  CTCATACTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACC    /AAAAEEEEEAEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEA    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9334,76M,1;
NS500355:NS500355:HHYN5BGXB:2:13301:2681:6114   99  chr1    569888  9   76M =   569952  140 CTAGGCATACTAACCAACACACTAACAATATACCAATGATGGAGCGATGTAACACGAGAAAGCACATACCAACGCC    AAAA/E/EE//AAEEEEE/EAAEAA//EEAEEE/E/EEEE/////A/EE/EE/EE/E//EE/EEEE/E/AEA####    NM:i:4  MD:Z:6C19C15C29G3   MC:Z:76M    AS:i:57 XS:i:52 RG:Z:../H2O XA:Z:chrM,+9340,76M,5;
NS500355:NS500355:HHYN5BGXB:2:23110:10125:9334  99  chr1    569888  51  76M =   569904  92  CTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCC    AAA/AEE6EEEE6E///E/A/AAE/E/AEEEEEEE/EE/EEEAEEEEAEE//EEEEEEAEEEAEEAE/AEE/AEE<    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9340,76M,1;
NS500355:NS500355:HHYN5BGXB:1:13209:10692:11045 83  chr1    569888  53  18S58M  =   569888  -58 CGTGTGCTCTTCCGATCTCTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAG    EE<A/EE//AE/AA///AAE///EEEEE//6///E//EE/6EE/EEEEEEEA/A/A/E/E//AEAEEEEE6A6AAA    NM:i:0  MD:Z:58 MC:Z:59M17S AS:i:58 XS:i:53 RG:Z:../H2O XA:Z:chrM,-9340,18S58M,1;
NS500355:NS500355:HHYN5BGXB:1:21304:26480:15673 83  chr1    569890  11  76M =   569796  -170    AGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCAC    AEEEEEEEAAEEEEEEEEEEEE/EEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAAAA    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,-9342,76M,1;
NS500355:NS500355:HHYN5BGXB:4:11602:21237:3028  83  chr1    569900  51  76M =   569900  -76 ACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCAC    EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAAAA    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,-9352,76M,1;
NS500355:NS500355:HHYN5BGXB:3:13504:4369:8812   99  chr1    569901  51  76M =   569901  76  CCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACC    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE6EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9353,76M,1;
NS500355:NS500355:HHYN5BGXB:2:11203:19774:18704 99  chr1    569904  11  76M =   569952  124 ACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGT    6AAAAEEEEEEEEEEEEEAEEEEEEEEEEEEEEEE6EEEEAEEEAEEEEEEEEEEEEEEEEEEAE/E6EEAEEEEA    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9356,76M,1;
NS500355:NS500355:HHYN5BGXB:2:11309:13954:17567 99  chr1    569905  11  76M =   569972  143 CACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGTC    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEE    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9357,76M,1;
NS500355:NS500355:HHYN5BGXB:1:21202:13252:14049 99  chr1    569905  11  76M =   569927  98  CACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGTC    AAAAAEE6EEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEE6EEEEEEEEEEEEEEEEE/EEEE<EEEEEEEEEEE    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9357,76M,1;
NS500355:NS500355:HHYN5BGXB:1:22311:21964:12072 83  chr1    569908  52  5S71M   =   569908  -71 GATCTACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTG    <AEEEA6EEE6EEEEEEEEEEEEEE/EEEEEEEEEEEEEEEEE/EEEEEEEEEEEEEEEEEEAAEEEAEEEAA6AA    NM:i:0  MD:Z:71 MC:Z:71M5S  AS:i:71 XS:i:66 RG:Z:../H2O XA:Z:chrM,-9360,5S71M,1;
NS500355:NS500355:HHYN5BGXB:4:12602:7858:19418  99  chr1    569909  15  76M =   569917  84  CTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGTCCAAA    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEE/AAEEEEEEEEEEEEEEE    NM:i:1  MD:Z:74G1   MC:Z:76M    AS:i:74 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9361,76M,1;
NS500355:NS500355:HHYN5BGXB:2:12205:19456:7598  83  chr1    569927  0   76M =   569907  -96 TGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGTCCAAAAAGGCCTTCGATACGGGA    EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAA/EEEEEEEEEAAAAA    NM:i:1  MD:Z:56G19  MC:Z:76M    AS:i:71 XS:i:76 RG:Z:../H2O XA:Z:chrM,-9379,76M,0;

Based on mapq and positions (respect to 569930) I see that the read with cigar=18S58M is skipped by mpileup.
WHY? The soft clip is way far from position chr1:569930

Examining the mpileup output for a single read on its own is a good way to figure out why bases are not appearing as expected.
So if you suspect that the cigar=18S58M read is the missing one (you are correct), prepare a SAM file containing just that read:
@SQ SN:chr1 LN:248956422
NS500355:NS500355:HHYN5BGXB:1:13209:10692:11045 83  chr1    569888  53  18S58M  =   569888  -58 CGTGTGCTCTTCCGATCTCTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAG    EE<A/EE//AE/AA///AAE///EEEEE//6///E//EE/6EE/EEEEEEEA/A/A/E/E//AEAEEEEE6A6AAA    NM:i:0  MD:Z:58 MC:Z:59M17S AS:i:58 XS:i:53RG:Z:../H2O  XA:Z:chrM,-9340,18S58M,1;

and run it through mpileup:
$ samtools mpileup -Q0 -B oneread.sam

This will show you that this read's base quality dips below 20 at the position you are interested in. Thus, as finswimmer suspected, your -Q20 is removing it at this position.

Why do you believe it's the read with cigar=18S58M? The mapping qualities given in sam files are different to those in the mpileup output. I guess samtools is doing some recalculation.
You are using the parameters -q20 -Q20. So you skip all reads with mapping quality <20. This will result in 7 reads got to mpileup. Furthermore you are skipping bases with base quality <20. Have a look at the base qualities for the read starting at position 569904. If I count it correct you end up with a quality of / at position 569930 which means 14. So I guess this read is dropped.
BTW: samtools mpileup is deprecated (only when creating bcf/vcf)  and should be replaced by bcftools mpileup.

Thank you, I missed the '/' quality for a 'E'. Now I understand.

The mpileup subcommand has been split between the two commands. Textual “mpileup output” stays in samtools (and always will), but VCF/BCF output has moved to bcftools. This groups the latter alongside the other VCF-related facilities and in particular means that both halves of the calling pipeline bcftools mpileup … | bcftools call … are in the same program, which simplifies updates and improvements to this pair of commands that work together.


text:
How to upload and process a data set to Galaxy instance and retrieve results programmatically?
How to upload and process (de novo assembly pipe) a data set to Galaxy instance and retrieve the results programmatically? Are there any packages beyond BioBlend?

Hi @player welcome. This question is way too general. There is a lot of expertise on de assembly, whether there is knowledge of you API is an issue. Please try and think of how to tie down your problem and describe it in more detail

You best bet is bioblend, since it's a more convenient wrapper around the Galaxy API, but if you really want to you could use that directly. I should note that large files should probably still be uploaded via ftp or sftp rather than via the http(s) interface.

Perhaps you could include in the answer how to do this to improve the question?


text:
c++ program to find length string an char*
I write under code but I can't run it when I run by g++ . erro is about  strlen.
errors:
a.cpp: In function ‘bool checker(char*, int)’:
a.cpp:185:26: error: ‘str’ is not a member of ‘std’
 185 |  for(int i = 0; i < std::str.length(word); i++)
  |                          ^~~
a.cpp: In function ‘int status(char*, int, int)’:
a.cpp:192:40: error: ‘str’ is not a member of ‘std’
   192 |  for(int i = 0; i < 4 && nw*4+i < std::str.length(pmt[pn]); i++)
  |                                        ^~~
a.cpp:199:10: error: ‘str’ is not a member of ‘std’
  199 |  if(std::str.length(pmt[pn]) > nw*4)
  |          ^~~
a.cpp: In function ‘int main()’:
a.cpp:225:22: error: ‘str’ is not a member of ‘std’
  225 |     len_words = std::str.length(pmt[i])/2;
    |                      ^~~

Isn't this a better question for stackoverflow?

how i can ask best question

@fatemehkhoramdel please format your question properly, check your English to make it easier to understand (capitalize the beginning of sentences, ...). Explain what this code it is supposed to do and why it is relevant that you need bioinformaticians to help you (there are other sites to get help that might have people with the right skills to help you, but they might not be here). Good luck

I think this is what you're looking for: error strlen not a member of std only on linux not on solaris.

it not usefull for me.i treid about it

I'm voting to close this question as off-topic because it belongs on Stack Overflow.

@AlexReynolds Respectfully, I would caution against simply closing as off-topic simply because this is a programming question. While I agree this question is not directly related to bioinformatics, the expectation of prior knowledge is going to be different for bioinformaticians and computer scientists. The SO C/C++ community is not really for beginners, and I think it would be more productive for bioinformatics to have a source for entry-level programming questions, as for them it's a tool and not the end result.

I think you're confusing C++ and C strings, which are very different. When you have a char *, that's a C-style "string". It's really just a bunch of bytes with a zero byte at the end, demarcating the end of the string. To find the length of a C string you include <cstring> (or <string.h> in C) and call strlen, like this:
/** Print every character in the string */
void print_chars_in_string(char* word) {
    for (size_t i = 0; i < strlen(word); ++i) {
        putchar(word[i]);
    }
}

strlen is actually pretty simple though, you can write your own version as an academic exercise like this:
size_t deluxe_bespoke_strlen(char* str) {
    size_t length = 0;

    while (*str++)
        ++length;

    return length;
}

Since the while-loop condition will be true until you hit the null terminator (the zero byte at the end), this is how you would implement the strlen function without a for-loop, otherwise you would need the strlen function for the strlen function.
In C++, strings are actual objects, with methods you can call. To get the length of a C++ string simply do this:
std::string my_name = "Jose";
std::cout << my_name.length() << std::endl;

If you have a C string, you can convert it to a C++ string with the built-in constructor:
const char* my_name = "Jose";
std::string my_name_cppstring = std::string(my_name);

You can also get the C version of a C++ string by calling the c_str() member function like this:
std::string sequence = "GATTACA";
size_t length = strlen(sequence.c_str());

My recommendation is to not mix C and C++. C is my favorite programming language, but that doesn't mean everyone should write everything in C. Since you're using C++ I'm guessing that's your preference, so I would suggest using C++ strings instead of C strings, just remember to include <string>.

As for your specific question, the compiler is giving you the correct message: str is not a member of std. You don't add the namespace to variables. This is how that works:
#include <cstdlib>
#include <iostream>

namespace CustomDeluxeBespokeStandardNamespace
{
    size_t custom_bespoke_strlen(char* str)
    {
        size_t length = 0;

        while (*str++)
            ++length;

        return length;
    }
}

int main(int argc, char *argv[])
{
    if (argc == 1) {
        std::cerr << "No input(s)." << std::endl;
        std::exit(EXIT_FAILURE);
    }

    while (*++argv) {
        std::cout << *argv << ": " << CustomDeluxeBespokeStandardNamespace::custom_bespoke_strlen(*argv) << std::endl;
    }
}

So if you run it like this:
testprogram horse house mouse cat

You should get this output:
horse: 5
house: 5
mouse: 5
cat: 3

So to recap, you were adding the namespace to your variables, and the compiler was giving you an error because it correctly could not find anything by that name in the std namespace. Even if it had though, C style strings do not have member functions. There are no such things as objects in C, at least not in the sense that you can define member functions for an object, so that was the second error.

As an aside, you will note all my examples using strlen both custom and actual return size_t. This is not an accident; in both C and C++ the return type of a call to strlen is size_t. You have to be careful because in your code you compare an int to a size_t. Comparing a signed value to an unsigned value is Not Good because it will work only if the signed value can be faithfully represented after a conversion to an unsigned value. If it can't though, you're in for a bad time.

ok, it did, but no results are shown when i create the file

Can you post your code?

my code have 2000 line...i think i can't post code here..

how i can post it?

Can you link to a github?

No,it's not possible.i'm iranian.github has boycotted iran:(

Damn, if the code doesn't contain any sensitive information you can email it to me at jflopezfernandez@gmail.com and I can help you over email, since the question got closed lol

Thanks, I sending you


text:
Patient-sample mapping in GSE72056 dataset
I want to use the single-cell data from the following expression profiling which concerns the RNA-seq of metastatic melanoma:
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE72056
The data was also analyzed here,
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4944528/
However I cannot find the mapping between patient and single-cell samples anywhere, neither in the paper nor in the supplementary materials. Anybody has an experience with this data?
Note that in the Supp1 (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4944528/bin/NIHMS791837-supplement-Supp1.doc) from the paper, at Table S1 (almost at the end of the file) that represents the mapping between tumor and patient, which is a one-to-one mapping. But no mapping between samples and tumors found.

No answer yet?!

If there is no answer, it is likely because no one has found an answer. Clearly you have had some trouble with this, perhaps it is a difficult question to answer.

Hey i briefly checked the publication. It is really hard to make a connection between the sequenced sample names and actual sample name. What I usually do is to write to the authors directly..

OK Thanks. I will go for it.


text:
Efficiently concatenating disjoint, ordered bigBed files
Given a set of bigBed files, each of which are disjoint to each other and in lexicographical ("sort-bed") sort order, I would like to know if there is a way to concatenate these archives efficiently.
Presuming I know the order in which to specify files (to constrain the order of intervals), is there a way to concatenate them directly, i.e., without first extracting them to bed, concatenating the extracted files, and remaking a bigBed file?

No, there's no simple way to do that without someone writing a program tailored toward it that handles the footer section and recreates the zoom levels that follow the compressed blocks of intervals. If you wanted to write such a program, modifying bigWigCat from Kent's tools would be a reasonable tack.

Seems reasonable, thanks.


text:
How do I create a filtered gene list using expression medians
Forgive the simple noob question
I have TPM data of ~50k genes (rows) across ~1k cell lines (columns). In R, I would like to output an "intermediate expression" gene list for each cell line, like:
>head(finalExpressedGeneList, n=3)

Celltype1    Celltype2    CelltypeN
gene1        gene5        gene3
gene3        gene6        gene6
etc          etc          etc

I am defining intermediate expression as "gene expression >= median value within each cell line." I have already replaced all 0 values with NA using:
>data[data == 0] <- NA

I then got medians for each column using:
>colMed <- apply(data,2, FUN = median, na.rm = TRUE)

Now I need to somehow filter the rows from each column and copy the associated gene into a new matrix. I started building the code below:
>if (data[,1]>=colMed [1])
   {
   finalExpressedGeneList <- row.names(data)
   }

But I realized that even if I got it working, I would still have to iterate over all of my columns. I'm sure that there is a cleaner way to do this.
EDIT (adding example data)
#create example data
>j = 5
>M <- matrix(NA,j,j)
>M <- sapply(1:j, function(i) `length<-`(1:(j-i+1), j))
>rownames(M) <- c("gene1", "gene2", "gene3", "gene4", "gene5")
>colnames(M) <- c("cellLine1", "cellLine2", "cellLine3", "cellLine4", "cellLine5")
>view(M)

       cellLine1   cellLine2   cellLine3   cellLine4   cellLine5
gene1          1           1           1           1           1
gene2          2           2           2           2          NA
gene3          3           3           3          NA          NA
gene4          4           4          NA          NA          NA
gene5          5          NA          NA          NA          NA

I would like to get an output of:

cellLine1   cellLine2   cellLine3   cellLine4   cellLine5
    gene3       gene3       gene2       gene2       gene1
    gene4       gene4       gene3
    gene5

because the median of cellLine1 has three genes greater than or equal to the median of that column (3), cellLine2 has two genes that are greater than or equal to the median of that column (2.5), etc...

This below should work. It returns a list,
## example data
j = 5
M <- matrix(NA,j,j)
M <- sapply(1:j, function(i) `length<-`(1:(j-i+1), j))
rownames(M) <- c("gene1", "gene2", "gene3", "gene4", "gene5")
colnames(M) <- c("cellLine1", "cellLine2", "cellLine3", "cellLine4", "cellLine5")
#####
genelist = apply(M,2,function(i)rownames(M)[which(i>=median(i,na.rm=T))])

####
#you can call the genes for cellline 1 using:
genelist$cellLine1

Hope I got it correct!


text:
Create a Grouped Barchart
I’m trying to create a Group Chart like this:

I have 7 classes for each of 4 species, not 5 like in the diagram. I want to build i diagram to show the distributions of the classes in each species
Here’s example provided online:
  library(ggplot2)

    # create a dataset
    specie <- c(rep("sorgho" , 3) , rep("poacee" , 3) , rep("banana" , 3) , rep("triticum" , 3) )
    condition <- rep(c("normal" , "stress" , "Nitrogen") , 4)
    value <- abs(rnorm(12 , 0 , 15))
    data <- data.frame(specie,condition,value)
    # Grouped
    ggplot(data, aes(fill=condition, y=value, x=specie)) +
        geom_bar(position="dodge", stat="identity")

------------------------------------------
    #Here’s code:

    > Species <- c(rep("Cames1" , 5) , rep("Came2" , 5) , rep("Cams3" , 5) , rep("cams4" , 5) )
    > Class <- rep(c("EVV" , "MaLR" , "psy" , "Nargo" , "Pao" , "Unknown" , "Undefiend") , 8)
    > value <- abs(rnorm(20 , 0 , 15 ))
    data <- read.table(text="Class Cames1 Came2 Cams3 cams4
    + EVV 2194 5435 2706 2902
    + MaLR 503 108 184 231
    + psy 284 187 413 295
    + Nargo 274 1977 413 2195
    + Pao 66 29 78 97
    + Unknown 286 283 273 377
    + Undefined 76 77 80 95")
   > ggplot(data, aes(fill=Class, y=7, x=Species)) + geom_bar(position="dodge", stat="identity")

Error: Aesthetics must be either length 1 or the same as the data (7): x, fill

i am able to build basics chart, but to group them all into one like the example provided challenging for me , could someone help out ?

I'm voting to close this question as off-topic because it is not about bioinformatics.  It is general plotting

I have noticed a few errors:
1) You define a bunch of variables (Species, Class, ...) but then instead of creating your data frame data with these you create it by reading from text. I don't think this is good practice, you should do so like in your example code.
2) With the way you define, your data frame does not contain columns like Class or Species but columns V1 to V5, you should use argument header = TRUE in your read.table() call.
3) Even when you read your data as in step 2, you would need to transform your data table into "long format", please see below.
4) In your ggplot() call, you assign 7 to your y aesthetic, this does not make sense. You should either provide a vector of length equal to your other aesthetics (x and fill) or better a column name of values of interest.
# This is your code minus the redundant variables that you have created but are
# not used

library(tidyr)
library(ggplot2)

# Pay attention to the header argument
data <- read.table(
  text = "Class Cames1 Came2 Cams3 cams4
     EVV 2194 5435 2706 2902
     MaLR 503 108 184 231
     psy 284 187 413 295
     Nargo 274 1977 413 2195
     Pao 66 29 78 97
     Unknown 286 283 273 377
     Undefined 76 77 80 95",
  header = TRUE
)

# data in wide format
> data
      Class Cames1 Came2 Cams3 cams4
1       EVV   2194  5435  2706  2902
2      MaLR    503   108   184   231
3       psy    284   187   413   295
4     Nargo    274  1977   413  2195
5       Pao     66    29    78    97
6   Unknown    286   283   273   377
7 Undefined     76    77    80    95

# Making the data wider: Keeping Class column as it is but condensing the data
# within the rest of the columns into two newly created columns, Species_ and
# Value_

data_long <- gather(data,
                    key = "Species_",
                    value = "Value_",
                    -Class)

# First few rows of your data in long format
> head(data_long)
    Class Species_ Value_
1     EVV   Cames1   2194
2    MaLR   Cames1    503
3     psy   Cames1    284
4   Nargo   Cames1    274
5     Pao   Cames1     66
6 Unknown   Cames1    286

# x and y aesthetics are the newly created columns Species_ and Value_
ggplot(data_long, aes(fill=Class, y=Value_, x=Species_)) + geom_bar(position="dodge", stat="identity")

oh thank you a lot i really appreciate you help, one more thing am i able to edit the colours or perhaps make it looks better than this , to report it in a journal

ggplot2 offers many options for customized plots but the easiest would be to use a built-in ggplot2 theme (just do "+ theme_void() after the geom_* line for example). Here is the link for all ggplot2 themes: https://ggplot2.tidyverse.org/reference/ggtheme.html . Moreover the ggthemes package offers many more themes.

ggplot(data_long, aes(fill=Class, y=Value_, x=Species_)) + geom_bar(HERE i ADD MY THEM????")

Please go through the link I have provided in my previous comment, if you do so you will see that a basic ggplot2 plot consists of a call to ggplot() function and this is followed by a call to your geom function (for example geom_bar() ) and then you can add other other ggplot2 functions (such as theme_void(), etc). All of these functions are piped with +.

dear sorry to disturb again but im wondering if im able to change the font style to italic for the species name ( Came2 , Cames3 ...etc)


text:
Inferring a phylogenetic tree from BLASTn
I am trying to infer a phylogenetic tree from a Blastn output and from what I have understood, what I should do is 1) extract the alignments and re-align them using Muscle, then 2) feed the .aln file into PhyML.
Is this correct?
Or should I be extracting the hsps from each alignment as suggested here https://biopython.org/wiki/Phylo? I don't quite understand why the hsps are extracted here, shouldn't the whole sequences be compared in the multiple alignment?

If you look into the get_seqrecs() function, you will see that this function is not extracting hsps to use in the following steps but the hsps are used in order to filter out low score BLAST alignments. The function accepts a blast alignment result and a threshold, and if the hsp score is above the supplied threshold the whole sequence (not just the hsp) would be the output. Your intuition is correct, the whole sequence should go in the MSA step.

You are right, my mistake, thanks!

Firsly an .aln file feed into PhylML is fine, unless you can program in Python.
Precisely what level of tree building you need depends on its end purpose. Here's past example of a similar question here using NCBI Blast.
Building trees (nj and parsimony) is one of the features of NCBI's Blast. Here's how to do it

Go to Blast here, https://blast.ncbi.nlm.nih.gov/Blast.cgi
Enter your sequence into the box (it doesn't accept PDB codes alone)
Enter the protein database - when I first did this calculation I used SwissProt, thinking there would be alot of sequences - I then used "nr"
Under the algorithm parameters enter "50" (default is too many)
Hit "Blast"
Once the search is complete at the top of the page are the hyperlinks: "Other reports: Search Summary [Taxonomy reports] [Distance tree of results]"
Click on "Distance tree of results"
This will automatically aligning your sequences and produce, in this case a parsimony based tree, but there is also the option of a nj tree (recommended)
Click "Tool", "Download", "PDF" ...

Notes

If you want probabilistic trees, PhylML is certainly one of them,  Blast remains quite central to building the alignment.
With regards building partial genome trees, this will depend on what organism you are studing
With regards Blast HSPs being treated individually, again it depends on what you question is and generally refers to advanced phylogenetics


text:
How can I find out if my gene of interest is duplicated?
To identify homologous genes, I performed tBLASTn analyses against my genome of interest resulting in sequences of a closely related species. Using this method, I found two (and in some cases up to five copies) of my different genes of interest in the target genome.
How can I determine if these gene copies correspond to true duplication events?
I aligned the newly found sequences to the reference genes and I found good alignments in more than one region (20% - 50% identity), but I am unsure if this is enough to consider them as duplicates.
Is there a threshold to help me determine if they are indeed duplications? How can I distinguish between a true duplication event and a pseudogene or a paralogous gene? The target genome belongs to a polyploid organism, although the genome assembly type is listed as haploid.
Any help will be greatly appreciated.

Thank you all very much for your responses, I am now beginning to understand how to find out if my genes are duplicated.
I performed some of the analyses that you suggested and I would very much appreciate your input on the results.
First, I followed @terdon's advice and I took the best hit for one of the query (source) genes in the target genome and I used this sequence as a query against the query (source) genome, i.e., I performed a reciprocal BLAST analysis (in my case it was a reciprocal BLASTp analysis). The first hit was my original query gene, which had an E-value of 0.0. For other sequences the reciprocal BLASTp hits were more variable, ranging from 3e-42 to 2e-154, but the first hit was always the query sequence. I also looked at the protein domains and they are conserved, i.e., both the query and the target sequences have the same domains. As for the exonic structure, it is fairly similar, for e.g., the query has 12 exons while the target sequence has 10. Exons are similar in length. So, in summary, when the potential duplicates from my target genome are blasted against the original source genome, the best hit is the original query sequence.
As for @zhanf's suggestion, I did a basic phylogenetic analysis using Clustal Omega (neighbour-joining tree without distance corrections) using four sequences: the original query sequence and the three gene copies (seq_5-1, seq_5-2 and seq_5-3) that I found in the target genome. This is the result:

I also compared each of the sequences against the original query:

seq_5-1 vs query: 34% identity (query length 780 aa, seq_5-1 length 400 aa)
seq_5-2 vs query: 42% identity (query length 780 aa, seq_5-2 length 580 aa)
seq_5-3 vs query: 40% identity (query length 780 aa, seq_5-3 length 490 aa)

So in summary, these sequences would appear to be paralogues, right?
Again, thank you very much for your help, I greatly appreciate it!
Leah

I am not sure if the phrase "true duplication" would be correct, a pseudogene could very well be the result of a duplication prior to losing its function. If the question would be if some of the newly identified orthogonal genes are pseudogenes, you might want to check the length of the resulting ORFs from the orthogonal genes, I would expect strikingly shorter ORFs for pseudogenes (depending on the time of the duplication and the selection pressure on the species).

Homologues vs. paralogues ... you've got to be quite serious about the question for someone else to describe this in detail. Are you familiar with concerted evolution?

I don't have enough reputation to comment, but you could also look at where they are located in the genome. duplicated genes are often arranged in a tandem array (right next to each other), although not always. If the genes are right next to each other in the genome, and they share a high sequence identity, they are likely duplicate genes. A much more specific way to tell, would be to make a phylogeny of the all of the supposed gene copies from the species you know. Depending on how they group together you may be able to tell if they are paralagous or separate genes.

Just to mention its easy to solve, but we'd need your cooperation. Understanding the gene and the organism would be a help.

Thanks for your answer @Leah. The BLAST hits may not be the same gene, because thats alot of distance at 34-40% identity. The "duplication" event would not be what I would call a paralogue and I suspect the genes have different functions. The duplication event could have occurred (probably did occur) before the species has evolved, in which case I personally wouldn't consider that paralogy, but it depends at what taxonomic level you are looking at (would have to be pretty major taxonomic grouping).

The way to resolve this is complicated and if you are happy with your analysis I would simply say seq_5-1 to 5-3 are the closest proteins measured by their percentage identity and may represent ancestral duplication events. Paralogues are heavily associated with concerted evolution and these proteins will not show mosaic gene recombination. I don't know if that means much to you, but I would steer clear of calling these paralogues.

@MichaelG., Leah, please don't post answers answering each other. These sites are strictly Q&A and don't work like forums. There should be one question and all answers should be answering that question. I tried to tidy up a little and included one of Leah's answers in the question and converted some of Michael's  answers to comments, but I had to delete some. Please take a moment and visit ou [help] or take the [tour] to understand how this site works.

You might want to have a look at What is the difference between orthologs, paralogs and homologs?. Paralogs are duplicates by definition, and whether a gene is a pseudogene or still active has no bearing on whether it comes from a duplication.
In any case, to figure out what's going in in your target genome, you should use the target's genes. So, first get the best hit for each of your query genes in the target genome. Then, to investigate its evolutionary history, use that sequence (the one from the target genome) as  query against the target genome. That will let you assess how similar the target's hits are to one another and that's how you can attempt to figure out their evolutionary history.
Note that 20% identity is very low for closely related species. While there is no golden threshold, I would expect recent duplication events (and in closely related species, a duplication event that is in one genome and not the other will be recent) to produce much more similar sequences. In any case, you can't do this based on sequence identity alone.
Are the protein domains of the gene products conserved? Is the exonic structure of the gene conserved? You should also try taking the putative duplicate of geneA from genome2 and blast it against genome1 (your original source genome). Is the best hit the geneA from genome1, the query you used to find the putative duplicate in genome2? Or is there a better hit for it against another gene?
These are all questions to investigate and take into account when reconstructing the evolutionary history. You can't just look at identity and use some threshold to decide this sort of thing.
If you blast the potential dupicate from your target genome against your original source genome, is the best hit the duplicate?

I don't have enough reputation to comment but if a sequence is duplicated then there has obviously been a duplication event. Doesn't matter how scientists call it now.It can be older or more recent but it's yet a duplication.


text:
Finding drugs based on gene expression data
I've got a list of genes and I want to assess their association with drugs based on gene expression data sets.
A comprehensive tutorial would also be helpful.
Any help would be appreciated.

I think the question needs elaborating, what drugs, what genes, what species, what disease are the drugs targeting. Without this I'm not sure the question is viable

It would be great if you can provide more information about the gene expression data set. If you have response to single drug, the general approach is to find the deferentially expressed genes with tools like DESeq2, or make a heat map based on the expression of your listed genes to see the drug effect before and after the treatment. But, if you want to study the effects of combination of drugs it is necessary to consider the additive and synergistic effect of drug combinations.

I agree with this


text:
Are there databases to annotate non-coding mutations?
I would like to ask if there are some known databases that allows to annotate non-coding mutations in a given region for the human genome?
Preferably, if the database is downloadable so I am able to automatize the processes of annotation.

What do you mean with non-coding mutations? SNPs? Check dbSNP.

What sort of annotations are you looking to include? Ensembl's VEP is great for annotation, and can include many different annotations applicable to both coding and non coding regions, e.g. minor allele frequency, conservation, pathogenicity predictions through CADD (via a plugin) and many other things

@benn yes SNPs, I know about dbSNP, I wanted to know if there are other known databases !

@JenG I have VEP installed on my cluster ! but CADD also comes with dbSNP plugin, but thanks I will have a look on it

We could give you better answers if you explained what sort of information you are looking for. What kind of "annotating" do you need?

Agreed. dbSNP can easily be downloaded and integrated into an automated annotation pipeline. Please provide more details about what precisely you're looking for.

Sorry for the late reply. I am looking to annotate non-coding mutations (introns, UTRs, etc... but also exons if possible) coming into active enhancers of the human genome of the pancreatic islets. For example, if a mutation comes into an active enhancer in an intron, what kind of predicted impact would that mutation do ?

VarSome allows you to annotate any variant, including non coding ones. While ACMG classification will not be available if the variant isn't in a gene's transcript (the ACMG criteria are only applicable to variants in transcribed regions), non-coding variants in transcripts will be fully annotated. For example, if you look up variant chr17-41197559-G-T which is a non-coding variant found in the UTR of one of the transcripts of BRCA1, you will get a lot of information (probably best to click on the link, the image has a lot of data):

For non-coding variants that don't fall in any transcript, you can still get frequency data, when available, as well as DANN pathogenicity predictions, GERP score, an overview of the genomic regions and other variants and/or srtuctural variants that can be found around the region, and any community annotations that may have been submitted by other users.
There is also an API you can use to automate the annotation, and VarSome Pro that lets you annotate an entire VCF file, but those are both paid services.

Disclaimer: I work for the company behind VarSome.

This looks interesting but yikes I don't think it would be much of use in my case since automating the annotation is in my needs. Thanks though.

@Law yeah, the API and/or the Pro version where you can send VCFs would probably be helpful for you, but those aren't free. I thought I'd mention them in case you do have a budget.


text:
alternatives to PHASTER phage annotation server
I would like to use PHASTER (http://phaster.ca) from the University of Alberta, but it's currently down for maintenance, as is its predecessor, PHAST.
What alternatives to PHASTER are there for searching and annotating phage sequences within assembled genomes?
EDIT:
I have seen that PROKKA has some entries mentioning the possibility of using PHASER (currently discarded as an option) but also people using their own HMMs to annotate phage sequences using PROKKA, e.g.:
https://github.com/tseemann/prokka/issues/375
Any ideas?

Phigaro https://github.com/bobeobibo/phigaro, Virsorter

Thanks @Chris_Rands , I had a look at phigaro, and I like the fact that the pVOG families give an annotation for things like tail fiber proteins and the rest. Useful.

The answer is to use phastaf, a recently developed software tool developed by Torsten Seemann (the same developer as prokka), which comes with a packaged phaster.ca db.  Get it here

Thanks, I will have a look.


text:
Adding an attribute to GFF3 file
I failed to add Note=Gene description to mRNA attribute with the below code:
#!/usr/bin/python3
import click
import gffutils
import gffutils.gffwriter as gffwriter

@click.command()
@click.option('--gff3', help="Provide GFF3 file", required=True)
@click.option('--keep', help="Keep GFF3 file", required=True)
@click.option('--reject', help="Reject GFF3 file", required=True)
def run(gff3, keep, reject):
    #db = gffutils.create_db(gff3, dbfn='data/test.db', force=True, keep_order=True, merge_strategy='merge',
    #                        sort_attribute_values=True)

    db = gffutils.FeatureDB('data/test.db', keep_order=True)
    for gene in db.features_of_type('gene', order_by='start'):
        print(gene)
        for i in db.children(gene, featuretype='mRNA', order_by='start'):
            print(i)
            print(i.id)
            print(i.attributes['Note']=["Gene description"])
            db.update([i])

I got the following error:
File "/projects/bioinf-scripts/mod_braker_gff3.py", line 20
    print(i.attributes['Note']=["Gene description"])
         ^
SyntaxError: keyword can't be an expression

What did I miss?
Thank you in advance

Please [edit] your question and add a few lines of your input gff file that we can use for testing.

But why are you printing it? Shouldn't you first add the attribute (i.attributes['Note']=["Gene description"]) and then print it print(i.attributes['Note'])? I don't know gffutils, but it seems strange to print something before setting it.

Thank you it works now.

@user977828 Please post the corrected code as an answer.

You are printing the attribute before setting it, that's why it's complaining. So instead of this:
    print(i.attributes['Note']=["Gene description"])

Set it first and then print:
i.attributes['Note']=["Gene description"]
print(i)


text:
Modifing a GFF3 file and writting to a new file
I modified the GFF3 file with the below code:
#!/usr/bin/python3
import click
import gffutils
import gffutils.gffwriter as gffwriter

@click.command()
@click.option('--gff3', help="Provide GFF3 file", required=True)
@click.option('--keep', help="Keep GFF3 file", required=True)
@click.option('--reject', help="Reject GFF3 file", required=True)
def run(gff3, keep, reject):
    #db = gffutils.create_db(gff3, dbfn='data/test.db', force=True, keep_order=True, merge_strategy='merge',
    #                        sort_attribute_values=True)

    out_keep = gffwriter.GFFWriter(keep)
    db = gffutils.FeatureDB('data/test.db', keep_order=True)
    for gene in db.features_of_type('gene', order_by='start'):
        print(gene)
        for i in db.children(gene, featuretype='mRNA', order_by='start'):
            print(i)
            i.attributes['Note'] = ["Gene description"]
            print(i)
            out_keep.write_rec(gene)

Unfortunately, the above code saved only this:
#GFF3 file (created by gffutils (v0.9) on 2019-08-28 22:14:13)
NbV1Ch08        AUGUSTUS        gene    7015    29794   0.01    -       .       ID=g1;
NbV1Ch08        AUGUSTUS        gene    60876   63944   0.03    +       .       ID=g2;
NbV1Ch08        AUGUSTUS        gene    64722   65524   0.32    -       .       ID=g3;
NbV1Ch08        AUGUSTUS        gene    66706   67484   0.12    -       .       ID=g4;
NbV1Ch08        AUGUSTUS        gene    115451  122754  0.16    -       .       ID=g5;
NbV1Ch08        AUGUSTUS        gene    147786  155025  0.24    +       .       ID=g6;
NbV1Ch08        AUGUSTUS        gene    155152  161104  0.28    -       .       ID=g7;
NbV1Ch08        AUGUSTUS        gene    237126  238098  0.05    +       .       ID=g8;
NbV1Ch08        AUGUSTUS        gene    265796  297830  0.82    +       .       ID=g9;
NbV1Ch08        AUGUSTUS        gene    265796  297830  0.82    +       .       ID=g9;

Rather than:
NbV1Ch08        AUGUSTUS        gene    7015    29794   0.01    -       .       ID=g1;
NbV1Ch08        AUGUSTUS        mRNA    7015    29794   0.01    -       .       ID=g1.t1;Parent=g1;Note=Gene description;
NbV1Ch08        AUGUSTUS        transcription_end_site  7015    7015    .       -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        three_prime_utr 7015    8531    0.2     -       .       ID=g1.t1.3UTR1;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    7015    8747    .       -       .       ID=g1.t1.exon1;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        stop_codon      8532    8534    .       -       0       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     8532    8747    0.31    -       0       ID=g1.t1.CDS1;Parent=g1.t1
NbV1Ch08        AUGUSTUS        intron  8748    9191    0.49    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     9192    9342    0.66    -       1       ID=g1.t1.CDS2;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    9192    9342    .       -       .       ID=g1.t1.exon2;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  9343    9915    0.58    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     9916    10006   0.71    -       2       ID=g1.t1.CDS3;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    9916    10006   .       -       .       ID=g1.t1.exon3;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  10007   10101   0.74    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     10102   10201   0.78    -       0       ID=g1.t1.CDS4;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    10102   10201   .       -       .       ID=g1.t1.exon4;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  10202   10712   0.8     -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     10713   11107   0.11    -       2       ID=g1.t1.CDS5;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    10713   11107   .       -       .       ID=g1.t1.exon5;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  11108   11569   0.07    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     11570   12151   0.09    -       2       ID=g1.t1.CDS6;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    11570   12151   .       -       .       ID=g1.t1.exon6;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  12152   12588   0.34    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     12589   12717   0.39    -       2       ID=g1.t1.CDS7;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    12589   12717   .       -       .       ID=g1.t1.exon7;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  12718   12789   0.42    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     12790   13075   0.39    -       0       ID=g1.t1.CDS8;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    12790   13075   .       -       .       ID=g1.t1.exon8;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  13076   14832   0.51    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     14833   15009   0.39    -       0       ID=g1.t1.CDS9;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    14833   15009   .       -       .       ID=g1.t1.exon9;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  15010   15278   0.59    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     15279   15415   0.56    -       2       ID=g1.t1.CDS10;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    15279   15415   .       -       .       ID=g1.t1.exon10;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  15416   15487   0.58    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     15488   15612   0.96    -       1       ID=g1.t1.CDS11;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    15488   15612   .       -       .       ID=g1.t1.exon11;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  15613   15706   0.96    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     15707   15957   0.98    -       0       ID=g1.t1.CDS12;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    15707   15958   .       -       .       ID=g1.t1.exon12;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        start_codon     15955   15957   .       -       0       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        five_prime_utr  15958   15958   0.99    -       .       ID=g1.t1.5UTR1;Parent=g1.t1
NbV1Ch08        AUGUSTUS        five_prime_utr  27458   28250   0.37    -       .       ID=g1.t1.5UTR2;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    27458   28250   .       -       .       ID=g1.t1.exon13;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        five_prime_utr  29272   29794   0.08    -       .       ID=g1.t1.5UTR3;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    29272   29794   .       -       .       ID=g1.t1.exon14;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        transcription_start_site        29794   29794   .       -       .       Parent=g1.t1;

What did I miss?
Thank you in advance,
Update This is base on @daniel-standage code.
import click
import tag
def add_description(instream):
    for feature in instream:
        for mrna in tag.select.features(feature, type='mRNA', traverse=True):
            mrna.add_attribute('Note', 'Gene description')

        yield feature

@click.command()
@click.option('--gff3', help="Provide GFF3 file", required=True)
@click.option('--keep', help="Keep GFF3 file", required=True)
@click.option('--reject', help="Reject GFF3 file", required=True)
def run(gff3, keep, reject):
    reader = tag.GFF3Reader(infilename=gff3)
    desc_stream = add_description(reader)
    writer = tag.GFF3Writer(desc_stream)
    writer.write()

if __name__ == '__main__':
    run()

which caused this error:
mod_braker_gff3.py --gff3 data/NbV1Ch08-augustus.hints_utr.gff3 --keep data/keep.gff3 --reject data/reject.gff3
Traceback (most recent call last):
  File "/Users/lorencm/projects/bioinf-scripts/mod_braker_gff3.py", line 26, in <module>
    run()
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/click/core.py", line 764, in __call__
    return self.main(*args, **kwargs)
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/click/core.py", line 717, in main
    rv = self.invoke(ctx)
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/click/core.py", line 956, in invoke
    return ctx.invoke(self.callback, **ctx.params)
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/click/core.py", line 555, in invoke
    return callback(*args, **kwargs)
  File "/Users/lorencm/projects/bioinf-scripts/mod_braker_gff3.py", line 22, in run
    writer.write()
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/tag/writer.py", line 57, in write
    for entry in self._instream:
  File "/Users/lorencm/projects/bioinf-scripts/mod_braker_gff3.py", line 8, in add_description
    for mrna in tag.select.features(feature, type='mRNA', traverse=True):
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/tag/select.py", line 24, in features
    for feature in entry_type_filter(entrystream, tag.Feature):
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/tag/select.py", line 99, in entry_type_filter
    for entry in entrystream:
TypeError: 'Directive' object is not iterable

Process finished with exit code 1

As haci, said you're filtering for genes. When you have this sort of issue, the first step should always be to print everything. For example, if you tried print db.children(gene, featuretype='mRNA', order_by='start'), you'd presumably see that it is empty.

If you cross post your question please leave a note, so that other can see what is already answered!

I guess you want to take the features for a gene?
As @haci and @terdon pointed out, it's a problem with the selected. Below I included an example, with mouse_extra_comma.gff3 obtained from example file provided by the package gffutils.
#gff.py
#!/usr/bin/python3

###
def addNote(x):
    x.attributes['Note'] = ["Gene description"]
    return(x)

###

import click
import gffutils
import gffutils.gffwriter as gffwriter

fn = gffutils.example_filename('mouse_extra_comma.gff3')
db = gffutils.create_db(fn, dbfn='test.db', force=True, keep_order=True,
    merge_strategy='merge', sort_attribute_values=True)

db = gffutils.FeatureDB('./test.db', keep_order=True)

## here's to test so you can see how it works
gene = next(db.features_of_type('gene', order_by='start'))
items =[addNote(i) for i in db.children(gene, featuretype=['exon','CDS'],order_by='start')]

items[0:5]
#[<Feature mRNA (chr17:6797760-6818159[+]) at 0x103659c50>,
# <Feature exon (chr17:6797760-6797769[+]) at 0x1037cee10>,
# <Feature five_prime_UTR (chr17:6797760-6797769[+]) at 0x103660898>,
# <Feature exon (chr17:6806513-6806553[+]) at 0x103659f98>,
# <Feature five_prime_UTR (chr17:6806513-6806526[+]) at 0x1036595f8>]
#####
## here's following your... code...
for gene in db.features_of_type('gene', order_by='start'):
    print(gene)
    for i in db.children(gene,order_by='start'):
        print(addNote(i))

Once you define the feature you want, should be ok

Thank you but I do not know how to do it?

Could you elaborate on which step are you stuck? Would be glad to guide you along

Where could I add something like `mrna.attributes['Note'] = "Gene description"?

You would write a function. I have edited my code above, although i am still a bit confused about what you need. If I get you correct, you would like to:

For every gene, print the gene first
Then print all features (mRNA, exon, CDS) associated with the gene, and attached a gene description note
Write to this a file

Is this correct? If so, could you please formulate the question in this way.

Thank you for the code. I need only to append to the original GFF3 file in the mRNA attribute Note=Gene description. For example  from NbV1Ch08 AUGUSTUS mRNA 22418839 22429909 0.17 + . ID=g1168.t1;Parent=g1168; to  NbV1Ch08 AUGUSTUS mRNA 22418839 22429909 0.17 + . ID=g1168.t1;Parent=g1168;Note=Gene description;. Next I would like to save the changes.

You might want to check this line out:
for gene in db.features_of_type('gene', order_by='start'):
You are specifically filtering for genes hence other types of annotations are filtered out.

Thank you but I do not know how to do it?

Please see this about the features_of_type method: https://pythonhosted.org/gffutils/autodocs/gffutils.FeatureDB.features_of_type.html . It basically tells to provide a string or tuple of feature(s) that you would like to retrieve, in your case all the features in the original file. The easiest would be take the unique values in the "feature" column of your input and feed it into the "features_of_type" method.

Where could I add something like `mrna.attributes['Note'] = "Gene description"?

If you're willing to give another GFF3 processing library a shot, I'd recommend checking out the tag library. The following code should add the note attribute to mRNA records as you're attempting to do in you're original code, but will also print out all other records.
import sys
import tag

def add_description(instream):
    for record in instream:
        if not isinstance(record, tag.Feature):
            yield record
            continue
        for mrna in tag.select.features(record, type='mRNA', traverse=True):
            mrna.add_attribute('Note', 'Gene description')
        yield feature

reader = tag.GFF3Reader(instream=sys.stdin)
desc_stream = add_description(reader)
writer = tag.GFF3Writer(desc_stream)
writer.retainids = True
writer.write()

I didn't use click or argparse to put any CLI on the script, so unless you want to add those just use STDIN and STDOUT for I/O.
$ python3 myscript.py < mydata.gff3 > mydatawithdescription.gff3

P.S. The add_description can be simplified if you don't care about retaining comments, ##sequence-region pragmas, and so on.
def add_description(instream):
    for feature in tag.select.features(instream):
        for mrna in tag.select.features(record, type='mRNA', traverse=True):
            mrna.add_attribute('Note', 'Gene description')
        yield feature

Thank you for your solution. Unfortunately, I ran into TypeError: 'Directive' object is not iterable as shown in my updated question. What did I miss?

@user977828 I wrote the initial code with the (bad) assumption that the input data didn't include any comments, pragmas, directives, etc. I've updated my answer to show how to handle input that includes these types of records, or below that is a simpler solution that ignores them.

Thank you it works.

How is it possible to avoid to change mRNA's ID?

At the moment, there is no way to change tag's default behavior: it rewrites all ID attributes, since according to the spec these attributes are only guaranteed to be consistent and unique within a single file, not across files.

Revisiting this a few months later: as of January, tag.GFF3Writer objects now have a .retainids attribute to keep the original ID values. Example updated reflect.


text:
Create heat map that groups genes by expression within cluster in Seurat
I have a set of cells that I am performing Drop-seq on to look at cell expression.  Among my heat maps for gene expression I want to be able to graph them similar to the graph below:

Where the cells are sorted by cluster on the left axis, and have the genes across the bottom.  However I also want to cluster the genes by the expression within the cluster (like how this graph does).  Is that possible in Seurat?
Currently I have created a heat map that looks like this:

It has the clusters listed on the left axis, and small sub-selection of genes on the bottom.  However I want to be able to order those genes by their expression within those clusters, without having to manually order the genes when I insert them into the "feature" argument.
Here is a sample of the code I am currently using to plot the graph:
(Note: I have an object "Combined" that has the non-averaged cell expression, but for now I do want to use average expression)
png("ExpressionHeatmap%03d.png", res=600, width=14, height=7, units="in")
top10 <- Combined.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
avg <- AverageExpression(Combined, features = NULL, add.ident = NULL, return.seurat = TRUE, verbose = TRUE)
Combined.markers %>% group_by(cluster) %>% top_n(2000, avg_logFC) -> top

cc.features <- readLines(con = "/home/alex/lab/drugres/DropSeqFiles/regev_lab_cell_cycle_genes.txt")
s.features <- cc.features[1:43]
g2m.features <- cc.features[44:97]
CCS <- CellCycleScoring(object = Combined, s.features = s.features, g2m.features = g2m.features, set.ident = TRUE)
DoHeatmap(avg, features = c("FTH1P3", "MET", "CYR61", "NFKBIA", "RIPK2", "IFI16", "PLK2", "TNFSF10", "CEBPD", "PNRC1"), group.by = "ident", size = 3, angle = 0, combine = TRUE, draw.lines = FALSE) + theme(axis.text.x = element_text(size = 7), axis.line = element_line(colour = "#ffffff")) + scale_fill_gradient2(low = '#1000ff', mid = "#ffffff", high = '#aa0101', space = "Lab", na.value = "#ffffff", midpoint = 0, guide = "colourbar", aesthetics = "fill") + coord_flip(expand = TRUE, clip = "on")

I tried to use the "FindAllMarkers", "group_by" and "top_n" functions to sort the gene data, and while it helped it still didn't fully sort the genes the way I wanted.
markers <- FindAllMarkers(Combined, features = c(genes-i-am-looking-for), only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25, test.use = "biomod")
geneorder <- markers %>% group_by(cluster) %>% top_n(n = number-of-genes, wt = avg_logFC)

I would then replace the 'features' argument in the DoHeatMap function with the 'geneorder' object.
I know that doing a top10 <- combined.markers style function can get me the effect I'm looking for, but I want to do it, with specifically selected genes, rather than just the top 10 expressed per cluster.
Note: this is for v3.1 of Seurat; the latest version of Seurat at the time of asking.


text:
length of 'dimnames' [1686] must match that of 'dims' [3]
Please if anyone has experience with the use of the BSEQ-SC package for the deconvolution of bulk RNA sequencing data with single cell RNA sequencing data I will be very grateful for your suggestion.
I wish to deconvolve my bulk RNA seq data with single cell RNA seq data from the same tissue. I run in to the following error when running the code line below and I do understand the error but not how I can overcome it. Can any one please provide any useful hints?
csfit <- bseqsc_csdiff(ilc_bulk[genes, ] ~ treatment | fibroblast + macrophage,
                       verbose = 2, nperms = 100, .rng = 12345)

Then comes the error:
   Groups: a_nacl1=4L | a_nacl2=4L | a_nacl3=4L | a_nacl4=4L | a_nacl5=4L | NA0L
Cell type(s): 'fibroblast', 'macrophage' (2 total)
Fitting mode: auto
Data (filtered): 1684 features x 20 samples
Model has factor effect with more than 2 levels: fitting lm interaction model
Fitting model with nonnegative effects
Model with more than 2 groups: switching to version 2
 Fitting linear interaction model ...  OK
 Computing FDR using 100 permutations ... 101/100
  Alternative 'two.sided' ...   OK
  Alternative 'greater' ...   OK
  Alternative 'less' ...   OK
 OK
Timing:
   user  system elapsed
  2.136   0.153   2.267
Error in array(NA_real_, dim = c(nrow(b), ncol(cc), length(lev)), dimnames = c(rownames(b),  :
  length of 'dimnames' [1687] must match that of 'dims' [3]

> sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=de_DE.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=de_DE.UTF-8       LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] edgeR_3.26.6                limma_3.40.6                forcats_0.4.0               stringr_1.4.0
 [5] purrr_0.3.2                 readr_1.3.1                 tidyr_0.8.3                 tibble_2.1.3
 [9] tidyverse_1.2.1             SoupX_0.3.1                 ggplot2_3.2.0               dplyr_0.8.3
[13] DropletUtils_1.4.3          SingleCellExperiment_1.6.0  SummarizedExperiment_1.14.0 DelayedArray_0.10.0
[17] BiocParallel_1.18.0         matrixStats_0.54.0          GenomicRanges_1.36.0        GenomeInfoDb_1.20.0
[21] RColorBrewer_1.1-2          xbioc_0.1.17                AnnotationDbi_1.46.0        IRanges_2.18.1
[25] S4Vectors_0.22.0            BisqueRNA_1.0               Seurat_3.0.2                preprocessCore_1.46.0
[29] e1071_1.7-2                 bseqsc_1.0                  csSAM_1.4                   Rcpp_1.0.2
[33] openxlsx_4.1.0.1            Biobase_2.44.0              BiocGenerics_0.30.0

loaded via a namespace (and not attached):
  [1] reticulate_1.13        R.utils_2.9.0          tidyselect_0.2.5       RSQLite_2.1.2          htmlwidgets_1.3
  [6] grid_3.6.1             Rtsne_0.15             devtools_2.1.0         munsell_0.5.0          codetools_0.2-16
 [11] ica_1.0-2              future_1.14.0          withr_2.1.2            colorspace_1.4-1       rstudioapi_0.10
 [16] ROCR_1.0-7             gbRd_0.4-11            listenv_0.7.0          NMF_0.22               Rdpack_0.11-0
 [21] labeling_0.3           GenomeInfoDbData_1.2.1 bit64_0.9-7            rhdf5_2.28.0           rprojroot_1.3-2
 [26] vctrs_0.2.0            generics_0.0.2         R6_2.4.0               doParallel_1.0.14      rsvd_1.0.2
 [31] locfit_1.5-9.1         bitops_1.0-6           assertthat_0.2.1       SDMTools_1.1-221.1     scales_1.0.0
 [36] gtable_0.3.0           npsurv_0.4-0           globals_0.12.4         processx_3.4.1         rlang_0.4.0
 [41] zeallot_0.1.0          splines_3.6.1          lazyeval_0.2.2         broom_0.5.2            modelr_0.1.4
 [46] BiocManager_1.30.4     yaml_2.2.0             reshape2_1.4.3         backports_1.1.4        tools_3.6.1
 [51] usethis_1.5.1          gridBase_0.4-7         gplots_3.0.1.1         sessioninfo_1.1.1      ggridges_0.5.1
 [56] plyr_1.8.4             zlibbioc_1.30.0        RCurl_1.95-4.12        ps_1.3.0               prettyunits_1.0.2
 [61] pbapply_1.4-1          viridis_0.5.1          cowplot_1.0.0          zoo_1.8-6              haven_2.1.1
 [66] ggrepel_0.8.1          cluster_2.1.0          fs_1.3.1               magrittr_1.5           data.table_1.12.2
 [71] lmtest_0.9-37          RANN_2.6.1             fitdistrplus_1.0-14    pkgload_1.0.2          hms_0.5.0
 [76] lsei_1.2-0             xtable_1.8-4           readxl_1.3.1           gridExtra_2.3          testthat_2.2.1
 [81] compiler_3.6.1         KernSmooth_2.23-15     crayon_1.3.4           R.oo_1.22.0            htmltools_0.3.6
 [86] Formula_1.2-3          lubridate_1.7.4        DBI_1.0.0              MASS_7.3-51.4          Matrix_1.2-17
 [91] cli_1.1.0              R.methodsS3_1.7.1      gdata_2.18.0           metap_1.1              igraph_1.2.4.1
 [96] pkgconfig_2.0.2        registry_0.5-1         plotly_4.9.0           xml2_1.2.1             foreach_1.4.7
[101] dqrng_0.2.1            rngtools_1.4           pkgmaker_0.28          XVector_0.24.0         rvest_0.3.4
[106] bibtex_0.4.2           callr_3.3.1            digest_0.6.20          sctransform_0.2.0      tsne_0.1-3
[111] cellranger_1.1.0       dendextend_1.12.0      curl_4.0               gtools_3.8.1           nlme_3.1-141
[116] jsonlite_1.6           Rhdf5lib_1.6.0         desc_1.2.0             viridisLite_0.3.0      pillar_1.4.2
[121] lattice_0.20-38        httr_1.4.0             pkgbuild_1.0.3         survival_2.44-1.1      glue_1.3.1
[126] remotes_2.1.0          zip_2.0.3              png_0.1-7              iterators_1.0.12       bit_1.1-14
[131] class_7.3-15           stringi_1.4.3          HDF5Array_1.12.1       blob_1.2.0             caTools_1.17.1.2
[136] memoise_1.1.0          irlba_2.3.3            future.apply_1.3.0     ape_5.3

Answer from @stupidwolf, @benn, @llrs, converted from comments:
The error complains about the 1686 rownames, while you have 1684 features. I think something is wrong with the naming of your rows.
This may be a bug that you should report to the package maintainer. I cannot really tell without more information about your eset. Looking at this line "Cell type(s): 'fibroblast' (1 total)", is it possible you only have one group under fibroblast? It will not make sense to run a random effect for only 1 group.
[please edit to improve this answer, if possible]


text:
Biplot scores from canonical correspondence analysis
I'm using the R package vegan to perform canonical correspondence analysis (CCA). As input we have two matrices, one being (sites)x(species) and the other being (sites)x(conditions).
Sample data (and source of plot) are here.
Species loadings are easily accessed with summary(cca_model)$species. What I'm trying to find is the loadings for the explanatory variables, the conditions. The only summary I can find is biplot scores. Looking through the documentation for vegan I can't find any description of how they are computed. Can I sum them across CCA components to get an idea of how much they influence the data?
This is a biplot of two CCA components. The scores are used as coordinates for the arrows.

So two questions:

What are biplot scores in the context of CCA?

and

Can biplot scores be used to determine how much of an effect conditions have on the response variables?

Canonical correlation analysis (CCA) not correspondence canonical analysis. Biplots are the same as in PCA: the representation in a new space of the samples and the relative position of the variables on the new space.
However, CCA seeks the values of one block that maximize the correlation with the other block. Here you don't have an independent variable that you maximize from the variables you introduced on the CCA. So the interpretation is much trickier as there isn't a direction: WaterC and FallenTwigs show the opposite effect than BarSand and CoverMoss, but there is a group of variables that explain what is the effect of water on the plants(?)

CCA refers to canonical correspondence analysis, at least in vegan. I've updated the question


text:
How to run function "convolve" in parallel in R?
I have a vector of 200 million numbers. It has been running for over 24 hours and still not done.
x <- data   # the target sequence
k <- exp(-(-win:win)^2/(2*sig^2))/(sig*sqrt(2*pi))  # gaussian kernel
ptm <- proc.time()
y <- convolve(x,k,type="open")
proc.time()-ptm

How to parallelize the convolve function?
To parallelize, I tried splitting up the vector into pieces:
> x <- data[1:1e7]
> ptm <- proc.time()
> y <- convolve(x,k,type="open") # smoothed sequence
> proc.time()-ptm
   user  system elapsed
 40.743   0.206  40.940
> x <- data[1e7+1:2e7]
> ptm <- proc.time()
> y <- convolve(x,k,type="open") # smoothed sequence
> proc.time()-ptm
   user  system elapsed
320.193   0.708 320.871

However, splitting it up and then adding back the convolutions will be tricky.
Does anyone have other ideas?
This data is based on biological data. We want to cluster values in the vector based on nearby vector indices.
Reproducible example with 20,000,000 random numbers should finish in about 5.5-6 minutes. But if you try running 200,000,000 locally or on HPC, it will take a long time (my HPC has limit 99 hours and it timed out at 99 hours unfinished).
win <- 100        # window length of gaussian kernel (units acc to orig sequence)
sig <- 5        # std deviation of gaussian
k <- exp(-(-win:win)^2/(2*sig^2))/(sig*sqrt(2*pi))  # gaussian kernel
l <- 20000000
x <- round(runif(l, 0.0, 1.0),1)   # the target sequence with peaks
t <- seq(0,1-1/length(x),1/length(x))
ptm <- proc.time()
y <- convolve(x,k,type="open")[(win+1):(win+length(x))] # smoothed sequence
proc.time()-ptm

Is the function vectorized? Did you try the same function with less numbers? If you need to make a cluster, why do you need to make the convolution? Perhaps you will get the attention of those that could help you better on stackoverflow.com.

Hi IIrs: The convolve function is a function in R. The input data (x <- data) is a vector. We need to group based on nearby indices in the vector, so we cannot use other clustering methods that were built, for example, clustering gene expression. I can re-post in stackoverflow, but I have heard re-posting the same question is not good (and I am new to the stackoverflow community). I can try re-posting though. Thank you.

What kind of vector is it? A list or a numeric one? Could you please post a minimal reproducible example (with faked data if needed)?

Hi IIrs, It is a numeric vector. I've troubleshooted and number of decimal places does not affect time to finish. I have edited the post above with a reproducible example with 20,000,000 random numbers. The time to finish with random numbers matches the time to finish with subsetted real data.

I don't know what is supposed to do but maybe the problem is that there isn't any continuous periodicity to use fft on it??


text:
How to apply RDP, Greengenes and other special taxonomies in Krona?
I used Kraken 2 to classify my 16S metagenomic data using both RDP and Greengenes database. As these are special databases the taxonomic ids assigned  do not match with the NCBI taxonomic ids. I found Krona uses NCBI taxonomy as its default. Krona github page has the following instruction regarding RDP database:

Taxonomic classifications of the Ribosomal Database Project 16S
  Classifier can be imported using ktImportRDP. To download
  classifications:
Run the classification tool.
Click "show assignment detail" for the root of the hierarchy.
Click "download allrank result" or "download fixrank result".

Problem is the online RDP classifier (http://rdp.cme.msu.edu/classifier/classifier.jsp) doesn't work and as far as I understood the above mentioned help applies only for web version of the tool.
And for Greengenes taxonomy id file found here , it is in different format than the one Krona uses. Krona taxonomy file looks like this:
1   0   1   no rank root
2   2   131567  superkingdom    Bacteria
6   7   335928  genus   Azorhizobium
7   8   6   species Azorhizobium caulinodans
9   8   32199   species Buchnera aphidicola
10  7   1706371 genus   Cellvibrio
11  9   1707    species Cellulomonas gilvus
13  7   203488  genus   Dictyoglomus
14  8   13  species Dictyoglomus thermophilum
16  7   32011   genus   Methylophilus
17  8   16  species Methylophilus methylotrophus
18  8   213421  genus   Pelobacter
19  9   18  species Pelobacter carbinolicus
20  7   76892   genus   Phenylobacterium
21  8   20  species Phenylobacterium immobile
....

And Greengenes taxonomy file mentioned above looks like following:
228054  k__Bacteria; p__Cyanobacteria; c__Synechococcophycideae; o__Synechococcales; f__Synechococcaceae; g__Synechococcus; s__
844608  k__Bacteria; p__Cyanobacteria; c__Synechococcophycideae; o__Synechococcales; f__Synechococcaceae; g__Synechococcus; s__
178780  k__Bacteria; p__Cyanobacteria; c__Synechococcophycideae; o__Synechococcales; f__Synechococcaceae; g__Synechococcus; s__
198479  k__Bacteria; p__Cyanobacteria; c__Synechococcophycideae; o__Synechococcales; f__Synechococcaceae; g__Synechococcus; s__
187280  k__Bacteria; p__Cyanobacteria; c__Synechococcophycideae; o__Synechococcales; f__Synechococcaceae; g__Synechococcus; s__
179180  k__Bacteria; p__Cyanobacteria; c__Synechococcophycideae; o__Synechococcales; f__Synechococcaceae; g__Synechococcus; s__
175058  k__Bacteria; p__Cyanobacteria; c__Synechococcophycideae; o__Synechococcales; f__Synechococcaceae; g__Synechococcus; s__
176884  k__Bacteria; p__Cyanobacteria; c__Synechococcophycideae; o__Synechococcales; f__Synechococcaceae; g__Synechococcus; s__
228057  k__Bacteria; p__Proteobacteria; c__Alphaproteobacteria; o__Rickettsiales; f__Pelagibacteraceae; g__; s__

I am stuck here. Where can I find taxonomy databases that use the format Krona uses? Is there any tool I can use to convert between these formats?

Same problem here!


text:
Calculating Unifrac Distance and missing data
I am trying to calculate unifrac distance for my OTU table, which has about 3000 taxa, and 300 samples. My phylogenetic tree only considers the most important 500 taxa. My questions are:

Do I ignore the 2500 taxa that were not considered in the phylogenetic tree?
Is their a package in R that deals with this kind of missing data?

Welcome, your request is a bit unusual. Could you provide details about the way you have built the tree please? I don't know of any tree package that would remove OTUs unless you are collapsing clades.

Sorry I am also a bit confusd by this "which has about 3000 taxa, 300 samples" ... that seems to imply 10 taxa per sample .. is this metagenomics?

To calculate the unifrac distance you need the phylogenetic tree. If you only have it for 500 taxa then you can only calculate it for those taxa.
It is impossible to impute the phylogenetic tree without knowing the sequence on which the taxa is based. But if you know it then you can calculate the phylogenetic tree for the 3k taxa you have.


text:
Duplicated genes and genome fragmentation
This is a follow up question to this post: How can I find out if my gene of interest is duplicated?
Briefly, to identify homologous genes, I performed tBLASTn analyses against my genome of interest using various sequences of a closely related species as queries for the BLAST analyses. Using this method, I found two (and in some cases up to five copies) of my different genes of interest in the target genome. These genes are involved in different organelle functions (i.e., they are not surface antigens) and I want to find out whether these sequences represent duplication events or not. These analyses are being performed on protozoans parasites that belong to the same genus. The query genome is properly annotated, but the target genome is not (it is a whole genome shotgun project that has been assembled at the scaffold level).
Following the advice of @Michael G., I re-analyzed my sequences and blasted them against the entirety of Genbank. If there are two sequences for the same protein (for e.g. protein 3-1 and protein 3-2) and I compare them among themselves they are very similar (over 70% identity at the protein level). When I blast each sequence against the whole protein RefSeq database, the first hit is my query and the rest of the hits correspond to the same protein in other organisms. I believe that these sequences represent duplication events and it is possible that they are paralogues, but more evidence is needed to classify them as such.
However, when there are more than two sequences for the same protein, I found something different. For e.g., say I have three sequences for the same protein: sequence 5-1, sequence 5-2 and sequence 5-3. Sequences 5-2 and 5-3 are very similar (78% identity at the protein level) but when I compare sequences 5-1 and 5-2 I find that they are overlap in the middle, i.e., they are part of the same sequence. Thus, it appears that I don't really have three sequences for this protein but two. Blastp analyses against the whole RefSeq database give the expected result (first hit is my query sequence, the rest of the hits are same sequence in other organisms). So, once again I have what appears to be a duplication event, but in this case one of the sequences appears to be split in two -- could this be because the target genome is not completely assembled? Why could this be happening?
Thank you very much for your help.
-Leah
Edit: The new results are posted in this thread: Duplicated genes and genome fragmentation - Part 2

Thanks @leah, let me think. The "duplication event" will depend on the sample of other protist genomes available.

Thanks for your reply @MichaelG. I am sorry, but I don't quite understand... What do you mean by "the sample of other protist genomes available"?

The tBlastn looks good and looks like some sort of parology, although the hits you are getting need to be checked to make sure they are not psuedo-genes. The shot-gun genome however appears poorly annotated because it had misannotated 5-1 and 5-2. What I don't understand is why you haven't performed a blast back over the original highly annotated genome? You must have done this? There are a few avenues to go down here, what I'd like to know which group of protists are you working with? I don't need to know the species, but I do need to know the family please.

Can you two please convert this to a Q&A instead of a discussion thread? Leah, please take all the information and make it into one, comprehensive question. @MichaelG. please post one answer that answers that comprehensive question. I don't really know what else to do here without deleting everything. Please don't use this site like a forum.

In the mean time, I have closed this as too broad since it is very much too broad for this site in its current form.

Hi @terdon, the issue is that it is different questions each time. The current question regards long branches and long branch attraction (LBA) in a phylogeny. This could be posted as a separate question. What is a Q&A thread?

@MichaelG. the idea here is to have one question and then everything else would be an answer to that question. If the question requires discussion like this, then it is too broad for the SE model and should be closed, and the question posted on SeqAnswers or BioStars instead since those to discussion. SE tries to have a high signal to noise ratio, which means that questions and answers need to be clearly separate and easy to find. When you have multiple back and forths like this, everything is confusing.

[cont] For instance, because of the voting system, there is no guarantee that the order of the answers will remain the same, so it's impossible to figure out who's answering whom and which answer is posted in answer to which question. The whole thing becomes a mess and ends up helping nobody except perhaps the people involved in the discussion. But we want to set up a resource that will help everyone, including the user who stumbles upon this question in 5 years time where the order of the various answers will be totally different. :/

In summary, apologies for the rolling nature of the question and answer and I will know in future. It is an interesting problem. If the OP targets specific questions, thus 1) long branch attraction (which is a big topic), 2) state-of-the-art tree amino acid tree construction (which is applicable to other studies), it should be compliant with forum rules?

The strategy from here is to determine duplication/parology through phylogeny, notably amino acid phylogeny using statistical support.
The basic result is whether the genes under investigation are:

A. monophyletic to the given species under investigation
B. Paraphyletic against protist species which are not under investigation,

Point A is compatible with paraology/duplication within the species. Point B infers the duplication is ancient.
The skill/difficulty is building a comprehensive alignment.
Before we start,

Parology is usually underestimated within a genome because repetitive sequences are often collapsed in the assembly process. Thus a de novo assembly around the gene of interest in both annotated and shot-gun assemblies might yield better results, ie. more duplicates. You have already seen what appear to be N-terminal and C-terminal ends of a related gene fragmented in the assembly process.
The tBlastn results looked quite cool and need to be included.

The core approach is to blastp your query sequence against the full NBCI database. Blast has an option to automatically download all hits (or all selected hits), e.g. in fasta format. This will provide a multi-species data set to understand the evolution of the protein of interest. This feature of blast is very cool and saves a lot of work. You want to combine the output of point 1 and point 2 in the fasta blastp list and align everything via ClustalO or Muscle (any alignment program). Then manually edit the AA alignment via an alignment editor (I use Seqotron [Mac]) and thereon produce a very basic tree using neighbor-joining just to check everything looks ok. The NCBI blast trick might need a better explanation so let me know if you have trouble figuring this.
If you post the tree as a PHYLOGRAM (I don't need to see the sequence names - just the shape and branch lengths) we can then see whether it looks reasonable before commencing a heavy duty algorhitm to identify support for the competeting hypotheses.
It is not ideal doing this via a message server, but its doable.

Thank you very much for the detailed explanation! I understand the basics of what you described, so I will give it a try. I hope that I will be able to download all the hits from blast... I will let you know.  Again, many thanks for all your help.

The tBlastn looks good and looks like some sort of parology, although the hits you are getting need to be checked to make sure they are not psuedo-genes. The shot-gun genome however appears poorly annotated because it had misannotated 5-1 and 5-2. What I don't understand is why you haven't performed a blast back over the original highly annotated genome? You must have done this?
There are a few avenues to go down here, what I'd like to know which group of protists are you working with? I don't need to know the species, but I do need to know the family please.

Many thanks for your response!

What I don't understand is why you haven't performed a blast back over the original highly annotated genome? You must have done this?

Yes, I have performed a blast back over the original highly annotated genome and, as expected, the first hit is the sequence that I used as query in the first blast analysis.

There are a few avenues to go down here, what I'd like to know which group of protists are you working with? I don't need to know the species, but I do need to know the family please.

I am working with Amoebozoans -- specifically, with members of the group Discosea - Centramoebida.
Once again, thank you very much for your help.

That is fine. From here it's difficult stuff and you need to build trees, providing advice over a forum like this is tricky but we'll have a go.

So this is what I did to produce the phylogenetic tree:

I performed a BLASTp using the query sequence (lets call it sequence 5 from species A, which is thoroughly annotated) against full NCBI protein database (ref-seq database)
I downloaded all the hits in fasta format (very cool feature indeed!)
Here I was a bit confused, so what I did was to take the three hits I got in my target genome (sequences 5-1, 5-2 and 5-3 from species B, which is poorly annotated) and added them to the fasta blastp list
I took the fasta blastp list and loaded it into Seqotron.
I aligned everything using the Compute - Alignment - MUSCLE option -- I didn't do any manual AA alignment (not too sure about this part)
I produced the tree using Compute - Phylogenetic Tree (default parameters: distance matrix: K83, resampling: none).
This is the tree I got. Sequences from the target genome (species B) are shown with a red rectangle:

Please let me know if I missed something or if I made a mistake somewhere. Many thanks!
-Leah

Cool @leah, this looks good. There is an issue with 2 of the taxa, the ones in red, otherwise if this was correctly rooted it looks fine. I'll post further things to check much later today

The first thing is the Sequence 5-2 and 5-3 are on very long branches. This is often indicative of sequence problems. You mentioned 5-2 and 5-3 are possibly two fragments of the same gene. My advice is to do a psuedo-annotation and reassess the phylogeny to see whether that super long branch length falls back in line
The branches close to 5-1 are of interest and their genomes should be further investigated.
You need to find the root of the phylogeny using external criteria. I can't do this because I've no idea was the species are.

Long branches are bad news in phylogeny, particularly if artefactual because of a phenomena known as 'long branch attraction'. Don't use your tree to understand duplication until this is resolved etc ... because the results can be distorted by LBA. Thus you must either delete them from the tree or perform a psuedo-annotation and redo the tree.
I was reluctant to forward the idea of phylogeny, because it is not just a bit of coding but an entire dialogue. Anyway by correcting the presumed misannotation you have sufficient results to move the investigation forward. You cant delineate your hypotheses until the root is correctly determined and then you set this as your outgroup under a more robust tree-building algorhithm.
I think we are out of space on this thread (you're allowed 1 post, I'm allowed 2), but clearly things have moved forward. Again, please note once the long-branch is gone you can end up with a different tree.

Thanks for the detailed explanation. I will annotate sequences 5-2 and 5-3 and run the analysis again to see if the branch length is shortened. I will also look into the branches close to sequence 5-1 and see what I can make of them. I will post those results in a new thread. I will also look into finding the root of the phylogeny, but I am not so sure how to go about it... Which external criteria could I use? Thanks again for everything!

The new results are posted in this thread: https://bioinformatics.stackexchange.com/questions/10365/duplicated-genes-and-genome-fragmentation-part-2


text:
Calculate the percentage of the genome that is represented by some sequences
I have construct a repeat libraries using de-novo identification methods. All  results from different tools has been combined into one fasta file.
How can I get the percentage of how much these sequences are found in the original genome? That has been used earlier for repeat detections.
For example these identified candidates account for 30% from the total genome.
I used repeatmasker with a flag -lib combinLibrary.fasta mygenome.fasta but its very slow now almost one week and still far to finish.

Do you want to know how much of the genome is in your sequences in a fasta file? What have you tried to do? Any perl/python/R/bash or other tools you have used or are familiar with?

i combined all the library then i used repeatmasker with a flag - lib combinLibrary.fasta mygenome.fasta . but its very slow now almost one week and still far to finish, im not sure if its right but im waiting the result while waiting i thought of asking for help

Great, thanks for clarifying this. I'll edit your question to include this information. You can also [edit] the question to add and improve the question whenever you want

i thought its easy task for the expert by using any kind of scripting language , my run still far from half time. any help please

Sorry I am not an expert on genome coverage. Perhaps if you link around in social networks you get more views and hopefully an answer.

how many repeats do you have?  have you tried repeatmasker with a smaller set to validate your method?  This would allow you to estimate time to completion for the full dataset.  Have you dereplicated your repeats?  That may be necessary/useful depending on how collected them, since you didn't give any details.

@EdwardKirton I collected the repeats using a structure-based method, each class of TEs detected with a different tool, after validation and discarding false positive, all the results merged in one fasta file, now my question is to know the coverage of this repeat occured in the whole genome.


text:
Python module for fetching NCBI id for a list of species
I have a list of scientific names of species. Is there a python module that can fetch NCBI taxonomy IDs?

Biopython has a collection of NCBI Entrez query modules. Please check this StackOverflow thread.
https://stackoverflow.com/questions/16504238/attempting-to-obtain-taxonomic-information-from-biopython

Biopython has an efetch frontend.
from Bio import Entrez
Entrez.email = "mg@example.edu"
handle = Entrez.efetch(db="nucleotide", id="KR052012.1", rettype="gb", retmode="text")
print(handle.read())

OUT
LOCUS       KR052012               10714 bp    RNA     linear   VRL 09-MAY-2015
DEFINITION  Dengue virus strain Hb33/CHN/2014, complete genome.
ACCESSION   KR052012
VERSION     KR052012.1
KEYWORDS
.... etc
You will likely need to perform an EFetch search and then extract out the accession numbers for import (as done above)

This isn't a Python library, but it's easily invoked from Python code, easy to install, and has a lot of useful features.

I just discovered the amazing taxonkit library this week. In particular, the taxonkit name2taxid command seems to be what you're looking for.
$ cat NAMES
Polistes dominula
Sinorhizobium meliloti
Gossypium hirsutum
$ cat NAMES | taxonkit name2taxid
Polistes dominula       743375
Sinorhizobium meliloti  382
Gossypium hirsutum      3635

Easily installed using conda install taxonkit.

If you want to perform a species search, here's the raw efetch code. I use this within an alias and the complex quotes are needed for complicated reasons. The reg-ex may need to be adapted to your search. $virus here is your species of interest, e.g. set virus='some species' (csh or tcsh), `virus='some species' (bash). You'll just need to translate this to Biopython.
alias esummary-acc \
'esearch -db protein -query "$virus" |\
efetch -format docsum |\
xtract -pattern DocumentSummary -element Title,AccessionVersion |\
grep "$virus" |\
perl -p -e '"'"'s/^.*\[.*\]\s+(\w+)/$1/'"'"''