title: How I can reproduce this heat map

text:
# Question
Title: How I can reproduce this heat map
I have single cell RNA-seq on 8 time points. I have clustered each time point to 2-3 clusters by seurat. I also have a list of differentially expressed genes between all my 1562 cells in 8 time points. I know how to plot this list of genes in heat map on each cluster of cells or a big heat map of my 1562 cells, but I have also sub-clusters for each cluster (time point).

Could you please help me to figure out how can I have a heat map with sub-clusters for each time point on which I am showing a list of differentially expressed genes?

Something like the picture below where each day has been clustered to 2-3 sub-population then they plotted common differentially expressed genes between day points

**Edit**

Lets say 500 genes commonly are being differentially expressed in 8 time point. time point 1 has 3 sub-clusters of cells and time points 2 has 2 sub-clusters of cells. time point 3 and 4 both have 2 sub-clusters of cells. I want to heat map these 500 genes in cells in time point 1 (3 sub-clusters), time point 2 (2 sub-clusters), time point 23 (2 sub-clusters) and time point 4 (2 sub-clusters) in one unified heat map on which cells in each subclusteres (a,b and c) are in column and genes in rows but I don't have any idea how to do that. column annotation function in ComplexHeatmap likely works in this regard But I don't know how to define cells inside each sub cluster in different colours

I can heat map of 500 genes only on time points easily by seurat but without underlying sub clusters

```
> DoHeatmap(object = SubsetData(object = object, max.cells.per.ident =
> 100), genes.use = 50 genes, slim.col.label = TRUE, group.label.rot =
> TRUE)

```

Figure 2 in this article.

# Answer
> 3 votes
In case it is still helpful as the post is rather old, the code below would generate a heatmap with annotations thanks to the `ComplexHeatmap` package.

But before that I would like to stress out the importance of example data, I bet a lot of experts could not help you just because of they did not have single cell data at their disposal.

I am using some scRNA-seq data stored within a `Seurat` object. The expression data come from the `@scale.data` slot and the annotations are from the `@meta.data` slot. I am using a few general cell type markers as genes and 3 columns from the metadata slot (cell type, cell subtype, sample region). When using the `ComplexHeatmap` package, the trick is to have the same order in the annotation and actual data tables and this is how you should control what (in your case single cells) is plotted in what order.

```
library(Seurat)
library(dplyr)
library(ComplexHeatmap)
library(circlize) # for the colorRamp2() function

seurat_object <- readRDS(...)

# Annotations

meta_data <- seurat_object@meta.data

# order of annotations/colors are defined here
ordered_meta_data <- meta_data[order(meta_data$cell_type), ]

# OPTIONAL: YOU CAN PICK COLORS FOR EACH LEVEL OF ANNOTATION
# HERE I PROVIDE COLORS FOR ONLY TWO OF THE THREE LEVELS, COLORS FOR THE
# REMAINING LEVEL IS TAKEN CARE OF BY THE PACKAGE
annotation_colors <- list("cell_type"= c("Epithelial" = "red",
                                         "Immune" = "orange",
                                         "Other" = "gray",
                                         "Stroma" = "green"),
                          "sample_region" = c("B" = "blue",
                                              "C" = "yellow",
                                              "N" = "maroon",
                                              "T" = "darkblue"))
ha = HeatmapAnnotation(df = ordered_meta_data,
                       show_annotation_name = TRUE,
                       col = annotation_colors)

# Expression data

genes_to_use <- c("EPCAM",
                  "CD3D",
                  "NCR1",
                  "CD79A",
                  "LYZ",
                  "COL1A1",
                  "CLDN5",
                  "S100B",
                  "KIT")

seurat_object <- ScaleData(seurat_object,
                           genes.use = genes_to_use)

my_data <- seurat_object@scale.data

# COLUMN ORDER OF THE EXPRESSION DATA SHOULD MATCH THE ROW ORDER OF THE
# ANNOTATION TABLE
my_data <- my_data[, rownames(ordered_meta_data)]

# Heatmap

col_fun = colorRamp2(c(-2, 0, 2), c("blue", "white", "red"))
png(filename = "example_heatmap.png",
    width = 1000,
    height = 1000)
Heatmap(
  my_data,
  col = col_fun,
  cluster_rows = FALSE,
  cluster_columns = FALSE,
  column_order = NULL,
  show_row_dend = FALSE,
  show_column_dend = FALSE,
  show_row_names = TRUE,
  show_column_names = FALSE,
  use_raster = TRUE,
  raster_device = c("png"),
  bottom_annotation = NULL,
  top_annotation = ha
)
dev.off()

```

---
Tags: r, scrnaseq, seurat, heatmap, figure-reproduction
---


# Question
Title: 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?

# Answer
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.
> 1 votes

---
Tags: api, galaxy
---


# Question
Title: 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...

# Answer
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!
> 0 votes

---
Tags: r, rna-seq, filtering
---


# Question
Title: 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;
    |                      ^~~

```

# Answer
> 1 votes
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.

---
Tags: c++
---


title: Create a Grouped Barchart

text:
# Question
Title: 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 ?

# Answer
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")

```
> 4 votes

---
Tags: r, statistics, ggplot2
---


# Question
Title: 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?

# Answer
> 0 votes
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.

# Answer
> 0 votes
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

1. Go to Blast here, https://blast.ncbi.nlm.nih.gov/Blast.cgi
2. Enter your sequence into the box (it doesn't accept PDB codes alone)
3. 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)
4. Hit "Blast"
5. Once the search is complete at the top of the page are the hyperlinks: "Other reports: Search Summary \[Taxonomy reports\] \[Distance tree of results\]"
6. Click on "Distance tree of results"
7. 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)
8. 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

---
Tags: sequence-alignment, phylogenetics, phylogeny, blastn
---


# Question
Title: 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?

# Answer
> 2 votes
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.

---
Tags: ucsc
---


# Question
Title: How to create Seurat object while RNA expression and ADT combined into one matrix
I got an input matrix which RNA expression and ADT capture are combined into one file. I loaded the file into Seurat successfully, however, when I tried to create Seurat object, it threw out an error saying

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

```

How can I create Seurat object when RNA expression and ADT data are combine together?

# Answer
> 1 votes
Given the fact that you have RNA and ADT data (probably on the same set of cells), some colnames (which represent individual cells) would be duplicated. Normally duplicated colnames are tolerated with matrices and sparse matrices, however, Seurat apparently does not do so, below is from the `?CreateAssayObject`:

```
Details
Non-unique cell or feature names are not allowed. Please make unique before calling this function.

```

Possible solution to your problem: Seuart has a dedicated vignette for working with multimodal data and as you would see you will need to initiate your Seurat object with one matrix per assay: RNA and ADT. All you need to do is split your matrix into RNA and ADT, create your Seurat object with RNA data and then add the ADT data with:

`seurat_obj_with_rna_only[["ADT"]] <- CreateAssayObject(counts = your_adt_matrix)`

For efficiency, Seurat uses sparse matrices so don't forget to convert your data matrices to sparse.

# Answer
> 0 votes
I don't know if you solved the problem. Howerver,I think it's caused by multiple data types.

For output from CellRanger \>= 3.0 with multiple data types

```
data <- Read10X(data.dir = data_dir)

```

This step creates multiple matrices, depending on how many types of data there are.

```
seurat_object = CreateSeuratObject(counts = data$`Gene Expression`)

```

So (counts = data$`Gene Expression`),You should also be able to choose the type of research you want to do.

I hope my answer will help some people.

---
Tags: seurat, single-cell
---


# Question
Title: 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.

# Answer
> 1 votes
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.

---
Tags: database, sequence-annotation, mutations, non-coding
---


# Question
Title: 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.

# Answer
> 1 votes
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.

---
Tags: gene-expression, drugs
---


# Question
Title: How to compute LD for pairs of variants with Plink
I have a set of SNPs (single nucleotide polymorphisms) { S1, S2, ..., SN } from approximately 200 humans. I wish to determine the linkage disequilibrium (LD) for a defined subset of these SNPs. By this I mean e.g. I would like to compute e.g. LD between S1 and S2, S3, S4, but not the LD between S2 and S3.

Could this calculation be performed in PLINK?

# Answer
Yes, you can do this with plink.

You can use the flag `--ld-snp <variant_id>` to calculate LD between your SNP of interest and all other SNPs.
If you have multiple SNPs of interest, and want each of those compared with every SNP in your dataset, use `--ld-snps <snp1-snp100,snp102,snp104>`
A final option to calculate LD between a specific pair of SNPs is `--ld <variant_id1> <variant_id2>`

This the relevant plink manual page
> 2 votes

# Answer
Anyway, it would be fairly easy via the filtering options here

The `--exclude-snp` is what you are wanting.

---

I have worked in LD and subsampling a genome is open to forcing a result. The alternative approach is the subsample is the population, or examine disequilibrium with respect to the physical chromosome. A local shift in disequilibrium or heterozygosity (sliding window along a chromosome) is an indication of selective pressure. This approach is frequently used in eukaryotic parasites to identify drug resistance. I don't have experience of human genetics but there are certainly enough here who do. However, I think you need to define the biological rationale for the subsample. Just having loads of data is easily resolved because any algorithm can be written to handle automation.
> 0 votes

---
Tags: vcf, snp, variation, plink
---


title: Adding an attribute to GFF3 file

text:
# Question
Title: 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

# Answer
> 2 votes
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)

```

---
Tags: python, sequence-annotation, gff3, gffutils
---


# Question
Title: 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

```

# Answer
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
> 2 votes

# Answer
You might want to check this line out:

`for gene in db.features_of_type('gene', order_by='start'):`

You are specifically filtering for `gene`s hence other types of annotations are filtered out.
> 1 votes

# Answer
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

```
> 1 votes

---
Tags: python, sequence-annotation, gff3, gffutils
---


# Question
Title: Getting a stretch of genomic ranges from a dataframe/granges object based on metadata column
I have a "test" `data.frame` object in R, which is basically a small subset of a 66000 row dataframe, which looks as follows:

```
chr   start    end      PC1     comp
chr1  3360001  3400000  -0.009  B
chr1  3400001  3440000  -0.004  B
chr1  3440001  3480000  -0.001  B
chr1  3480001  3520000  0.003  A
chr1  3520001  3560000  0.002  A
chr1  3560001  3600000  -0.001  B
chr1  3600001  3640000  0.002  A
chr1  3640001  3680000  0.004  A
chr1  3720001  3760000  0.003  A
chr1  3760001  3800000  0.003  A
chr1  3800001  3840000  0.007  A
chr1  3840001  3880000  0.006  A
chr1  3880001  3920000  -0.004  B
chr1  3920001  3960000  -0.017  B

```

As you can see, the interval between each start and end position is `40kb`. The column `comp` can have only two values, `A` or `B`. Going from `row 1` to the last row of my dataframe and based on the values of the column `comp`, I want to get the whole stretch/block of the genomic region, once the compartment(character in column `comp`) is changed, and perform a `mean` of the `PC1` values in that stretch.

So in this example, I basically want to shrink my test dataframe as follows:

```
chr1  3360001  3480000  mean(c(-0.009,-0.004,-0.001))  B
chr1  3480001  3560000  mean(c(0.003,0.002))  A
chr1  3560001  3600000  mean(c(-0.001))  B  <---- #notice this is a 'singleton'
chr1  3600001  3880000  mean(c(0.002,0.004,0.003,0.003,0.007,0.006))  A
chr1  3880001  3960000  mean(c(-0.004,-0.017))  B

```

I tried with a simple code, but it only runs until the first compartment value change(so row 4 in this test), and I am unable to figure out how to make it run for the whole dataset of 66000 rows. Also the use of multiple `for` loops, which is what I tried, may not be the prettiest way to do this.

I guess there exists a `Bioconductor` package which performs this task very easily. If anyone can point me to that direction, it'd be great.

# Answer
The GenomicRanges package in bioconductor will do this. It lets you do 'set' operations on genomic ranges, such as concatenating them, finding overlaps, etc. Hopefully something in the following code will be inspiring, even if I haven't solved the exact problem for you here:

```
library(GenomicRanges)

# make a GenomicRanges set corresponding to the original table
test.gr = makeGRangesFromDataFrame(df = test, start.field = "start", end.field = "end", seqnames.field = "chr")

# make separate GenomicRanges sets for each of the two partitions of the genome
test_A.gr = makeGRangesFromDataFrame(df = test[test$comp=="A",], start.field = "start", end.field = "end", seqnames.field = "chr")
test_B.gr = makeGRangesFromDataFrame(df = test[test$comp=="B",], start.field = "start", end.field = "end", seqnames.field = "chr")

# apply the reduce() function to concatenate segments which are adjoining each other into single segments
reduced_test_A.gr = reduce(test_A.gr)
reduced_test_B.gr = reduce(test_B.gr)

# group PC1 values by segment classification
test[queryHits(findOverlaps(test.gr, reduced_test_A.gr)), "PC1"]
test[queryHits(findOverlaps(test.gr, reduced_test_B.gr)), "PC1"]

```
> 2 votes

# Answer
Take a look at the `rle` function:

> Compute the lengths and values of runs of equal values in a vector – or the reverse operation.

In this way you can mark the moment in which A changes to B and the length of the streches of A and B.

With that I think you can reduce the number of for loops you have to use (I can not think a way of doing this without at least one loop, but for sure it may exist)
> 0 votes

# Answer
What is proposed by @Jonathan works, just adding another possible solution here.

library(GenomicRanges)

test.gr = makeGRangesFromDataFrame(test,keep.extra.columns=TRUE)
# as long as the previous comp does not agree with current, we make it a new segment
test.gr$seg = c(0,cumsum(diff(as.numeric(test.gr$comp)) != 0))

# simple function to collapse the range
fn = function(u){
  out = reduce(u)
  out$comp = unique(u$comp)
  out$mean = mean(u$PC1)
  out
}

reduced_range = unlist(GRangesList(lapply(split(test.gr,test.gr$seg),fn))))

```

> 0 votes

---
Tags: r, bioconductor, data-management, genomicranges
---


# Question
Title: 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?

# Answer
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(?)
> -1 votes

---
Tags: r, statistics
---


# Question
Title: 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?

# Answer
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
> 1 votes

---
Tags: sequence-annotation
---


# Question
Title: How to efficiently get human gene names from NCBI based on a large list of SNPs
I found a good answer related to my question here: How to get a list of genes corresponding to the list of SNPs (rs ids)?

But it says about small number of SNPs. I want to get gene names based on approximately 12000 to 19000 SNPs. NCBI says about 10 read per second with API key. Is there any efficient way to get gene names for thousands of SNPs? Is that database available to store locally? Any advice please?

# Answer
> 1 votes
Modify SNPs and gene annotations for your genome assembly of interest. This example shows how one might do this for `hg19`, for instance:

1. Download SNPs and convert to BED:

```
$ wget -qO- ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh37p13/VCF/common_all_20180423.vcf.gz \
    | gunzip -c \
    | convert2bed --input=vcf --output=bed --sort-tmpdir=${PWD} - \
    > hg19.snp151.bed

```
2. Filter this BED file for entries of interest:

```
$ grep -wFf snps-of-interest.txt hg19.snp151.bed | cut -f1-6 > hg19.snp151.filtered.bed6

```
3. Download gene annotations and convert to BED:

```
$ wget -O - ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz \
    | gunzip -c \
    | awk '($3 == "gene")' \
    | gtf2bed \
    > gencode.v19.genes.bed

```

***Note****: In the future, newer versions of Gencode annotations may be available for your reference assembly. Always check their site to confirm what is most recent.*
4. Map genes to SNPs:

```
$ bedmap --echo --echo-map-id hg19.snp151.filtered.bed6 gencode.v19.genes.bed > snps_with_associated_gene_names.bed

```

A list of genes that overlap a SNP of interest will be in the seventh column of the output.

If this must be done within Python, you could use something like `subprocess.call(...)`, where `...` includes commands wrapped in `bash -c`.

I don't recommend this because of the need to escape quote marks and deal with Python's quirks in running command-line tasks, which makes a Python-based approach very difficult to set up and maintain. I'd suggest learning some basic `bash` scripting to make this easy and fast.

# Answer
> 0 votes
Run them through the Ensembl VEP. You'll get the gene name and all the information about how the genes are affected by the variants.

---
Tags: gene, snp, python
---


# Question
Title: 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:

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

# Answer
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.
> 1 votes

---
Tags: phylogenetics, phylogeny
---


# Question
Title: RNA-seq: How to get new expression count after normalization
I've RNA seq, Human, Paired-end data, Sample size is \<40. These are aligned using STAR, RSEM processed. With RSEM I've TPM and expected counts, that is two files columns as individual IDs and row as gene names.

I'm interested to normalize gene data. With edgeR tutorial (link in the end) and few other online resources I see that after following steps there's an R object that contains `norm.factors` (Page 15) value for each individual.

I'm unable to wrap my head around it, now. If I'm interested to get normalized gene counts, can I go ahead and multiple each individual's norm.factor into its gene counts? For example, the expected count for `IND1` for `Gene1` is 100 and it's `norm.factor` is `0.80`, can I say that the normalized gene count is `100*0.80=80`?

I'm not interested to perform voomm, or differential expression analysis.

Follow up question: If I'm to sample few individuals' data from these RNA-seq should I
a) normalize all together and extract interested individuals, or,
b) extract interested individuals and then perform normalization.
I think since normalization uses per person library size, each person's normalization is independent of the others. Can someone please correct me if I'm wrong?

https://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf

# Answer
> If I'm interested to get normalized gene counts, can I go ahead and multiple each individual's norm.factor into its gene counts? For example, the expected count for IND1 for Gene1 is 100 and it's norm.factor is 0.80, can I say that the normalized gene count is 100\*0.80=80?

No; you divide, but if you poke around, you can probably find a way to get edgeR to return the normalized counts.

> I think since normalization uses per person library size, each person's normalization is independent of the others.

Not true for the default normalization used in edgeR, TMM. CPM normalization is independently determined for each sample.

In general, normalize everything together, unless you think that the algorithm assumptions might be violated by including all of them together.
> 2 votes

---
Tags: rna-seq, normalization, edger
---


# Question
Title: Show presence of known mutation in RNA-seq data
We have RNA-seq fastq data from control (WT) patients and a patient with a point mutation at a known location in one gene. I'd like to retrieve the reads aligning to that gene and show the presence of the mutation.

I can think of 2 strategies:

## A)

Alignment and visualization, based on GATK:

1. Run STAR 2-pass
2. Extract alignments of interest, e.g.:

`samtools view input.bam "Chr10:18000-45500" > output.sam`
3. Visualize:

**How can I count the reads showing mutation vs wild type?**

---

## B)

The second approach would simply look for the count numbers mapping (e.g. salmon) to the transcript ID(s) corresponding to mutation and WT, based on this SE answer.

**Is there a better way to do this?**

# Answer
For counting reads I use mpileup, e.g. `samtools mpileup --reference hg38.fa -r Chr10:18000-45500 input.bam`, which will give base-resolution coverage for a BAM file.

I've written my own script to process mpileup output and make it easier to understand. By default it reports read coverage as a proportion of total coverage, but this can be modified by using the `-counts` command line argument:

`$ samtools mpileup -r tig00018708_tig00000379:210665-240664:2107-2130 --reference trimMmerged_tig00018708_other.fasta local_Sampled_50M_vs_trimMmerged_tig00018708_other.bam | /bioinf/scripts/readstomper.pl -counts`

```
[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000
Assembly,Position,Coverage,ref,cR,pR,A,C,G,T,d,i,InsMode
tig00018708_tig00000379:210665-240664,2107,10,C,10,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2108,11,G,11,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2109,11,T,11,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2110,11,T,11,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2111,12,G,12,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2112,14,G,14,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2113,13,C,13,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2114,13,G,13,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2115,13,C,2,0.15,0,0,0,11,0,0
tig00018708_tig00000379:210665-240664,2116,15,T,14,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2117,15,G,14,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2118,18,A,17,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2119,19,T,18,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2120,19,C,18,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2121,18,A,7,0.41,0,0,10,0,0,0
tig00018708_tig00000379:210665-240664,2122,18,G,18,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2123,18,T,18,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2124,18,C,18,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2125,18,C,9,0.50,0,0,0,9,0,0
tig00018708_tig00000379:210665-240664,2126,18,C,14,0.78,2,0,1,1,0,0
tig00018708_tig00000379:210665-240664,2127,18,C,18,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2128,18,T,15,0.83,0,0,3,0,0,0
tig00018708_tig00000379:210665-240664,2129,17,A,13,0.76,0,0,4,0,0,0
tig00018708_tig00000379:210665-240664,2130,16,A,9,0.56,0,6,1,0,0,0

```

Note that only *non-reference* coverage is reported in the A/C/G/T/d/i columns. I have found that this makes it easier for subsequent data analysis, but your mileage may vary.
> 3 votes

# Answer
If you only have one gene and you only need to do this once then the simplest possible workflow is to generate the aliment using STAR (optimally with the two pass method) and open the two resultant `.bam` files in IGV, the coverage column at the stop above your alignment track should clearly show the counts of reference and non-reference supporting reads. The IGV documentation shows a screen shot the pop-up coverage indicator and allele specific read count.
> 3 votes

# Answer
SCmut is a method to detect cell-level mutation from single-cell RNA-sequencing, published by Vu et al. (2019)

https://github.com/nghiavtr/SCmut
> 0 votes

---
Tags: rna-seq, sequence-alignment, snp, rna-alignment, mutations
---


# Question
Title: Can I use GSEA to look for enrichment in data that is grouped into sets but are not necessarily genes?
I have yeast growth data in various experimental conditions that are classified by the strain type. Can I use the GSEA tool from the Broad Institute to test for enrichment of specific groups of strains within the total population (i.e. gene sets)? If not, would something like a hypergeomtric test be useful for this kind of analysis?

# Answer
Yes, if you have a ranked list, you can use GSEA to calculate the enrichment of the strains along the ranked list.

Note that hypergeometric test is for count data, so you would test if there are more strains of one type that fit on certain conditions.

I would recommend this article comparing different enrichment methods to understand well what they do.
> 1 votes

---
Tags: gsea
---


# Question
Title: PDB file downloading: pymol automation vs. manual
I automated a PDB download using a Pymol script (below)

```
python
pdb_lists = ['3LZG', '6HJN'] # lots more pdbs
for x in pdb_lists:
    cmd.fetch(x)
    cmd.select(x)
    cmd.save(x + '.pdb', x)
    cmd.delete(x)
cmd.quit()
python end

```

When the script is run within Pymol it pulls down a collection of pdbs, however this flopped when I plugged it through a Biopython PDB parser, but was okay when I manually downloaded the same pdb file from RCSB.org

A manual downloaded PDB was perfect, meaning the amino acids which cannot be observed by X-ray crystalography (e.g. they not sufficiently stationary) were recognised as being absent from the parser (resulting in a protein ~390 amino acids), but the Pymol scripted download contained all amino acid residues, i.e. even those which couldn't be observed in the X-ray (resulting in a protein of ~440 amino acid residues). Any ideas why the script resulted in a pdb that flopped and how to correct it?

# Answer
> 4 votes
To expand on the comment by marcin:

1. `fetch` downloads files in mmCIF format by default (https://pymolwiki.org/index.php/Fetch). Not all PDB entries have PDB format files, e.g. due to too many chains. Presumably this is why the change was made, though mmCIF files tend to be larger and hence download slower.
2. When calling `save` with the `.pdb` file extension the structure is converted to PDB format. Perhaps the converter uses the SEQRES records to define the sequence, and hence gives bad files when there are some residues not resolved.

To correct it you could try:

1. Saving as `.cif` and reading in with the Biopython mmCIF parser.
2. Providing the `type=pdb` argument to `fetch`, which will force PDB format downloads. Some PDB files may not be available.
3. Don't do it in PyMol and instead use the Biopython structure downloader (http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc187).

Also, if download speed is limiting you then consider using the MMTF file format, which is a binary format and hence smaller. Biopython can read in and write out MMTF files.

# Answer
> 1 votes
So this was the code I used and it worked (as far as I could see) and refers to point 2 of @jgreener's advice and was based on the comment by @marcin

```
python
pdb_lists = ['3LZG', '6HJN'] # lots more pdbs
for x in pdb_lists:
    cmd.fetch(x, type='pdb')
    cmd.select(x)
    cmd.save(x + '.pdb', x)
    cmd.delete(x)
cmd.quit()
python end

```

There are some Bio SE members who have serious skills in Pymol scripting, and constructive critique would be very welcome.

---
Tags: python, phylogenetics, phylogeny, pdb, pymol
---


title: Mutation annotation to the heatmap

text:
# Question
Title: Mutation annotation to the heatmap
I have mutation data of patients.How do i add that mutation info or annotate those as column with my annotation information expression and mutation

To help myself i did try this Zuguang Gu but couldn;t able to do it, i would like similar to the one to show patient mutations . any help or suggestion would be highly appreciated .

I wanted something like this , as here they annotated with batch ,type and cell clone. I would like to have something similar where i would like to see how mutation in patients affects expression profile .

```
df <- read.csv("HSC_pHSC_LSC_data/heatma_cluster.txt",header = TRUE,sep = ",")
        head(df)
        tail(df)
        dim(df)
        heat <- t(scale(t(df[,3:ncol(df)])))

        head(heat)

        rownames(heat) <- df$Gene
        ColAnn <- data.frame(colnames(heat))
        Sample = c(rep("HSC",4),rep("pHSC",12),rep("LSC",8))
        anno <- as.data.frame(Sample)
        categories <- data.frame(anno)
        colnames(anno) <- c("Condition")
        ColAnn <- HeatmapAnnotation(df=anno,which="col", col=list(Condition=c("HSC"="black","pHSC"="red","LSC"="dodgerblue4")))
        RowAnn <- data.frame(df$Family)
        colnames(RowAnn) <- c("Gene family")
        colours <- list("Gene family"=c("Cluster_one"="red","Cluster_two"="blue","Cluster_three"="brown","Cluster_four"="yellow",
                                       "Cluster_five"="green","Cluster_six"="maroon","Cluster_seven"="pink","Cluster_eight"="navyblue",
                                       "Cluster_nine"="darkgrey","Cluster_ten"="violet"))

        RowAnn <- HeatmapAnnotation(df=RowAnn, col=colours, which="row")
        boxAnnCol <- HeatmapAnnotation(boxplot=anno_boxplot(heat, border=TRUE, gp=gpar(fill="#CCCCCC"),
                                                            pch=".", size=unit(2, "mm"),
                                                            axis=TRUE,  axis_param = list(side = "left")),
                                       annotation_width=unit(c(1, 6.5), "cm"))

        boxAnnRow <- rowAnnotation(boxplot=row_anno_boxplot(heat, border=TRUE, gp=gpar(fill="#CCCCCC"),pch=".",
                                                            size=unit(3, "cm"), axis=TRUE, axis_param = list(side = "bottom")),
                                   annotation_width=unit(c(3), "cm"))
        pdf("ALL1_CLUSTER_DISEASE.pdf",width = 10, height = 20)
        #png("HSC_CMP_DIFF_HEATMAP.png", units="in", width=20, height=10, res=100)

        #dev.off()

        myCol <- colorRampPalette(c("navyblue", "white", "red"))(100)
        myBreaks <- seq(-2, 2, length.out=100)
        hmap <- Heatmap(heat,
                        name="Z-score",
                        col=colorRamp2(myBreaks, myCol),
                        heatmap_legend_param=list(color_bar="continuous",
                                                  legend_direction="vertical", legend_width=unit(5,"cm"),
                                                  title_position="topcenter", title_gp=gpar(fontsize=10, fontface="bold")),

                        #Split heatmap rows by gene family
                        split=df$Family,
                        #split = NULL,

                        #Row annotation configurations
                        cluster_rows=TRUE,
                        show_row_dend=FALSE,
                        row_dend_side = "left",
                        #row_title="Transcript", #overridden by 'split' it seems
                        row_title_side="left",
                        row_title_gp=gpar(fontsize=8,fontface="bold"),
                        row_names_gp = gpar(fontsize = 8),
                        show_row_names=TRUE,
                        row_names_side="right",
                        row_title_rot=0,

                        #Column annotation configuratiions
                        cluster_columns=TRUE,
                        show_column_dend=TRUE,
                        column_title="Samples",
                        column_title_side="top",
                        column_title_gp=gpar(fontsize=10, fontface="bold"),
                        column_title_rot=0,
                        column_names_gp = gpar(fontsize = 20, fontface="bold"),
                        show_column_names=TRUE,

                        #Dendrogram configurations: columns
                        clustering_distance_columns="euclidean",
                        clustering_method_columns="complete",
                        column_dend_height=unit(10,"mm"),

                        #Dendrogram configurations: rows
                        clustering_distance_rows="euclidean",
                        clustering_method_rows="complete",
                        row_dend_width=unit(15,"mm"),

                        #top_annotation_height=unit(0.5,"cm"),
                        top_annotation=ColAnn,

                        #bottom_annotation_height=unit(2, "cm"),
                        bottom_annotation=boxAnnCol)
        hmap + RowAnn

        dev.off()

```

The image i get after running the code is something like this sample numbers vary, as i have only HSC ,pHSC and LSC samples

# Answer
You are trying to "map" mutation data to your heatmap, however, there is no mutation data in your code as far as I can see! You specify sample names within a data frame called `anno` but this data frame has only one column, `Condition`. All you have to do is add another column, `mutation`, to the `anno` data frame. Moreover, you have to pay attention to the order of your data and your meta data. The rownames of the meta data (`anno` data frame in your case) should be the same as the colnames of your gene expression data frame.

I also see that you create a row annotation object `RowAnn`, which is not used in the `Heatmap()` call. Similarly, there is a `categories` data frame which is not used either. Your code would be easier to read without these unused variables.
> 1 votes

---
Tags: r, visualization
---


title: Help Adjusting Bedtools Output in R

text:
# Question
Title: Help Adjusting Bedtools Output in R
I have been asked to write a script in R to edit, using data.frame, this .fa file that is the result from extracting gene sequences from a repeat masker coordinate file on the genome I am working with. Each input looks as such:

```
>rnd-5_family-5445_Unspecified(-)
TTTCACCGTAAATTAACTTTGAGAGGAGCTAATTCCTAAAAGAATTATACCGGCCTATTTG
>rnd-5_family-2614_Unspecified(-)
AATCTGAGATTAGATTTATTTATCTGTTATCCGCTCAGGTTGAAAAGTTTGCGCAATGTATAACTATA
AAAATCTTTCTACCACTTTTGATTCATTTTTTAAATCTGGGGTCACATTTTATTCACGGTAGAAATTG
TAATTTACAAAATGAATTACTTGAAGGCAACACGAATCCAGAGTGATGCTTTACATAAATCTGCTTCT
ACCGATGCCAAAAATTGACGATATTCTATTATTTAATCTAAATGTTAGTCTTTACATACCCTCCCCTA
ATTGTTAGAATTTTATGAAATTTGATTTCAGGGGTCAGTTTAGCATGCTAAATCTAATTCAATAGATT
GATATTTTTCTTCAGGTAAAGAAAATTTTTGCGTCAAAGTAATCATATTCCTCCACGATTGCATATAA
CTATGGTATATAATTTAAAAGATTACACTTTACGTAATGAAAAATCGGCCAATCATTCAAAAGTTATG
AAAGTGATC
>rnd-5_family-2614_Unspecified(-)
TTTCTAGTTTGTATTTGAGGAAAGAAAATCCTACTATCACTTTTTAAAAAATTCACCAATCTGAGATT
AGATTTATTTATCTGTTATCCGCTCAGGTTGAAAAGTTTGCGCAATGTATAACTATAAAAATCTTTCT
ACCACTTTTGATTCATTTTTTAAATCTGGGGTCACATTTTATTCACGGTAGAAATTGTAATTTACAAA
ATGAATTACTT
>rnd-4_family-798_Unspecified(-)
ATGCTTAATTTGCTAAGCATTTAATAACCTTTTTGGGTTTTGTGATAAAGGATGCTGTAGACATTAAA
ATAAACCTTATACTGCTAT

```

This is the repeat masker coordinates that the above were derived from:

```
>KB824701.1 417 478 rnd-5_family-5445_Unspecified . -
>KB824701.1 587 1072 rnd-5_family-2614_Unspecified . -
>KB824701.1 914 1129 rnd-5_family-2614_Unspecified . -
>KB824701.1 1138 1225 rnd-4_family-798_Unspecified . -

```

and ideally I would like it to be in the format of this, for example using the first one:

```
>rnd-5_family-5445_KB824701.1_417_478

```

I was told to use data.frame in R but I am not sure exactly how this will work? Can someone help please

# Answer
> 1 votes
In case an answer would still be relevant after a year and a half, below is a solution with the `Biostrings` package.

```
library(Biostrings)

# loading data.table for the fread() function.
# fread() is not only fast but is quite good at guessing how
# the file is formatted
library(data.table)

# readDNAStringSet() function creates a "DNAStringSet" object
my_fasta <- readDNAStringSet("fasta.fa")

> my_fasta
  A DNAStringSet instance of length 4
    width seq                                                                names
[1]    61 TTTCACCGTAAATTAACTTTGAGAGGAGCTAATTCCTAAAAGAATTATACCGGCCTATTTG      rnd-5_family-5445...
[2]   485 AATCTGAGATTAGATTTATTTATCTGTTATCC...GGCCAATCATTCAAAAGTTATGAAAGTGATC rnd-5_family-2614...
[3]   215 TTTCTAGTTTGTATTTGAGGAAAGAAAATCCT...TAGAAATTGTAATTTACAAAATGAATTACTT rnd-5_family-2614...
[4]    87 ATGCTTAATTTGCTAAGCATTTAATAACCTTT...GTAGACATTAAAATAAACCTTATACTGCTAT rnd-4_family-798_...

# the fasta header is stored at the @NAMES slot of the @ranges slot of the
# DNAStringSet object created above
> my_fasta@ranges@NAMES
[1] "rnd-5_family-5445_Unspecified(-)" "rnd-5_family-2614_Unspecified(-)"
[3] "rnd-5_family-2614_Unspecified(-)" "rnd-4_family-798_Unspecified(-)"

# removing the trailing "_Unpecified..."
my_fasta@ranges@NAMES <- sub("_Unspecified.*",
                              "",
                              my_fasta@ranges@NAMES)

my_coord <- fread("repeat_masker.txt")
> my_coord
           V1   V2   V3                            V4 V5 V6
1: KB824701.1  417  478 rnd-5_family-5445_Unspecified  .  -
2: KB824701.1  587 1072 rnd-5_family-2614_Unspecified  .  -
3: KB824701.1  914 1129 rnd-5_family-2614_Unspecified  .  -
4: KB824701.1 1138 1225  rnd-4_family-798_Unspecified  .  -
#removing the ">" sign

my_coord$V1 <- sub(">",
                   "",
                   my_coord$V1)

# appending first 3 columns of the coordinate file to "NAMES" as requested by the OP
my_fasta@ranges@NAMES <- paste(my_fasta@ranges@NAMES,
                               my_coord$V1,
                               my_coord$V2,
                               my_coord$V3,
                               sep = "_")

# writing the modified DNAStringSet object into a fasta file
writeXStringSet(my_fasta, "edited_fasta.fa")

>rnd-5_family-5445_KB824701.1_417_478
TTTCACCGTAAATTAACTTTGAGAGGAGCTAATTCCTAAAAGAATTATACCGGCCTATTTG
>rnd-5_family-2614_KB824701.1_587_1072
AATCTGAGATTAGATTTATTTATCTGTTATCCGCTCAGGTTGAAAAGTTTGCGCAATGTATAACTATAAAAATCTTTCTA
CCACTTTTGATTCATTTTTTAAATCTGGGGTCACATTTTATTCACGGTAGAAATTGTAATTTACAAAATGAATTACTTGA
AGGCAACACGAATCCAGAGTGATGCTTTACATAAATCTGCTTCTACCGATGCCAAAAATTGACGATATTCTATTATTTAA
TCTAAATGTTAGTCTTTACATACCCTCCCCTAATTGTTAGAATTTTATGAAATTTGATTTCAGGGGTCAGTTTAGCATGC
TAAATCTAATTCAATAGATTGATATTTTTCTTCAGGTAAAGAAAATTTTTGCGTCAAAGTAATCATATTCCTCCACGATT
GCATATAACTATGGTATATAATTTAAAAGATTACACTTTACGTAATGAAAAATCGGCCAATCATTCAAAAGTTATGAAAG
TGATC
>rnd-5_family-2614_KB824701.1_914_1129
TTTCTAGTTTGTATTTGAGGAAAGAAAATCCTACTATCACTTTTTAAAAAATTCACCAATCTGAGATTAGATTTATTTAT
CTGTTATCCGCTCAGGTTGAAAAGTTTGCGCAATGTATAACTATAAAAATCTTTCTACCACTTTTGATTCATTTTTTAAA
TCTGGGGTCACATTTTATTCACGGTAGAAATTGTAATTTACAAAATGAATTACTT
>rnd-4_family-798_KB824701.1_1138_1225
ATGCTTAATTTGCTAAGCATTTAATAACCTTTTTGGGTTTTGTGATAAAGGATGCTGTAGACATTAAAATAAACCTTATA
CTGCTAT

```

---
Tags: r, ngs, phylogenetics, sequence-annotation
---


# Question
Title: 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

# Answer
> 1 votes
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,

1. 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.
2. 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.

# Answer
> 1 votes
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.

# Answer
> 1 votes
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.

# Answer
> 1 votes
So this is what I did to produce the phylogenetic tree:

1. 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)
2. I downloaded all the hits in fasta format (very cool feature indeed!)
3. 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
4. I took the fasta blastp list and loaded it into Seqotron.
5. I aligned everything using the Compute - Alignment - MUSCLE option -- I didn't do any manual AA alignment (not too sure about this part)
6. I produced the tree using Compute - Phylogenetic Tree (default parameters: distance matrix: K83, resampling: none).
7. 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

# Answer
> 1 votes
* 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.

---
Tags: phylogenetics, phylogeny, duplications
---


# Question
Title: snakemake with --conda-prefix, can I use a pre-built conda environment?
I have read this question (Running Snakemake in one single conda env) but still have some doubt.

I am using the '--conda-prefix /some/dir' option and I have a rule in my snakemake file with the directive:

```
rule xxx:
  ...
  conda:
    "envs/some.yaml"
  ...

```

The first call of that rule will create the conda environment according to the yaml file in /some/dir/envs with some random name (in my case 'bf55a79e'). If I run the workflow in a different directory, the same conda environment (/some/dir/envs/bf55a79e) gets activated directly instead of being built again. So this is what I expected.

My question is:

* is there some way of skipping the building and specifying that rule xxx should use an already existing conda environment somewhere?
* does anyone know where the mapping between rule and directory gets stored?

# Answer
1. If the environment already exists it will simply be used. You cannot, however, tell it to use an environment with a normal name (e.g., "my\_env").
2. I assume you mean, "where the heck does the random looking conda environment name come from?!?", which is a very good question. The answer to that is that it's a hash of the `--conda-prefix` setting and the contents of the environment yaml file you specify. The actual code for doing this (I'm "stealing" it from snakePipes) is:

```
import hashlib
md5hash = hashlib.md5()
md5hash.update("what you gave to --conda-prefix".encode())
f = open("your environment yaml file", 'rb')
md5hash.update(f.read())
f.close()
h = md5hash.hexdigest()

```

`h` is then the directory name. There are actually two lengths of those, the shorter 8 character version that snakeMake uses by default or the full length version that the code above will produce. snakeMake will actually use either of these. If you're wondering, no I don't think any of this is documented anywhere, you have to dig through the snakeMake code to find it.
> 3 votes

---
Tags: snakemake, conda
---


# Question
Title: How to take into account alternative bwa mem mapping when computing coverage
when mapping short reads with bwa mem, if a read has alternative mapping positions they are reported by bwa mem in the X0 and XA tags. Now, let's say I want to compute the coverage of my bam file. How would I do for the alternative hits to be taken into accounts. What I mean is that, if a read QUERY has a selected best HIT position and 3 alternate HIT positions, I would like that QUERY read to be counted as 1 unit of coverage for the best HIT position, but also as 1 unit of coverage for each of the alternate HIT position.

My understanding is that samtools depth will only consider the best hit (which makes sense in 99.99% of the cases).

I guess I could edit the bam file with custom scripts so that the alternative hits are also reported as primary hits but I was wondering if there was an already existing tool (as I am not that confident with my ability to mess around with sam files).

Thanks a lot.

# Answer
> 1 votes
You would want to instead use the `-a` option with `bwa mem` to output all alignments. You would be best off filtering these with a custom script to only keep those with a close enough score to the primary hit. `samtools depth` can't be told to include secondary alignments (it's hard-coded), so you'll need to use something else like `bamCoverage` or `samtools mpileup`.

---
Tags: samtools, bwa, coverage
---


# Question
Title: DMPfold error while running protein structure prediction
I am using DMPfold to predict protein structure of a 382 residue protein sequence. All goes wel, and even the files are produced but during the "seq2maps.csh", I get this error in between:

```
Running PSIPRED & SOLVPRED
Running PSI-BLAST with sequence c600m2-tail-before-lysin-sequence.temp.fasta ...
Segmentation fault (core dumped)
[blastpgp] WARNING: Unable to open BLOSUM62
[blastpgp] WARNING: BlastScoreBlkMatFill returned non-zero status
[blastpgp] WARNING: SetUpBlastSearch failed.
BLASTP 2.2.17 [Aug-26-2007]

Reference: Altschul, Stephen F., Thomas L. Madden, Alejandro A. Schaffer,
Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997),
"Gapped BLAST and PSI-BLAST: a new generation of protein database search
programs",  Nucleic Acids Res. 25:3389-3402.

Reference for composition-based statistics:
Schaffer, Alejandro A., L. Aravind, Thomas L. Madden,
Sergei Shavirin, John L. Spouge, Yuri I. Wolf,
Eugene V. Koonin, and Stephen F. Altschul (2001),
"Improving the accuracy of PSI-BLAST protein database searches with
composition-based statistics and other refinements",  Nucleic Acids Res. 29:2994-3005.

Query= C600M2_00100 tail protein [Salmonella phage BPS17S6]
         (377 letters)

Database: c600m2-tail-before-lysin-sequence.a3m
           46,978 sequences; 11,323,713 total letters

SearchingFATAL: Error whilst running blastpgp - script terminated!
Running PSICOV
Running FreeContact
Running CCMpred
 _____ _____ _____               _
|     |     |     |___ ___ ___ _| |
|   --|   --| | | | . |  _| -_| . |
|_____|_____|_|_|_|  _|_| |___|___|
                  |_|

using CPU (10 thread(s))

Reweighted 30404 sequences with threshold 0.8 to Beff=22917.9 weight mean=0.753779, min=0.0232558, max=1

Will optimize 62686429 32-bit variables

iter    eval    f(x)        ║x║         ║g║         step
1       1       1.50338e+07 152357      3.021629e+11    1.55e-06
2       1       1.47369e+07 152356      2.1447226e+11   1.18e-06
3       1       1.44292e+07 152355      1.4646483e+11   1.22e-06
4       1       1.41268e+07 152353      1.0388224e+11   1.5e-06
5       1       1.38295e+07 152358      7.337667e+10    1.92e-06

```

Can you please let me know if I need to rerun seq2maps.csh because of this error? And, how this error can be solved? Thank you!

# Answer
> 1 votes
The comment by Devon Ryan is correct - looks like your `blastpgp` command is not working properly. For a protein of that length it's more likely to be a problem with the BLAST executable running on your system than a memory issue.

Try running the `blastpgp` command (https://github.com/psipred/DMPfold/blob/master/bin/runpsipredandsolvwithdb#L26) on its own and see what happens:

```
~/path/to/blastpgp -a 12 -b 0 -j 3 -h 0.001 -d yourfile.a3m -i yourfile.fasta

```

# Answer
> 0 votes
Thanks, I did this but got this erro:

```
blastpgp -a 12 -b 0 -j 3 -h 0.001 -d c600m2-tail-before-lysin-sequence.a3m -i c600m2-tail-before-lysin-sequence.temp.fasta
[blastpgp] WARNING: Unable to open BLOSUM62
[blastpgp] WARNING: BlastScoreBlkMatFill returned non-zero status
[blastpgp] WARNING: SetUpBlastSearch failed.
BLASTP 2.2.26 [Sep-21-2011]

Reference: Altschul, Stephen F., Thomas L. Madden, Alejandro A. Schaffer,
Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997),
"Gapped BLAST and PSI-BLAST: a new generation of protein database search
programs",  Nucleic Acids Res. 25:3389-3402.

Reference for compositional score matrix adjustment: Altschul, Stephen F.,
John C. Wootton, E. Michael Gertz, Richa Agarwala, Aleksandr Morgulis,
Alejandro A. Schaffer, and Yi-Kuo Yu (2005) "Protein database searches
using compositionally adjusted substitution matrices", FEBS J. 272:5101-5109.

Reference for composition-based statistics starting in round 2:
Schaffer, Alejandro A., L. Aravind, Thomas L. Madden,
Sergei Shavirin, John L. Spouge, Yuri I. Wolf,
Eugene V. Koonin, and Stephen F. Altschul (2001),
"Improving the accuracy of PSI-BLAST protein database searches with
composition-based statistics and other refinements",  Nucleic Acids Res. 29:2994-3005.

Query= C600M2_00100 tail protein [Salmonella phage BPS17S6]
         (377 letters)

Database: c600m2-tail-before-lysin-sequence.a3m
           46,978 sequences; 11,323,713 total letters

SearchingSegmentation fault (core dumped)

```

Could you please let me know how this could be solved? Thanks

---
Tags: protein-structure
---


# Question
Title: Does this FASTQ data contain single or paired end calls?
I have this fastq data from GEO:

```
zcat SRR1658526.fastq.gz | head -n 20
@SRR1658526.1 HWI-ST398:296:C1MP4ACXX:1:1101:1093:2094 length=102
GATCTCTATTACTTTTTGAAGGATTNNNNNNNNNNAANTTTTGAATCANNNNNNNNNNNNNNNNNNNNNNNNNNNNNTNNNNNNNNNNNNNNNNNNNNNNNN
+SRR1658526.1 HWI-ST398:296:C1MP4ACXX:1:1101:1093:2094 length=102
@<@FFDFDHHFHHIIIIGHBGGGIG##########10#0:BGDDHHII######################################################
@SRR1658526.2 HWI-ST398:296:C1MP4ACXX:1:1101:1167:2107 length=102
ATAATATTGTAGATATAAATGTTATCTAATCTTATCTGATCAGCTTGCTNNATANNNNNNNNNNNNNNNNNNACNTATGNNNNNNNNNNNNNNNNNNNNNCC
+SRR1658526.2 HWI-ST398:296:C1MP4ACXX:1:1101:1167:2107 length=102
CCCFFFFFHHFHHJJIIIIJIJJHJHJJJJIJJIJJJJIJHIJJJJGII#####################################################

```

It is supposed to be paired-end sequences. Is the prefix `@` and `+` are the R1 and R2? What's the convention here?

# Answer
> 4 votes
Entries in a fastq file occupy 4 lines each and R1 and R2 are typically in different files. Since that SRA project is 2x51 you seem to have run `fastq-dump` without the `--split-3` option, so both R1 and R2 are merged together. Make your life easier and never use SRA, but instead just download the individual files from ENA.

# Answer
> 2 votes
The data shown in the fastq files are paired end data.

In fastQ files @ and + do not correspond to R1 and R2.

The wikipedia page for fastq files describe the following conventions:

> A FASTQ file normally uses four lines per sequence.
>
> ```
> Line 1 begins with a '@' character and is followed by a sequence identifier and an optional description (like a FASTA title line).
> Line 2 is the raw sequence letters.
> Line 3 begins with a '+' character and is optionally followed by the same sequence identifier (and any description) again.
> Line 4 encodes the quality values for the sequence in Line 2, and must contain the same number of symbols as letters in the sequence.
>
> ```

**When using fastq-dump ".1" and ".2" stand for read suffices assuming u used the -I flag.**

In this case you can see that R1 (".1") and R2 (".2") are within the same file.

---
Tags: fastq
---


title: Understanding Single Cell RNAseq Plots

text:
# Question
Title: Understanding Single Cell RNAseq Plots
I have just started some single cell RNAseq analysis with 500 cells. I have used Seurat for analysis. Following the guideline I have generated few plots. Now I am trying to understand how to read those plots. Can the community provide their input in the following plots. To give a brief background of the experiment. This data was processed in MinION. Its from human brain sample from our collaborator. In their shortseq data they identified 15 clusters. However, both the libraries were processed independently. The main question I am trying to answer is celltype specific isoforms. but in my dataset I can find only 6 Cell clusters. So, I am trying to understand why I am getting such a small cluster? Is it because I stared with only 500 Cells? Also I want to know what these two PCA plot is indicating. Thanks

# Answer
Identifying 15 clusters from 500 cells is challenging. It will depend on the relative abundance of each cell type within the tissue. You can use the cell type frequencies found in the shortSeq data set to see what to expect from a 500 cell dataset using the howmanycells tool

About the pc plots: The first is an elbowplot which explains how many variation is explained by each PC. The knee basically defines when components become less informative.

**EDIT**: I noticed that the y-axis of the elbowplot starts at 4 suggesting that the dimensionality of your data is most likely ~10 dimensions. But I expect many of the clusters that don't show up on your Umap plot to be present at very low frequencies.

The second plot is a jackstraw plot and provides a visualization of a statistical test that calculates the significance of your principal components. This one shows that all your components are significant. You can combine this with the elbowplot to estimate the dimensionality (10-15 dimensions) of your data and decide which components to use for downstream analysis. In general the more they dropoff to the dotted line the less significant they become.

While these plots are useful as a first impression, always rely on biological information and common sense to decide what information (here: components) to include in downstream analysis such as Umap, clustering, etc. For example: usually you don't want to include components mainly defined by ribosomal, mitochondrial or hsp genes.
> 3 votes

# Answer
The first two plots are described on the Seurat PBMC tutorial. Specifically:

> ‘Elbow plot’: a ranking of principle components based on the percentage of variance explained by each one (ElbowPlot function).

Based on the tutorial, you should use at least 10 PCs.

> The JackStrawPlot function provides a visualization tool for comparing the distribution of p-values for each PC with a uniform distribution (dashed line). ‘Significant’ PCs will show a strong enrichment of features with low p-values (solid curve above the dashed line).

Based on the tutorial, you should use more than 15 PCs.

No one can tell you the right number of clusters. If you already identified 15 clusters in a different dataset, then you should have cluster markers for those genes. You can check those genes in your dataset.
> 2 votes

# Answer
The first two plots are used in order to estimate the number of PCs to be used in later stages, for example clustering (this was split in to two functions in Seurat v3). The number of PCs selected would have an impact on the number of clusters obtained. The more PCs the more information for downstream applications. It would be a good idea to go over the genes within each PC, a trained eye can tell what kind of biology (or technical artefacts) contribute to the PCs.

Moreover, `FindClusters()` has a `resolution` parameter, which enables the number of communities (clusters) that is generated. In Seurat v3, there is a wide selection of different clustering algorithms but recommending one over the other would be beyond my expertise.

If you have a comparable dataset generated with some other technology/method, you can give a try with "integrating" the two datasets. An alternative non-Seurat approach that aims to do the same is `MNNCorrect` from the `scran` package.
> 1 votes

---
Tags: scrnaseq, seurat, figure-explanation
---


title: blast nt vs Genbank

text:
# Question
Title: blast nt vs Genbank
What is the difference between these two databases?

ftp://ftp.ncbi.nlm.nih.gov/blast/db/ and ftp://ftp.ncbi.nih.gov/genbank/

I understand that BLAST database has both nt for nucleotide and np for protein. But genbank is also for genomic database. So what is the difference between nt and Genbank databases?

Thanks

# Answer
> 2 votes
I think the BLAST folder contains pre-formatted databases for BLAST whereas the GenBank folder contains sequences.

---
Tags: blast, blastn
---


# Question
Title: How can we find the distance between all residues in a PDB file?
If we have a PDB structrure, how can we find residues physically interacting with each other in space? I know that we must find the distance between residues and if the distance is less than 5-6 Angstrom, we say that residues are physically interacting. But how can we find the distance between all residues and how can we finally determine the distances between all residues? Is there a software or webserver for that?

# Answer
If you need to process multiple files, you could use Biopython to parse a PDB structure.

```
from Bio.PDB import PDBParser

# create parser
parser = PDBParser()

# read structure from file
structure = parser.get_structure('PHA-L', '1fat.pdb')

model = structure[0]
chain = model['A']

# this example uses only the first residue of a single chain.
# it is easy to extend this to multiple chains and residues.
for residue1 in chain:
    for residue2 in chain:
        if residue1 != residue2:
            # compute distance between CA atoms
            try:
                distance = residue1['CA'] - residue2['CA']
            except KeyError:
                ## no CA atom, e.g. for H_NAG
                continue
            if distance < 6:
                print(residue1, residue2, distance)
        # stop after first residue
        break

```

If you need to look at one structure, using a viewer perhaps would be easier. You could try PyMOL: (how to measure distance). There are other PDB viewers, some of which can work even through a browser.
> 10 votes

# Answer
Could you use CCP4's NCONT program? There's a GUI and a command line interface, whatever suits. You can specify which chains you want to target and interact with and set a cut off for distance. The bonus here is once you're in you have a nice suite of other structural tools to use.

If you're just doing it once, the GUI is friendly enough to work things out, if you're doing a batch then you can run it across several files via the command line.

Download
> 4 votes

# Answer
As part of a project me and some teammates did a script that outputs visual maps of distances between residues. It uses Biopython.

The module contact\_map.py does what you are looking for. As an example, if you want to find the residues whose CA are below 5 you can run the following command:

```
python3 contact_map.py pdb1cd8.ent -a CA -CA 5

```

This will produce three files:

```
distance_map_pdb1cd8_CA.png # Heatmap of the distance between the residues
contact_map_pdb1cd8_CA.png # Black/White heatmap: If it is at that min distance
contact_map.log  # The actions taken

```

If you don't have downloaded already the pdb structue you can use the main module:

```
python3 cozmic.py real 1cd8 -a CA -CA 5

```
> 2 votes

# Answer
You can use MDtraj. The package is easy to install using Anaconda.

You can get the interacting residues using the following snippet (taken from http://mdtraj.org/1.6.2/examples/native-contact.html)

```
heavy_pairs = np.array(
    [(i,j) for (i,j) in combinations(heavy, 2)
        if abs(native.topology.atom(i).residue.index - \
               native.topology.atom(j).residue.index) > 3])

# compute the distances between these pairs in the native state
heavy_pairs_distances = md.compute_distances(native[0], heavy_pairs)[0]
# and get the pairs s.t. the distance is less than NATIVE_CUTOFF
native_contacts = heavy_pairs[heavy_pairs_distances < NATIVE_CUTOFF]
print("Number of native contacts", len(native_contacts))

```
> 0 votes

# Answer
Distance between C-alphas can also be found using the BioStructures.jl package in Julia:

```
using BioStructures
struc = read("1AKE.pdb", PDB)
calphas = collectatoms(struc['A'], calphaselector)
dm = DistanceMap(calphas)

```

Then `dm[5, 10]` gives `12.19677604943208` etc.

Other variations also work, see the docs.
> 0 votes

---
Tags: protein-structure, pdb
---


title: how to Know the longest ORF in Protein?

text:
# Question
Title: how to Know the longest ORF in Protein?
I have a protein Seq. and I want to know

* how to Know the longest ORF in this Protein?

and if there is any tutorial for this process ??

# Answer
Proteins don't have ORFs, DNA sequences do. What DNA sequences do you have? If cDNA (i.e. transcript/mRNA) then an ORF finder will give multiple ORFs and the coordinates on the full-length input DNA which indicate it's length.
> 2 votes

---
Tags: proteins, orf
---