title: 5' and 3' bias in Rna-seq data

text:
# Question
Title: 5' and 3' bias in Rna-seq data
I'm working with rna-seq samples. I see 5' bias and also 3' bias in the per-base sequence content plot. From this link I see that the bias at the start of the sequences appears to be the result of biased selection of fragments from the library. And this won't effect any downstream analysis. But In all of my samples I also see bias at the end of the sequences. What can be the problem? Samples are strand-specific sequenced with poly-A selection protocol.

In the Transcript coverage profile I see 5'-3' bias for all my samples is: 0.9

Showing the per-base sequence plot for one of the sample: For all the samples I see that 5' and 3' bias.

When I looked at Adapter content this is what I see:

No samples found with any adapter contamination \> 0.1%

# Answer
> 7 votes
I expect there was a sequencing problem during the last base, where some of the reagents were running low on the sequencer. This won't pose any real problem, RNAseq aligners like STAR will just soft-clip the last base or two if they're mismatches.

It's common to see a bit of bias toward the 5' or 3' ends in RNAseq, mostly due to whether poly-A selection was done or there was a bit of degradation. The primary thing to ensure is that the bias is similar across samples/groups. If that's the case then you don't have to worry that your downstream analyses are going be affected by some sort of per-group bias.

# Answer
> 3 votes
I know this is an old post, but if this plot is from *after* trimming I would suggest a different explanation: some trimming tools remove poly-A sequences from reads. If that's the case then any read ending with `A` will have that removed, this leads to a 0% `A` base content in the final base position (and a dropping %`A` in the final couple of bases).

This should have no effect on alignments or analysis.

---
Tags: rna-seq, r, fastq, rna, rna-alignment
---


# Question
Title: How to incorportate RIN values as covariate in the design matrix?
I have been following the last DESeq2 pipeline to perform an RNAseq analysis with a dataset with low rin samples in the experimental (or treated) and high rin on the control ones.

I read a paper in which they perform RNAseq analysis with time-course RNA degradation and conclude that including RIN value as a covariate can mitigate some of the effects of low RIN in samples.

My question is how I should construct the design in the DESeq2 object:

```
~conditions+rin
~conditions*rin
~conditions:rin

```

none of them... :)

I cannot find proper resources where explain how to construct these models (I am new to the field...) and I recognise I crashed against a wall with these kinds of things. **I would appreciate also some links to good resources to be able to understand which one is correct and why**.

Thank you in advance!

# Answer
In the DESeq2 manual there is a section titled "How can I include a continuous covariate in the design formula?" that deals with your question.
Basically the process is no different from using a discrete covariate.

```
DESeq <- DESeqDataSetFromMatrix(countData = counts, colData = metadata, design = ~ RIN + condition)

```

In this example condition is the experimental groups and RIN is the continuous covariate for the RNA-integrity.

For this to work properly there has to be a linear relationship between the RNA integrity and the gene expression values. As stated in the manual: "Continuous covariates might make sense in certain experiments, where a constant fold change might be expected for each unit of the covariate."
Therefore, it may be better to code the RIN variable as a factor like `RIN_high` or `RIN_low`. This is more robust than using the absolute values of the variable itself.

However, from your desription I get the impession that the RIN values are generally lower in the treatment group compared to the control group. If this is the case, using RIN as a covariate will by of limited use. DESeq2 will not be able to discriminate if a expression difference is due to a difference in RNA integrity or due to your treatment. Also, in this case coding RIN into high and low would result in an error because the RIN and treatment variable would define the same groups.

Generally, covariates work best when they are evenly distributed across the experimental groups.
> 4 votes

---
Tags: rna-seq, deseq2, edger, design
---


# Question
Title: How to specifically select cluster on Featurescatter on Seurat?
Nowadays, I am trying to calculate Pearson correlation values between two genes of my interest from single-cell RNA-data (features.tsv, barcodes.tsv, and matrix.mtx files) which are obtained from the 10xgenomics cellranger.

FYI, The filtered\_feature\_bc\_matrix includes an aggregated population of tumor-infiltrating immune cells that are derived from 6 different syngeneic mouse tumors (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE121861). I downloaded fastq files and ran "cellranger count" and "cellranger aggr" commands to obtain combined data.

My strategy includes the use of the "feature scatter" function deposited in the library(Seurat). Finally, I ended up getting results below.

My question is that can we select a specific cluster(s) or identity on this plot? For example, I want to see correlation values in cluster 2 only. The usage of Feature scatter is as follows.

FeatureScatter( object, feature1, feature2, cells = NULL, group.by = NULL, cols = NULL, pt.size = 1, shape.by = NULL, span = NULL, smooth = FALSE, combine = TRUE, slot = "data" )

Do you know how to get the Pearson correlation of my two genes using Seurat? This is the best I can do.

Thank you.

# Answer
You can use the `cells` argument to only display the data for a specific subset of cells.

For example to extract the values for a cluster 2:

```
# The cluster identity is stored in the meta.data
# In this example there is a clustering column called "RNA_snn_res.0.7"
# but any metadata column can be used

clust2.cells <- rownames(Seurat.data@meta.data[Seurat.data$RNA_snn_res.0.7 == 2, ])

FeatureScatter(Seurat.data, feature1 = "gene_A", feature2 = "gene_B", cells = clust2.cells)

```
> 3 votes

---
Tags: seurat, correlation
---


# Question
Title: How can I extract all known mutations of my BAM (or SNP/INDL files)?
I am using a Genome Explorer tool to see all the mutations on my own DNA.

My particular interest is on listing the variants and get their names/ids.

Here are a few screenshots:

You can see that they tell you the genomic position, the allele frequency of the mutation, the name of the disease, the Category, etc etc...

However, the filters on this tool work very badly so I cannot filter all my mutations only...

So I was trying to find out an alternative tool to get this data in such a nice way this tool provides. I tried IGV but they just provide raw readings, I didnt find a way to get a list of variants the same way this Genome Explorer does.

The closest thing I've found was Geneious:

but it crashes on my computer and is too complex to use (I wasn't able to reach that screen above, from my BAM file). And it also eats a lot of RAM that I don't have.

Are there any alternatives? Even paying?

# Answer
> 1 votes
If you have a BAM file, you could easily run a variant caller analysis software such as freebayes, though you would need a good amount of computing power like a cluster or a dedicated node with good memory. Freebayes outputs a vcf file that has all the variants listed. Then you can load it into the UCSC or Ensembl human genome browsers and look at what genes they are.

---
Tags: snp, ngs, variant-calling, phylogenetics, indel
---


# Question
Title: How can i modify FASTA headers in a multi fasta file using BioPython SeqIO
I have a multi fasta file similar to this (relatively new here so uncertain of best way to present this; I have gone for an output and the code i used to make it - belt and braces...):

```
# create a test dataset
with open('test_dat.txt', 'w') as test_dat:
    l1 = [">AAA a aa", ">BBB b bb", ">CCC c cc"]
    l2 = ["ATGATG", "GACGAC", "TTGTTG"]

    for i in range(len(l1)):
        test_dat.write(l1[i] + "\n")
        test_dat.write(l2[i] + "\n")

# which outputs something like this
>AAA a aa
ATGATG
>BBB b bb
GACGAC
>CCC c cc
TTGTTG

```

I want to modify the fasta headers in my fasta file by adding some text to the beginning and, ideally, stripping out spaces (as i was advised that white space in the headings could cause problems for me down the line). My plan was to do this using BioPython and SeqIO.

I started like this:

```
# my attempt at altering the fasta headers using SeqIO
from Bio import SeqIO

# open my test data and create an output file
with open('test_dat.txt', 'r') as inputs, open('test_dat_out.txt', 'w') as outputs:

    # create dictionary of the fasta sequences using SeqIO
    my_dict = SeqIO.to_dict(SeqIO.parse(inputs, "fasta"))

    # renaming my fasta headers - adding "TEST~~~TEST~~~" to front of header and remove white space
    for v in my_dict.values(): # using the values() to loop through...
        v.description = "TEST~~~TEST~~~" + v.description.replace(" ", "_")

    # write result to file
    SeqIO.write(my_dict.values(), outputs, 'fasta')

```

which gave me a result with some duplication in the fasta header:

```
#result i want
>TEST~~~TEST~~~AAA_a_aa
ATGATG

#result i get
>AAA TEST~~~TEST~~~AAA_a_aa
ATGATG

```

If i modify both the `.id` and `.description` in exactly the same way, as suggested here, but **without** removing white space:

```
# renaming my fasta headers - adding "TEST~~~TEST~~~" to front of header
    for v in my_dict.values(): # using the values() to loop through...
        v.id = "TEST~~~TEST~~~" + v.id
        v.description = "TEST~~~TEST~~~" + v.description

```

I get almost what i want, but the header still has spaces:

```
>TEST~~~TEST~~~AAA a aa
ATGATG

```

I have tried to make the `.id` and `.description` match by adding a "\_" to the end of `.id`, or by adding `.replace(" ", "_")` to both `.id` and `.description` (despite there being no spaces in `.id`), but I just cant seem to get exactly the output i want.

I would really appreciate any pointers on this.

# Answer
> 3 votes
You're very close, but unfortunately for legacy reasons this aspect of the Biopython parser is not very intuitive. The key point is the `.id` is anything up until the first white space character. One day I hope the `.description` and `.id` attributes are replaced by a single attribute like `.header`. Anyway, here is a solution for your case:

```
from Bio import SeqIO

to_add = "TEST~~~TEST~~~"

with open("out.fa", "w") as outputs:
    for r in SeqIO.parse("in.fa", "fasta"):
        r.id = (to_add + r.description).replace(" ", "_")
        r.description = r.id
        SeqIO.write(r, outputs, "fasta")

```

Output on your sample input:

```
>TEST~~~TEST~~~AAA_a_aa
ATGATG
>TEST~~~TEST~~~BBB_b_bb
GACGAC
>TEST~~~TEST~~~CCC_c_cc
TTGTTG

```

---
Tags: fasta, python
---


# Question
Title: about the normalization of RNAseq data for calculating distance for unsupervised learning
I have been working on clustering using RNAseq data. To compute distance, what kinds of normalization is optimal?

Can we use normalization using relative log normalization? I understand this should be used for differential gene expression analysis, not for unsupervised learning in which we do not use information of disease condition.

If you can show good papers, please share.

# Answer
> 2 votes
Typically some sort of variance stabilizing transform is used before clustering. Popular options are regularized log transformation or a vst transform, which are available in DESeq2. Note that these are NOT used for performing differential expression, just things like clustering and PCA. For differential expression one would use TMM, RLE or similar.

---
Tags: rna, clustering
---


# Question
Title: How to remove zero value of gene on FeatureScatter plot using Seurat?
Nowadays, I am trying to calculate Pearson correlation values between two genes of my interest from single-cell RNA-data (features.tsv, barcodes.tsv, and matrix.mtx files) which are obtained from the 10xgenomics cellranger. I am interested in tumor-infiltrating Tregs. I downloaded fastq files from a published paper that describes tumor-infiltrating Tregs with IL10+IL35+ (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE126184), and ran "cellranger count".

Yesterday, I learned how to get a Pearson correlation value from a kind guy. I made a graph below.

In general, Tregs are defined by CD4 and Foxp3 expression. As a trial, I want to find the correlation between two genes. The correlation value was 0.05, which might have been due to the presence of Tregs that have zero value of either CD4 or Foxp3. So, my question is that can we remove cells with 0 value of genes to get more high correlation value? My data looks like this.

```

> Treg.data[c("Cd4", "Foxp3"), 1:30]

------------------------------------- # results

2 x 30 sparse Matrix of class "dgCMatrix"
   [[ suppressing 30 column names ‘AAACCTGAGAAACCAT-1’, ‘AAACCTGAGAAACCGC-1’, ‘AAACCTGAGAAACCTA-1’ ... ]]

Cd4   . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Foxp3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

------------------------------------- # results

```

Each dot represents 0 value.

Meanwhile, I performed to standard pre-processing workflow, as suggested in Seural Tutorial. For example,

```
> Treg[["percent.mt"]] <- PercentageFeatureSet(Treg, pattern = "^MT-")

> Treg <- subset(Treg, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

```

Can you help me? Thank you very much.

# Answer
It will not improve your correlation most likely, unless the zeros are not missing at random:

```
library(Seurat)
data <- pbmc_small
data[["percent.mt"]] <- PercentageFeatureSet(data, pattern = "^MT-")
data <- NormalizeData(data)
genes <-c("HLA-DRA","LYZ")
countdata = pbmc_small[["RNA"]]@data
countdata <- t(as.matrix(countdata)[genes,])
plot(countdata)

```

Now you need to subset the matrix, and remove the rows where at least one of the entries is 0, you can use rowSums(..==0)==0, meaning no zeros in the rows:

```
plot(countdata[rowSums(countdata==0)==0,])

```

```
cor(countdata[rowSums(countdata==0)==0,])
          HLA-DRA       LYZ
HLA-DRA 1.0000000 0.2338371
LYZ     0.2338371 1.0000000

```
> 1 votes

---
Tags: r, scrnaseq, seurat, correlation
---


# Question
Title: Phylogenetics software that can represent collapsed nodes as triangles
I am trying to represent a large (100-300+ sequences) phylogenetic tree in a way that gets my point across. I figured that reducing the number of branches (by collapsing clades which include highly similar sequences) would be a good idea. The software that I'm currently using (iTOL) technically has an option for that ('collapse node'), but I'm not exactly a fan of the particular graphical representation, so I'd like to have a look at other options. Thank you for your time!

# Answer
For Mac Figtree here is highly recommended, it is also available for Windows and Linux here. The technique is called "collapsing clades" and is the only package I know that does it. I would recommend using it, because I would assume it is averaging the branch length of a given clade and the 'base' of the triangle is proportional to the number of taxa. Thus it recognised as being very reliable, albeit the precise algorhithm it uses ... I don't know :-O

With FigTree you get two different "triangle options" in the menu, which give different shaped triangles for collapsing the number of lineages within the select clade you chose. The big one, is the algorhithm appropriately described above. It has many options for producing publication quality diagrams, albeit I then import the output into a graphics program for further editing (there are shortcomings IMO in colouring and fonts. Then again for free is better than good.

If you are outputing a tree within a script, ETE3 (Python) will do a "collapsing of clade" option, but it isn't pretty triangle stuff. ETE3 will produce a publication quality tree diagram from this and reduce to 'low resolution', which is very useful in large analyses of loads of trees.

In summary, I don't think FigTree has any competition in the context you are seeking.

Moreover, you can configure your bin, or alias file so you can automatically open a treefile just typing,

```
figtree treefile.tre

```

Also you can reduce you font size to next point 1 or 2 (okay you can do that in anything), if you are attempting optical analysis on loads of taxa. Its also really stable. I just find it really useful, like a "jazz standard" its a "phylogeny standard".
> 2 votes

---
Tags: r, python, phylogenetics, phylogeny, visualization
---


# Question
Title: Calculating polarity and hydrophobicity of a ligand/small molecule
I have lot of SMILES IDs of the ligands. What I am interested in is the polarity and hydrophobicity score of each of the ligands.

In other words, does ligand have more polar groups and can it make more polar interactions or is it more hydrophobic?

Does anybody know the way/software to look at it?

Best

# Answer
> 1 votes
You can calculate their polar surface area or logo using rdkit or openbabel. Other molecular properties may be available as well.

---
Tags: database
---


# Question
Title: Where can I visualize Ramachandran plot?
I am using Rampage webserver to obtain a Ramachandran plot. However, the server does not reveal any plot just values (favored, allowed, outlier regions).

Can you recommend other similar web sites?

# Answer
Here is a page with several suggested methods. I agree that Rampage looks the best of the options, but it seems that it provides some suggested workflows for drawing it in R so that it is very customizable and you can probably get it to show any features you want if you're willing to customize a little.

Here is a Python tool for drawing Ramachandran plots.

There are additionally several other servers/tools.
> 4 votes

---
Tags: 3d-structure, validation
---


title: umap highlighting two different models

text:
# Question
Title: umap highlighting two different models
I'm trying to create a umap for single cell data from human samples and ptx samples. I can get the umap to where it shows the umap with the different clusters but I want to show where the ptx samples are and where the human samples are.

My code is as follows:

```
#create the Seurat object

OD_10K_HUMAN <- CreateSeuratObject(counts = HUMAN_OD_10K.data, min.cells = 0, project = "human")
SD_5K_HUMAN <- CreateSeuratObject(counts = HUMAN_SD_5K.data, min.cells = 0, project = "human")
BNL.5K <- CreateSeuratObject(counts = SD_BNL_5K.data, min.cells = 0, project = "ptx")
BNL.10K <- CreateSeuratObject(counts = OD_BNL_10K.data, min.cells = 0, project = "ptx")
BNM.10K <- CreateSeuratObject(counts = OD_BNM_10K.data, min.cells = 0, project = "ptx")
BNM.5K <- CreateSeuratObject(counts = SD_BNM_5K.data, min.cells = 0, project = "ptx")

#merge data
scData <- merge(BNL.10K, y = c(BNL.5K, BNM.10K, BNM.5K, SD_5K_HUMAN, OD_10K_HUMAN), add.cell.ids = c("A", "B", "C", "D", "E", "F"), project = "HTB2876")

mark the mito genes
mito.genes <- grep(pattern = "^MT-", x = rownames(x = scData), value = TRUE)
length(mito.genes)
scData[["percent.mt"]] <- PercentageFeatureSet(scData, pattern = "^MT-")

scData[["log_nCount_RNA"]] <- log2(scData[["nCount_RNA"]]+1)

# remove cells with <200 RNA molecules, or >6000 molecules, or >30% mito
scData <- subset(scData, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 30)
scData <- NormalizeData(scData, normalization.method = "LogNormalize", scale.factor = 10000)
scData <- FindVariableFeatures(scData, selection.method = "vst", nfeatures = 2000)

all.genes <- rownames(scData)
scData <- ScaleData(scData, features = VariableFeatures(object = scData), vars.to.regress = c("nCount_RNA"))
scData <- RunPCA(scData, features = VariableFeatures(object = scData))
#DimPlot(scData, reduction = "pca")

numPC = 20
scData <- FindNeighbors(scData, dims = 1:numPC)
scData <- FindClusters(scData, resolution = 0.4)
scData <- RunUMAP(scData, dims = 1:numPC)

DimPlot(scData, reduction = "umap", label = TRUE)

```

# Answer
> 1 votes
By default, `DimPlot` colors the cells by their "identity class", stored in `Ident(scData)` in your case. This indeed contains the clustering information.

However, you can you use the `group.by` argument for `DimPlot` to pass another column from `scData@meta.data` to color cells by. In your case you would need a metadata column that stores the *ptx*/*human* information.

The `merge` function automatically stores the original "project" information of the individual Seurat objects that you combine under `orig.ident` in the metadata. Since you set the project names to either `"ptx"` or `"human"`, the `orig.ident` column should contain the information you need.

Something like this should work:

```
DimPlot(scData, reduction = "umap", group.by = "orig.ident", label = TRUE)

```

---
Tags: scrnaseq, seurat, umap
---


# Question
Title: Outputting crunchable list of HBonds from Pymol
X-posted from the Main Stack Overflow because folks here are possibly more familiar with PyMol(?)

I have been trying to assess the strength of an interface I want to mutate using Pymol. What I am interested in is a list of contacts withing a specified interface (which I have calculated using the InterfaceResidues.py script from the PyMol wiki script library). My aim is to produce a bunch of mutants using the mutagenesis wizard and assess how the HBond network changes.

Actions\>Find\>Polar Contacts does the trick in visualising them, but doesn't output a file with numbers which I can then crunch to assess which mutations weaken my interface most. I have been trying the scripts list\_contacts.py and list\_hb.py from Prof. Robert Campbell of Queen's University at Kingston, but with limited success (nothing gets printed or written to the file). The scripts can be found here http://pldserver1.biochem.queensu.ca/~rlc/work/pymol/

I am going to contact Prof. Campbell directly re: scripts but, in the meantime, and in case he's too busy to deal with a Pymol noob, does anyone have any ideas of what else I could try?

I am using Pymol 2.3.1 and Python 3.7.1

# Answer
Thanks a lot to all who commented/replied, the author of the scripts provided a response. It was down to a syntax change in PyMol of which I wasn't aware that caused my selection to come up empty.
> 1 votes

---
Tags: python, protein-structure, protein-protein-interaction, structural-biology, pymol
---


# Question
Title: Why does the Genbank phylogenetic tree of 4000 sars-cov-2 sequences only display 200 nodes?
I manually selected all 4000 sequences on the 20 pages of the SARS-CoV-2 sequences on genbank here.

When I click on `Build phylogenetic tree`, it only has 200 nodes.

Is this expected? I was expecting to see 4000 nodes.

Is it simply a displaying issue where only 200 sequences per page are revealed while browsing the list of sequences? Is there a way to get the tree for all 4000 sequences?

# Answer
Okay, I'll provide a formal answer, Genbank phylogeny is to get a gist of the diversity. It isn't a formal phylogenetic analysis. You can down the tree as a pdf, I don't think you can download the treefile (in fact I'm pretty certain you can't). Genbank will perform a reasonable distance method and from recollection it will "collapse clades", so if you see little triangles on the tree it represents a large number of taxa represented as a triangle.

I was surprised to learn there are 4000 SARS-CoV-2 sequences.

To make a tree you need to download the data, align (muscle) and send through a formal phylogeny program. If you think Genbank is okay, then try MegaX its very point and clicky. Generally, maximum likelihood phylogeny has returned to vogue for at least 1 decade now (previously Bayesian). However, the Beast package is a Bayesian analysis which is still going strong and the main tool now used in molecular dating.

In summary, to perform a tree analysis beyond Genbank

* use the download option (there are two variants - the entire genome or the cut locus based on the query sequence)
* align using muscle or clustalo (muscle is really quick(
* perform a tree analysis using MEGA X, phylml, raxml or Beast depending on what you objectives are
> 2 votes

---
Tags: phylogeny, covid-19, sars-cov-2, genbank
---


# Question
Title: PCA plot did not work in single cell RNA-seq
I run plotPCA for single cell RNA-seq data, while I get error message (I use R 4.0). I attached the code and error message here. Did anyone know the reason and how to deal with it. Thanks!

```
library(pcaMethods)
library(SC3)
library(scater)
library(SingleCellExperiment)
library(pheatmap)
library(mclust)
set.seed(1234567)
deng <- readRDS("deng-reads.rds")
plotPCA(deng, colour_by = "cell_type2")

```

```
Error in value[[3L]](cond) : invalid subscript 'type' in 'reducedDimNames(<SingleCellExperiment>)'
'PCA' not in 'reducedDimNames(<SingleCellExperiment>)'

```

I also get error message when I run the code as below:

```
colData(deng)$tSNE_kmeans <- as.character(kmeans(deng@reducedDims$TSNE, centers = 8)$clust)

```

```
# Error in array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x), : 'data' must be of a vector type, was 'NULL'

```

# Answer
> 1 votes
The error is pretty obvious. There is no PCA data in the respective slot of the SingleCellExperiment that `plotPCA` could use, at least it is not called `PCA` which by default is expected. Use `runPCA` first, please read the manuals of `scater` and `scran` first, they have step-by-step code.

---
Tags: r, scrnaseq, bioconductor
---


title: Very low probability densities

text:
# Question
Title: Very low probability densities
I'm generating a density plot using the ggplot2 library:

```
ggplot(df_chr20, aes(x=start1)) + geom_density() + theme_bw()
ggplot(df_chr21, aes(x=start1)) + geom_density() + theme_bw()

```

The data frame for *df\_chr20* has 1573 rows, and I wanted to plot the probability density of genomic regions as x = *start1*. The same for *df\_chr21*, with 1271 rows. What I see in the plot is what I expect, but is there a reason why the y-axis probablilty densities are so low? The area under the curve is surely not equal to 1.

I expected the high localised high density in chr 21, and both the plots have similar scales in terms of y-axis values. But since they are so low, how can I interpret this plot?

I found some information here, but my values are much lower than the examples.

Any help would be very kind!

# Answer
It is not like your conventional frequency which adds up to one. The density is low because the width of your bins is huge, and the number of observations you have is low. From this by Wickham, the basic kernel is:

where K is the kernel and h is the bandwidth (Scott, 1992b)

If you don't specify it, by default h will revert to nrd0:

```
library(rtracklayer)
gr_narrowPeak <- import("https://www.encodeproject.org/files/ENCFF072FBH/@@download/ENCFF072FBH.bed.gz",format = "BED",extraCols = extraCols_narrowPeak)
gr_chr21 = gr_narrowPeak[seqnames(gr_narrowPeak)=="chr21"]

da = as.data.frame(gr_chr21[rep(1:length(gr_chr21),gr_chr21$signalValue)])
dens = density(da$start)
dens$bw
[1] 1133889

plot(dens)

```

```
ggplot(da,aes(x=start))+geom_density()

```

You can decrease h such that you get something smaller (but at the cost of looking coarse). Important thing is it doesn't sum up to one.

If you would like to explain it easier like a histogram, consider converting it into counts, like explained in this post:

```
binW = 50000

ggplot(da,aes(x=start))+
geom_density(aes(y = ..density..*(nrow(da)*binW)),bw=50000) +
geom_histogram(binwidth=binW,fill="steelblue",alpha=0.5) +
theme_bw()

```

```
ggplot(da,aes(x=start))+
geom_histogram(aes(y=..count../nrow(da)),binwidth=50000) +
geom_density(aes(y = ..density..*binW),bw=binW,fill="steelblue",alpha=0.5) +
theme_bw()

```
> 3 votes

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


# Question
Title: Why use Needleman-Wunsch if there is no way to evaluate the statistical significance
I thought that Needleman-Wunsch is the best approach to align sequences. However, I read that it is impossible to evaluate the statistical significance of the alignment if you do global alignment. So I would like to ask, what is the point of global alignment if we cannot evaluate statistical significance?

# Answer
> 0 votes
It's best to flip the question on its head and ask yourself instead, "Do we care about the statistical significance of an individual alignment?" The answer is "yes" if you're in a context like homology modeling, where you're aligning two protein sequences. The answer, however, is a very definite "no" in the context of NGS, where you just want the most likely alignment and a vague assessment of its likelihood of being wrong (for filtering purposes).

---
Tags: sequence-alignment, algorithms
---


title: Edit FASTA header using sed

text:
# Question
Title: Edit FASTA header using sed
I need to rename the following headers from a FASTA file. Something like:

```
>ENST00000488147.1::chr1:29570-30013(-)
>ENST00000405392.5::chr2:24919186-24920186(-)

```

into different format like this:

```
>chr1:29570-30013(-)_ENST00000488147.1
>chr2:24919186-24920186(-)_ENST00000405392.5

```

How can I rename the FASTA headers, perhaps using `sed` or `awk`?

# Answer
> 6 votes
```
awk -F "::" '{if($1~">"){gsub(">","");print ">"$2"_"$1} else {print $0}}' foo.fa

```

Basically the same as from Devon but using `-F` to indicate the initial field delimiter and then using `gsub` to remove the `>`.

# Answer
> 8 votes
```
awk 'BEGIN{FS="::"}{if($1~">"){printf(">%s_%s\n",$2,substr($1,2))}else{print $0}}' input_file.fa > output.fa

```

An explanation of the non-obvious bits:

* `BEGIN{FS="::"}` split columns using :: as the designator
* `$1~">"` If there's a \> in the line
* `substr($1,2)` Trim the `>` that was at the beginning of the line.

# Answer
> 3 votes
`sed -E 's/>(.+)::(.+)/>\2_\1/' in.fa > out.fa`

If you want to use sed.

# Answer
> 2 votes
I would use `perl` here instead of `sed` so you can use non-greedy patterns (e.g. `.*?`) and so ensure that you always match the *first* occurrence of `::` if there are more than one on the line. Perl also has `-i`, and in fact is where sed got the idea from, so you can edit the file in place just like you can with `sed`.

Using this example file:

```
$ cat file.fa
>ENST00000488147.1::chr1:29570-30013(-)
ACCTCTGAGAGCACCTCTGAGAGCACCTCTGAGAGCACCTCTGAGAGC
ATATATCGGCATATATCGGCATATATCGGCATATATCGGCATATATCG
>ENST00000405392.5::chr2:24919186-24920186(-)
ACTCTAGGATATCGCGATATCGCCTTAGGCTAGAATCGGCTTAGC
TCGCGATATCATCGCCTTAGGCTAGAATCGTCTAGGATATCGCG
CCTTAGGCTATCATCGCCTAGAATCGTCTAGGATATCGAATC

```

You can do:

```
$ perl -pe 's/>(.+?)::(.*)/>$2_$1/' file.fa
>chr1:29570-30013(-)_ENST00000488147.1
ACCTCTGAGAGCACCTCTGAGAGCACCTCTGAGAGCACCTCTGAGAGC
ATATATCGGCATATATCGGCATATATCGGCATATATCGGCATATATCG
>chr2:24919186-24920186(-)_ENST00000405392.5
ACTCTAGGATATCGCGATATCGCCTTAGGCTAGAATCGGCTTAGC
TCGCGATATCATCGCCTTAGGCTAGAATCGTCTAGGATATCGCG
CCTTAGGCTATCATCGCCTAGAATCGTCTAGGATATCGAATC

```

Or, to edit the file in place:

```
$ perl -i -pe 's/>(.+?)::(.*)/>$2_$1/' file.fa
$ cat file.fa
>chr1:29570-30013(-)_ENST00000488147.1
ACCTCTGAGAGCACCTCTGAGAGCACCTCTGAGAGCACCTCTGAGAGC
ATATATCGGCATATATCGGCATATATCGGCATATATCGGCATATATCG
>chr2:24919186-24920186(-)_ENST00000405392.5
ACTCTAGGATATCGCGATATCGCCTTAGGCTAGAATCGGCTTAGC
TCGCGATATCATCGCCTTAGGCTAGAATCGTCTAGGATATCGCG
CCTTAGGCTATCATCGCCTAGAATCGTCTAGGATATCGAATC

```

### Explanation

* `perl -pe`: **p**rint each line of the input file after applying the script given by `-e` to it.
* `s/>(.+?)::(.*)/>$2_$1/'`: the `s/old/new/` is the substitution operator and will replace `old` with `new`. Here, we look for the shortest stretch of one or more characters (`.+?`) from the `>` until the first occurrence of `::`, and capture it as `$1` (because the pattern is in parentheses). Everything after the first `::` is then captured as `$2` (`(.*)`). Then, we replace everything we matched with a `>`, then the 2nd captured pattern, an underscore and the first captured pattern.

---

Just for fun, here's a golfed awk approach:

```
$ awk -F'::' '$0=sub(/>/,"")?">"$2"_"$1:$0' file.fa
>chr1:29570-30013(-)_ENST00000488147.1
ACCTCTGAGAGCACCTCTGAGAGCACCTCTGAGAGCACCTCTGAGAGC
ATATATCGGCATATATCGGCATATATCGGCATATATCGGCATATATCG
>chr2:24919186-24920186(-)_ENST00000405392.5
ACTCTAGGATATCGCGATATCGCCTTAGGCTAGAATCGGCTTAGC
TCGCGATATCATCGCCTTAGGCTAGAATCGTCTAGGATATCGCG
CCTTAGGCTATCATCGCCTAGAATCGTCTAGGATATCGAATC

```

---
Tags: fasta, bash, linux
---


# Question
Title: Converting hmmer --tblout to pandas dataframe
Is there a way to convert a hmmer output to a pandas dataframe? I am also unsure how to load a hmmer tblout table into python via the Bio module. I believe you can call a hmmer format with SeqIO.parse or SeqIO.search. How can you convert the table to a pandas dataframe?

# Answer
The manual says fields in a tabular output file are space-delimited. It should be fairly straightforward to load a space-delimited file in pandas.

From the manual:

> `--tblout <f>` Save a simple tabular (space-delimited) file summarizing the per-target output, with one data line per homologous target model found.
> 2 votes

---
Tags: hmmer, pandas, dataframe
---


# Question
Title: snakemake - using output file not explicitly declared as input for subsequent rule
I am writing a snakemake that, up to now should

1. use raw sequencing reads to produce a draft assembly using canu
2. use minimap2 to align the raw reads against the draft assembly, generated by canu.

The following two rules are described below:

```
rule canu:
    input:
        fqdir + "{sample}.fastq"
    output:
        directory("{sample}/canu")
    shell:
        "/home/bhinckel/downloads/canu-1.9/Linux-amd64/bin/canu genomeSize=4.1m -p {wildcards.sample} -d {output} -nanopore-raw {input}"

rule minimap:
    input:
        reads = fqdir + "{sample}.fastq",
        draft_assembly = "{sample}/canu/{sample}.contigs.fasta"
    output:
        "{sample}/aligned/{sample}.sam"
    conda:
        "envs/general.yaml"
    shell:
        "minimap2 -ax map-ont {input.draft_assembly} {input.reads} > {output}"

```

However, if I do `snakemake -np` I get `Missing input files for rule minimap`. Which does make sense, after all I am not specifying the input for rule `minimap` (`{sample}/canu/{sample}.contigs.fasta`) anywhere i.e. snakemake cannot know that the file `{sample}/canu/{sample}.contigs.fasta` will be generated by canu.

The thing is that `canu` produces several files as output (amongst which the file `{sample}/canu/{sample}.contigs.fasta` is one of them), that's why I specified its output as a `directory`.

One way of solving this I guess would be to specify all output files produced by canu as an individual rule output (e.g. `{sample}/canu/{sample}.contigs.fasta`, `{sample}/canu/{sample}.report`, `{sample}/canu/{sample}.correctedReads.fasta.gz`), though I am looking for a more robust way to do so (if that's possible at all).

# Answer
Another way is to add one command after minimap: `touch done.txt`. This creates an empty file after minimap is finished. This can then be used in output rule for snakemake.

But specifying all the output files may not be a bad idea if these outputs are used in subsequent rules.
> 1 votes

---
Tags: ngs, nanopore, snakemake
---


# Question
Title: Error: Insufficient values in manual scale. 24 needed but only 1 provided?
folks.

Can you give me idea about an issue that I encountered? I ran single-cell RNA sequencing data using Suerat in R, and when I tried to draw violin plot, there is an error as below.

> VlnPlot(tumors, "Cd79a", "Foxp3", "CD14")

Error: Insufficient values in manual scale. 24 needed but only 1 provided.

Run `rlang::last_error()` to see where the error occurred.

However, I can draw the violin plot for Cd79a. The problem happened to only several genes. I haven't solved this problem yet. Do you know what should we do?

Thank you.

# Answer
You must supply genes (or features) as a vector, not individually:

```
VlnPlot(tumors, c("Cd79a", "Foxp3", "CD14"))

```

By the way it is probably `Cd14`, not `CD14`.
> 1 votes

---
Tags: scrnaseq, seurat
---


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

Much Appreciated

# Answer
Check out the tool `Intervene`. It's exactly designed for this. Github link.
> 1 votes

# Answer
Venn diagrams can be helpful, but for peak calling they can be somewhat misleading. This topic is addressed in the DiffBind documentation, in the "Comparison of occupancy and affinity based analyses" section.

Suffice to say, asking "which peaks are shared?" is not the same as asking "is there evidence for differential peaks, and where are they?". For this reason, performing a differential binding analysis might be a better option.

If you're just looking for Venn diagram plotting software, though, the VennDiagram (note, I have personally contributed to this tool in the past) and UpSet R packages are useful tools.
> 1 votes

---
Tags: sequence-alignment, sequence-analysis, chip-seq, peak-calling, macs2
---


title: Comparing phylogeny in R

text:
# Question
Title: Comparing phylogeny in R
So I want to compare the phylogeny created using two methods for example Maximum likelihood and maximum parsimony.Is there any way to compare the two phylogeny ?

I did read about phangorn but not sure if its the right R library for comparative analysis.

Any suggestion or helped would be highly appreciated

My data file

```
library(phangorn)
library(phytools)
library(dendextend)

data <- read.dna("abhi_seq/clean_dup_align_fast.fas", format = "fasta")
data
dat <- as.phyDat(data)

dm <- dist.ml(dat)
treeUPGMA <- upgma(dm)
treeNJ <- NJ(dm)

layout(matrix(c(1,2), 2, 1), height=c(1,2))
par(mar = c(0,0,2,0)+ 0.1)
plot(treeUPGMA, main="UPGMA")
plot(treeNJ, "phylogram", main="NJ")

dev.off()

parsimony(treeUPGMA, dat)
parsimony(treeNJ,dat)

tr.mp <- optim.parsimony(treeNJ, dat)

#tr.ml = optim.pml(treeNJ, dat)

fit <- pml(treeNJ, dat)
fit <- optim.pml(fit, rearrangement="NNI")

fit.ini <- pml(treeNJ, dat)
fit.ini

fit <- optim.pml(fit.ini, optNni=TRUE, optBf=TRUE, optQ=TRUE, optGamma=TRUE)
fit

tr.ml <- root(fit$tree,1)

tr.mp.ultra<-force.ultrametric(tr.mp)
tr.ml.ultra<-force.ultrametric(tr.ml)

is.ultrametric(tr.mp.ultra)
is.binary.tree(tr.mp.ultra)
is.rooted(tr.mp.ultra)

dd.ml.ultra<-as.dendrogram(tr.mp.ultra)

```

> Error in ape::as.hclust.phylo(object) : the tree is not rooted

I ran into this error this error which says trees is not rooted

**Comparing phylogeny code updated \** *Working* thanks to \***thomas duge de bernonville\*** thomas for putting codes together and fixing the errors

```
library(dendextend)
library(seqinr)
library(phytools)
library(phangorn)

a<-read.alignment("abhi_seq/clean_dup_align_fast.fas", format="fasta")
a.phydat<-as.phyDat(a)
dist.a.phydat<-dist.dna(as.DNAbin(a.phydat))
upgma.a<-upgma(dist.a.phydat)
parsimony(upgma.a,a.phydat)
pars.a <- optim.parsimony(upgma.a, a.phydat)
pars.a<-acctran(pars.a, a.phydat)
pars.a.rooted<-root(pars.a, outgroup="AAA64460", resolve.root=T)
pars.a.rooted.dd<-as.dendrogram(force.ultrametric(pars.a.rooted))

mt <- modelTest(a.phydat, tree=upgma.a,multicore = TRUE,mc.cores=10)

#ml.a = pml(upgma.a,a.phydat)
#fitJC <- optim.pml(ml.a, TRUE)
#ml.a.rooted<-root(midpoint(fitJC$tree), outgroup="AAA64460", resolve.root=T)
#ml.a.rooted.dd<-as.dendrogram(force.ultrametric(ml.a.rooted))

#########################################################################

######################################################################33

#mt <- modelTest(dat, tree=tree, multicore=TRUE)
mt[order(mt$AICc),]# choose best model from the table according to AICc
bestmodel <- mt$Model[which.min(mt$AICc)]
env = attr(mt, "env")
fitStart = eval(get(bestmodel, env), env)
fit = optim.pml(fitStart, rearrangement = "stochastic",optGamma=TRUE, optInv=TRUE, model="GTR")#tree
bs=bootstrap.pml(fit, bs=25, optNni=TRUE, multicore=TRUE)

ml.a.rooted<-root(midpoint(fit$tree), outgroup="AAA64460", resolve.root=T)
ml.a.rooted.dd<-as.dendrogram(force.ultrametric(ml.a.rooted))

```

# Answer
> 3 votes
phangorn is a really powerful package for phylogenies. But to compare the trees, I think you may convert them into dendrograms and calculate a correlation measure, such as the Fowlkes-Mallows Index or a distance measure such as Baker’s Gamma Index. These can be easily computed using the `dendextend` R package (https://cran.r-project.org/web/packages/dendextend/vignettes/dendextend.html#correlation-measures).

# Answer
> 4 votes
Another 2 tree compare tools: ggtree and blatic 3

# Answer
> 2 votes
**Bootstrapping** I agree correlation is one way in, but the classic way to do this is via bootstrapping, i.e. resampling the alignment with replacement and making a consensus phylogeny of of 100 to 1000 replicates. Values exceeding 75 to 80% define node robustness and you can calculate incongruence from there.

`Phangorn` certainly does bootstrapping, but the efficiency of the algorhithm for maximum likelihood (ML) I don't know (might take a very long time), just think how long one tree takes and times it by 100 (or even 1000).

If you just performed a correlation measure there is a risk the differences are not robust, but if you've got loads of taxa you could bootstrap collapse any value less than the robustness threshold into a polytomy and then perform a correlation method. Its a bit complicated, but it would work.

Bootstrapping parsimony isn't hard nor computationally expensive, bootstrapping maximum likelihood sometimes requires specific efficient algorhithms notably `RAxML`.

**Non-R** The other way in is just dump the trees into `Dendroscope3`, okay its not R but will work. Dendroscrope is hard for a non-tree person to read, however it will flag the differences for a specialist. Dendroscope will make a network tree for the topological differences, in this case between methods, i.e. non-bifurcating. You really then usually have to present the differences as a "mirror tree", i.e. parsimony and ML side by side, because usually non-tree people wonder what a Dendroscope output means. Technically there is nothing wrong with describing incongruence via networks.

Basically, the heart-beat is bootstrapping.

Generally speaking parsimony will give the same basic output as ML, except when you have some taxa undergoing rapid evolution against their sister group. This does happen in pathogen evolution and parsimony falls over into a phenomina known as 'long-branch attraction'.

# Answer
> 2 votes
If you only have to compare two phylogenies (as opposed to more), I would suggest that an alternative to using R would be using iqtree . One of its parameters (-rf) can calculate the Robinson-Foulds metric between two trees. (you can also use iqtree to calculate the phylogenies in the first place)

It is command-line based, but can be easily installed with conda. Having anaconda installed in general makes bioinformatics work easier, I feel.

Alternatively, differences in trees can be visualized (without numbers) in Dendroscope, or Splitstree. They have different strategies for it. With Dendroscope it is easier to compare side-by-side and, if having precisely the same taxon names, to draw lines between taxa in different trees. Splitstree can accept an input of multiple trees (File\>Tools\>Load Multiple trees) and build a network (which can be a way to visualize the consensus of the trees).

---
Tags: r, phylogenetics, phylogeny, phangorn
---


# Question
Title: Adding external annotation to the phylogeny tree
I was looking to overlay annotation on my phylogeny made,which was solved here question asked where i was looking to compare two different phylogeny which i could do.

Now i want to label the tips into various groups such as if its human or canine etc.One of the possible solution i found is this one not sure how to implement it.

alignment file annotation file

```
library(dendextend)
library(seqinr)
library(phytools)
library(phangorn)

a<-read.alignment("clean_dup_align_fast.fas", format="fasta")
a.phydat<-as.phyDat(a)
dist.a.phydat<-dist.dna(as.DNAbin(a.phydat))
upgma.a<-upgma(dist.a.phydat)
parsimony(upgma.a,a.phydat)
pars.a <- optim.parsimony(upgma.a, a.phydat)
pars.a<-acctran(pars.a, a.phydat)
pars.a.rooted<-root(pars.a, outgroup="AAA64460", resolve.root=T)
pars.a.rooted.dd<-as.dendrogram(force.ultrametric(pars.a.rooted))

```

i want to label according to the start column which is my source in my annotation file Any suggestion or help would be really appreciated

My working solution using this , this works as I did the alignment of sequences using full header name which kind of help to look for the pattern and grep it. The data file

I would still like to see the ggtree way instead of putting lots of label in the plot I would only like to see if I can annotate all the tips which have one organism source as one color rather putting all the labels

```
pars.a
# "Phylogenetic tree with 516 tips and 514 internal nodes"

# by using the class() function
class(pars.a)
# "phylo"

# or by using the str() structure function
str(tree)
str(pars.a)
# "List of 4"
# this list includes $edge, $Nnode, $ tip.label and $edge.length
# the tree$tip.label includes family designation
#tree$tip.label  # 516 of these
pars.a$tip.label
# from the Science paper, we have seven kinase families:
# kinase categories... TK, TKL, STE, CK1, AGC, CAMK, CMGC
# with the following colours
# "red", "green", "paleblue", "orange", "yellow", "purple", "pink", "green"

# by using the grep()function on the tree$tip.label part of the object
# we can find the tip labels that include "TK/" - i.e. tyrosine kinases
#grep("TK/", tree$tip.label)  # gives a list of numbers with "TK/" in tip label
#length(grep("TK/", tree$tip.label))

grep("Human",pars.a$tip.label)
length(grep("Human",ml.a.rooted$tip.label)
)
# thus there are 94 tip labels with that are designated TK (not TKL tyrosine kinase like)

# make a vector for each tip.label called tipcol with black on all of these...
#tipcol <- rep('black', length(tree$tip.label))
tipcol <- rep('black',length(pars.a$tip.label))

# make a vector with our list of kinase categories
#kinaseCats <- c("TK/", "TKL", "STE", "CK1", "AGC", "CAMK", "CMGC", "RGC")
kinaseCats <- c("Human","Rat","Murine","Bovine","Canine","Feline","Porcine","Gallus","Bat","Infectious","Transmissible")

# make a vector of color we want:
colorsList <-c("red", "darkolivegreen3", "blue", "orange", "blueviolet", "purple", "pink", "green","brown","navyblue","purple4")

# replace colours where grep gives "TK" as red, etc in a loop
#for(i in 1:length(kinaseCats)){
 # tipcol[grep(kinaseCats[i], tree$tip.label)] <- colorsList[i]
#}

for(i in 1:length(kinaseCats)){
  tipcol[grep(kinaseCats[i], pars.a$tip.label)] <- colorsList[i]
}

# plot with edge length false to see nodes better
#plot(tree,
 #    use.edge.length = FALSE,
  #   tip.color=tipcol,
   #  cex = 0.25)

plot(pars.a,
     use.edge.length = FALSE,
     tip.color=tipcol,
     cex = 0.5)
# slow to draw due to text - a bit annoying!
nodelabels(cex=0.4)

```

figure

# Answer
> 1 votes
there are nice features in R **ggtree** package. You may use the *geom\_cladelabel()* function for your purposes. But I haven't find yet a way to automatize clade annotation. To do it, you have to be sure of the clade homogeneity before removing tiplabs. here are some commands to do clade label annotation:

```
annot<-read.delim("annota.txt", header = T)

ggtr<-as.treedata(pars.a)
g<-ggtree(ggtr)
#find the node number corresponding to the common ancestor of the top clade in you tree
MRCA(ggtr, "AAL40400", "AAF25519")
#here, check that all accessions belong to the same organism/group
subset(annot, label=="AAL40400" | label=="AAF25519")
#in this case, sequences are related to bovine.
g<-g + geom_cladelabel(node=89, label="Bovine", align=F, fontsize=4)
plot(g)

#another example
MRCA(ggtr, "AAA46456", "AAF97738")
subset(annot, label=="AAA46456" | label=="AAF97738")
g<-g+ geom_cladelabel(node=73, label="Rodent", align=F, fontsize=4)
plot(g)

```

---
Tags: r, phylogeny
---


# Question
Title: How can I plot multiple sequence alignment in R as a heatmap
I have many amino acids sequences in a fasta format that I did multiple sequence alignment. I was trying to plot something like binary a code as heatmap. If it had a change in the AA compared to the reference sequence it would be red, if it did not change would be yellow, for example.

I came across to `msaplot` from `ggtree` package. I also checked `ggmsa` package for that. But so far, I did not get what I wanted. So basically I wanted to:

1. change the multiple sequence alignment to a binary matrix (if the amino differs from the reference sequence, plot x, if not y)
2. plot as a heatmap with two colors based if the amino acid changed or not

here's some data example of what I doing, just instal `ggmsa` and you can try this:

```
protein_sequences <- system.file("extdata", "sample.fasta", package = "ggmsa")
ggmsa(protein_sequences, start = 265, end = 300, color = "Chemistry_AA",font = NULL)

```

# Answer
@StupidWolf answered me on stackoverflow and he's answer solved my question. I will share it here as well.

@StupidWolf:

We read in the alignment:

```
library(Biostrings)
library(ggmsa)
protein_sequences <- system.file("extdata", "sample.fasta", package = "ggmsa")
aln = readAAMultipleAlignment(protein_sequences)
ggmsa(protein_sequences, start = 265, end = 300)

```

Set the reference as the 1st sequence, some Rattus, you can also use the consensus with consensusString() :

```
aln = unmasked(aln)
names(aln)[1]
[1] "PH4H_Rattus_norvegicus"

ref = aln[1]

```

Here we iterate through the sequence and make the binary for where the sequences are the same as the reference:

```
bm = sapply(1:length(aln),function(i){
as.numeric(as.matrix(aln[i])==as.matrix(ref))
})

bm = t(bm)
rownames(bm) = names(aln)

```

The plot you see above has sequences reversed, so to visualize the same thing we reverse it, and also subset on 265 - 300:

```
library(pheatmap)
pheatmap(bm[nrow(bm):1,265:300],cluster_rows=FALSE,cluster_cols=FALSE)

```

The last row, is Rattus, the reference, and anything similar to that is read, as you can see in the alignment above, last 4 sequences are identical.
> 5 votes

---
Tags: r, sequence-alignment, heatmap
---


title: How to get a MSA fasta from BAM/SAM?

text:
# Question
Title: How to get a MSA fasta from BAM/SAM?
I am working with paired end NGS miseq data from a viral genome with multiple timepoints. I have filtered and trimmed this data for quality and adapter sequences. I have then merged the filtered paired end reads and mapped them to a reference genome using bowtie2. I then used samtools to filter for mapping quality and visualized the resulting filtered bam file on IGV.

My goal is to now scan along the genome using various window sizes (e.g. 300 nt) across the reads that have mapped to the reference. I will then chop up the genome and build phylogenetic trees based on each window.

How can I extract a multiple sequence alignment of these 300 nt windows from a sam or bam file using the coordinates of a reference viral genome? I have tried using `samtools view` but this does not actually provide the subsection or alignment within the window but rather all reads that overlap the window in some way.

The virus is HIV1 and the immediate goal is to get an alignment of all reads that completely span a user specified coordinate range eg. 1-300 based on some reference sequence. Ideally reads that fit this criteria would be extracted and aligned. After ragged edges are removed, the resulting alignments subjected to phylogenetic analyses.

# Answer
> 2 votes
If need to get msa between every read and ref in the specific region , you can try this https://github.com/orangeSi/bam2msa.git

```
bam2sam ref.fa input.bam --regions chr1:200-300

```

---
Tags: sam, samtools, covid-19, sars-cov-2, bowtie2
---


title: separate genes with expression values

text:
# Question
Title: separate genes with expression values
I have a file like this with 2 columns. I want to separate the columns based on plus and minus signs into separate files. How do I do this with awk or sed. Kindly help.

# Answer
> 3 votes
What do you want to do with NA values?

The below scripts ignore them (h3ab74's answer will print them to positiveExpression)

```
awk '$2 ~ /^[0-9]+$/  {print $0}' expressionFile.txt > positiveExpression.txt
awk '$2 ~ /^-[0-9]+$/ {print $0}' expressionFile.txt > negativeExpression.txt

```

These do a regex match to respectively, positive and negative numbers. If there's anything else in field 2 the line won't be printed.

# Answer
> 1 votes
```
awk '{if($2>0) print $0}' expressionFile.txt >> positiveExpression.txt
awk '{if($2<0) print $0}' expressionFile.txt >> negativeExpression.txt

```

Assuming that the file with gene expressions has the title: expressionFile.txt, this should do the trick.

Explanation of first line of code: if values in the 2nd column - `$2` \- are greater than 0, print all columns (`$0`) and write to a file (`>>`) titled, positiveExpression.txt. The use of `awk '{}'` is just syntactical.

---
Tags: text-processing, awk
---


# Question
Title: How to calculate mutation rate and mutation sites in a genome using FASTA file?
I have 6 viral genome sequences of the same virus and 1 reference sequence in FASTA format.

1. How I can identify mutations and mutation sites in those genomes using FASTA sequences but how I can do this for FASTA file?
2. And how I can identify the mutation rate for one genome?

**Background** If I've FASTQ file then I'll simply align the reads to the reference and by using variant calling tool I will get the mutate sites

# Answer
> 1 votes
You can still align a FASTA file with a tool like bwa mem ( if they are short reads ) or minimap2 for long reads and run it through a variant caller like freebayes. Alternatively, if you have sufficient reads to create an assembly of all of the genomes, or if they already are you can create a MSA of them using a tool like Cactus (as it is built for genomes) or MAFFT if the sequences are short enough and then use your favorite viewer to find variable sites

---
Tags: sequence-alignment, covid-19, mutations, sars-cov-2, virus
---


# Question
Title: How to transform and sort the matrix to make a heatmap showing signatures?
I have a matrix `data` with cells as rows and samples as columns. Here I am giving the data with `dput`

```
dput(data)

structure(c(0, 0, 0.0095, 0.0091, 0, 0.0195, 0, 0.006, 0.0023,
0.0035, 0, 0.0306, 0, 0.0859, 0, 0, 0.0059, 0.0229, 0, 5e-04,
0, 0.0116, 0.0173, 0.0076, 0.0133, 0.0344, 0, 0.1263, 0, 0, 0.0378,
0, 0, 0.0198, 0, 0.0268, 0, 0, 0.0044, 0.0607, 0, 0.0073, 0,
0.1343, 0.0016, 0.0864, 0, 0, 0, 0.0296, 0.0122, 0.0329, 0.0204,
0.0582, 0, 0, 0.003, 0.1011, 0, 0, 0, 0, 0.0124, 0.0013, 0, 0.0025,
0.0076, 0.0498, 0, 0.0103, 2e-04, 0.0562, 0.0673, 0.069, 0, 0.0297,
0, 0.0207, 0.0329, 0.0072, 0.0028, 0.0589, 0.0533, 0.0139, 0,
0.0791, 0.0361, 0.026, 0, 0.0036, 0, 0.0319, 0.0155, 0.0401,
0.064, 0.0743, 0.0026, 0.012, 0, 0.0089, 0.0468, 0.0994, 0, 0.0646,
0, 0.0391, 0.0209, 0.0051, 0.007, 0.0689, 0, 0.0218, 0, 0, 0.0157,
0, 0, 0.0119, 0, 0.0136, 0, 0.0036, 0.0016, 0.0583, 0, 0, 0.0031,
0.2211, 0.6257, 0, 4e-04, 0.173, 0.0557, 0.1104, 0.0288, 0.0474,
0, 0.0078, 0.0914, 0.0289, 0.0293, 0.1067, 0, 0.0121, 0, 0, 0.018,
0, 0.021, 0.039, 0.0146, 0.043, 0.0367, 0.0245, 0, 0, 0, 0, 0,
0.0097, 0, 0.003, 0.0072, 0.017, 0.0061, 0.0453, 0, 0.0587, 0,
0, 0.1009, 0.0231, 0, 0.0279, 0, 0.048, 0.01, 0, 0.0083, 0.0355,
0, 0.1058, 0, 0, 0.0634, 0.0138, 0, 0.0116, 0, 0.0683, 0, 0,
0.0229, 0.064, 0, 0, 0, 0, 0.0222, 0.005, 0, 0.0176, 0, 0.0291,
0.0017, 0.008, 0.003, 0.0693, 0, 0.0156, 0.0038, 0.0577, 0, 0,
0, 0, 0.0521, 0, 0, 0.0017, 0, 0.0218, 0.1336, 0, 0, 0, 0, 0,
0, 0.0398, 0, 0.0099, 9e-04, 0.0046, 0.0055, 0.0361, 0.0382,
0.0168), .Dim = c(14L, 17L), .Dimnames = list(c("Cell1", "Cell2",
"Cell3", "Cell4", "Cell5", "Cell6", "Cell7", "Cell8", "Cell9",
"Cell10", "Cell11", "Cell12", "Cell13", "Cell14"), c("Sample1",
"Sample2", "Sample3", "Sample4", "Sample5", "Sample6", "Sample7",
"Sample8", "Sample9", "Sample10", "Sample11", "Sample12", "Sample13",
"Sample14", "Sample15", "Sample16", "Sample17")))

```

Annotation of samples is in a dataframe `df` look likes this:

```
dput(df)

structure(list(Samples = c("Sample2", "Sample6", "Sample7", "Sample8",
"Sample9", "Sample10", "Sample11", "Sample16", "Sample17", "Sample1",
"Sample3", "Sample4", "Sample5", "Sample12", "Sample13", "Sample14",
"Sample15"), Group = c("A", "A", "A", "A", "A", "A", "A", "A",
"A", "B", "B", "B", "B", "B", "B", "B", "B")), row.names = c(NA,
-17L), class = c("tbl_df", "tbl", "data.frame"))

```

Using the above data and annotation with below code I made a heatmap:

```
#Associated libraries
library(ComplexHeatmap)
library(RcolorBrewer)
library(circlize)

#Set annotation
ann <- data.frame(df$Group)
colnames(ann) <- c("Group")
colours <- list("Group"=c("A"="orange","B"="darkgreen"))
colAnn <- HeatmapAnnotation(df=ann, which="col", col=colours, annotation_width=unit(c(1, 4), "cm"), gap=unit(1, "mm"))

myCol <- colorRampPalette(c("#ffffff","#0030dd", "#3259e3", "#7f97ee"))(100)
myBreaks <- seq(0,1, length=100)
hmap <- Heatmap(data, name = "xCell scores", col = colorRamp2(myBreaks, myCol),
                width = unit(100, "mm"),show_row_names = TRUE, show_column_names = FALSE, cluster_rows = FALSE,
                cluster_columns = FALSE, show_column_dend = FALSE, show_row_dend = FALSE,
                row_dend_reorder = FALSE, column_dend_reorder = FALSE, clustering_method_rows = "ward.D2",
                top_annotation=colAnn)
draw(hmap,heatmap_legend_side="right", annotation_legend_side="right")

```

The heatmap looks like this:

But I want to sort the rows (cell1 to cell14) keeping the cell4 (interested cell) on the top and also keep the samples with values (coloured in blue in heatmap) together and samples with 0 values (white coloured) following that.

It should look something like the below heatmap which is from the Figure 4a in this research paper

As you see in this heatmap all the coloured samples are together.

This must be done with some sorting and using transform, but not able to make it.

Any help is appreciated. thank you in advance.

# Answer
> 1 votes
You are almost there except for a few tricks:

If you want columns to be ordered using the values in Cell 4, you will need to provide a custom column order either with the `column_order` argument of the `Heatmap()` function or by providing ordered data. Below is the code to perform the latter. Most of the code is yours, I just tweaked it a little.

```
library(tibble)
library(dplyr)
library(ComplexHeatmap)
library(RColorBrewer)
library(circlize)

# Since you want to use cell4 as reference
cell4 <- as.data.frame(t(data["Cell4",]))

df <- cbind(cell4,
            Group = as.character(ann$Group))

# ordering by Group and by expression value, minus denotes decreasing order
df <- df[with(df, order(Group,-Cell4)),]

# to make it compatible to your code in your post
df <- rownames_to_column(df, "Samples")

#Set annotation
ann <- data.frame(df$Group)
colnames(ann) <- c("Group")
colours <- list("Group" = c("A" = "orange", "B" = "darkgreen"))
colAnn <-
  HeatmapAnnotation(
    df = ann,
    which = "col",
    col = colours,
    annotation_width = unit(c(1, 4), "cm"),
    gap = unit(1, "mm")
  )

# rownames of the annotation df and the colnames of the data should be in
# the same order
data <- data[, df$Samples]

myCol <-
  colorRampPalette(c("#ffffff", "#0030dd", "#3259e3", "#7f97ee"))(100)
myBreaks <- seq(0, 1, length = 100)
hmap <-
  Heatmap(
    as.matrix(data),
    name = "xCell scores",
    col = colorRamp2(myBreaks, myCol),
    width = unit(100, "mm"),
    show_row_names = TRUE,
    show_column_names = FALSE,
    cluster_rows = TRUE, # you asked your cells to be ordered
    cluster_columns = FALSE, # column order is the same as the provided data
    show_column_dend = FALSE,
    show_row_dend = FALSE,
    row_dend_reorder = FALSE,
    column_dend_reorder = FALSE,
    clustering_method_rows = "ward.D2",
    top_annotation = colAnn,
    column_order = NULL
  )
draw(hmap,
     heatmap_legend_side = "right",
     annotation_legend_side = "right"

```

)

The resulting heatmap, still not quite what you want:

You will realize that the values for Cell 4 are ordered from high to low for each of the two groups. To make Cell 4 be drawn at the top, you will need to reorder the rows. To do so, you will need to extract the row orders with the `row_order()` function:

```
> row_order(hmap)
 [1]  3  2 14  4 12 11  9 10  7  1  5 13  6  8

```

Now you can re-draw the heatmap after altering the row order somehow. For this sample case, doing it manually.

```
# promoting Cell 4 to the first position
modified_row_order <- c(4, 3, 2, 14, 12, 11,  9, 10,  7,  1,  5, 13,  6,  8)

hmap <-
  Heatmap(
    as.matrix(data),
    name = "xCell scores",
    col = colorRamp2(myBreaks, myCol),
    width = unit(100, "mm"),
    show_row_names = TRUE,
    show_column_names = FALSE,
    cluster_rows = FALSE, # set to FALSE as a custom row order is supplied
    cluster_columns = FALSE,
    show_column_dend = FALSE,
    show_row_dend = FALSE,
    row_dend_reorder = FALSE,
    column_dend_reorder = FALSE,
    clustering_method_rows = "ward.D2",
    top_annotation = colAnn,
    column_order = NULL,
    row_order = modified_row_order # argument to supply custom row order
  )
draw(hmap,
     heatmap_legend_side = "right",
     annotation_legend_side = "right")

```

If I have misunderstood your question and you do not want columns to be ordered according to expression values in Cell4, forget about the first part with the ordered annotation data frame.

# Answer
> 0 votes
You may also use the excellent ***pheatmap*** package. It is highly customizable, so do not hesitate to look at its help page and examples found therein. The code may look like this: `library(pheatmap)
data<-as.data.frame(data)
df<-data.frame(df, row.names=1)
pheatmap(data, annotation_col = df, cluster_cols = F, color=colorRampPalette(c("white", "steelblue"))(100))`

---
Tags: r, rna-seq, heatmap, figure-reproduction
---


# Question
Title: What are the p01sum and p0prop parameters in godon?
When using branch-site model in Godon, model optimization results include the following parameters: `p01sum` and `p0prop`. What are those and how do they relate to `p0` and `p1` in the branch-site model?

# Answer
> 1 votes
**Disclaimer**: I'm the author of godon. The question above comes from the real communication with a software user.

In the godon implementation of the branch-site model I performed the following reparametrization:

$p\_{01sum} = p\_0 + p\_1$

$p\_{0prop} = \frac{p\_0}{p\_0 + p\_1}$.

This has the following advantages:

* Both parameters have a well defined range $(0, 1)$.
* There is less dependency between the two. E.g., $p\_2=1-p\_0-p\_1$ (proportion of sites under positivie selection) does not depend on $p\_{0prop}$.
* Both properties are very helpful for both likelihood maximization and MCMC.

In case you are interested in $p\_0$ and $p\_1$ there is a straightforward way to go back to the original parameters:

$p\_0=p\_{0prop}\*p\_{01sum}$

$p\_1=p\_{01sum}-p\_0$.

P.S. I also added this information to the tutorial.

---
Tags: codon, positive-selection, godon
---


# Question
Title: about the cpm normalization after using normalization factor
Is it okay to use CPM normalization (with/without log transform) after using TMM normalization? Why do we need both?

```
 library(edgeR)
 library(SummarizedExperiment)
 load(url("http://duffel.rail.bio/recount/SRP049355/rse_gene.Rdata"))
  counts <- assays(rse_gene)$counts
  y <- as.matrix((counts))
 y <- DGEList(counts = y, group=c(1,2,3,4,5,6,7,8,9,10))
  y <- calcNormFactors(y)
 z <- cpm(y, normalized.lib.size=TRUE)
scaledata <- t(scale(t(z))) # Centers and scales data.
hc <- hclust(as.dist(1-cor(scaledata, method="spearman")), method="complete") # Clusters columns by Spearman correlation.
 TreeC = as.dendrogram(hc, method="average")
 plot(TreeC,
 main = "Sample Clustering",
 ylab = "Height")

```

# Answer
> 2 votes
In fact the `cpm` function will use the size factors (TMM) that were calculated with `calcNormFactors`. In the absence of the size factors `cpm` would do a naive per-million scaling so only correct for differences in library size. The power of the size factors is to also correct for library composition. I recommend this video for details. In short: What you do is fine. I would suggest though that you use the `log = TRUE` option for `cpm` as one typically performs downstream applications (even when using the Z-score) on the log scale. Or alternatively `t(scale(t(log2(z+1))))`.

---
Tags: r, rna
---


# Question
Title: What is the full script to extract sequences from fasta file by using Ids in text file in python 3 and pycharm?
I have a fasta file (gene.fasta)with sequences with names patterns:

```
>Hsc_gene_9627.t1
ATGGCACGCATTTTCCTCATTCTTTTATTGCTGCACAACATCTGCTGTGCAGCCGCTTCA
TTGCTCATTTTGAATGCCGTTACATTGGAGAAGGATGCTAATGATTATGCCGTTGGCGAT

```

and I have the Ids in text file (seq.txt) which are not the exact as in the fasta file:

```
HSC_gene_996
HSC_gene_9734
and some of the names came as
HSC_gene_996|HSC_gene_9734

```

How can extract the sequences? I am not experienced in python so please use python for dummies language :) Thank you in advance

# Answer
For example, the fasta are like these:

```
gene.fasta:

>Hsc_gene_9627.t1
ATGGCACGCATTTTCCTCATTCTTTTATTGCTGCACA
ACATCTGCTGTGCAGCCGCTTCA
>Hsc_gene_9627.t2
GGGGTTTTCCCC
>Hsc_gene_962.t1
AGTCGTCAGTCAGTAGTCGC

seq.txt:

HSC_gene_9627
HSC_gene_9734
HSC_gene_996|HSC_gene_9734

```

We use a few modules, read in the fasta:

```
import re
from Bio import SeqIO
records = list(SeqIO.parse("gene.fasta","fasta"))

```

Read in the list:

```
seq = [i.rstrip().lower() for i in open('seq.txt').readlines()]

```

Here I use two functions, one to write fasta files, and the other to format the header for comparison:

```
def fasta_out(rec):
    return(">"+str(rec.description)+"\n"+str(rec.seq))

def format(rec):
    return re.sub("[.][^ ]*$","",rec.description).lower()

```

Define the records to keep:

```
keep = [fasta_out(i) for i in records if format(i) in seq]

```

Write them out:

```
f=open('matched.fasta','w')
f.writelines("\n".join(keep))
f.close()

```

We can check:

```
open('matched.fasta').readlines()
Out[50]:
['>Hsc_gene_9627.t1\n',
 'ATGGCACGCATTTTCCTCATTCTTTTATTGCTGCACAACATCTGCTGTGCAGCCGCTTCA\n',
 '>Hsc_gene_9627.t2\n',
 'GGGGTTTTCCCC']

```

Most likely there are some simple while, for loops that does not require modules..
> 1 votes

---
Tags: python, fasta
---


title: List of global ethnicities

text:
# Question
Title: List of global ethnicities
Looking at the 1000 genomes project I see the following list of population groups and super-groups: https://www.internationalgenome.org/category/population/

However, is there something like a standard list of global ethnicity? Like an ontology of population diversity or similar?

I'm trying to tie genome data to 'ancestry', and to do that I'd like to survey a set of people on their ancestral background (as far back as grandparents ideally).

How to capture diversity in terms of what people 'self report' in their backgrounds? e.g. middle / far east, north / south India, etc.

# Answer
It's not clear exactly what your goal is. It might help to have an idea of how you are trying to use this data; and it's probably not a bioinformatics question, as ethnicities are social/political rather than biological constructs.

However, there exist some controlled vocabularies used by LoC for ethnic groups; that looks a bit old so it might not be the best source.

I am by no means an anthropologist but it does look like there are anthropological databases devoted to e.g. ethnography that may have complementary information.

Wikipedia has a simple list of ethnicities that notes, instructively, "There has been constant debate over the classification of ethnic groups."

I think that you are unlikely to find a single authoritative answer that is safe to use as a source of truth. I would normally suggest to try another stackexchange site, though it looks like there is none for sociology or anthropology. There is no anthropology tag available on this site, but you might try posting in e.g. history.stackexchange.com and adding an anthro tag.
> 1 votes

---
Tags: human-genome, phylogenetics
---


title: Ranking metric in GSEA

text:
# Question
Title: Ranking metric in GSEA
I am using fgsea in r to calculate and plot a bunch of GSEA graphs.

My question is whether it is acceptable to use the Wald statistic from DeSeq2 to rank the gene list?

I have seen in the GSEA application that the signal to noise ratio is used but is the Wald statistic sufficient?

If not then is there an easyish way to calculate the signal to noise ratio in R given a normalised count matrix?

Thanks in advance.

# Answer
> 2 votes
I personally rank by signed nominal p-value. The fgsea author recommended on biostars.org to use either -log10(nominal p-value) or the F-statistics column (or whatever statistic the tool you use outputs) but not FDR/adjusted p-values as the latter produces a lot of ties for genes with low significances.

# Answer
> 3 votes
If using DESeq2 with GSEA, I'd recommend ranking by shrunk log2FC values.

It'd also be worth considering ranking positive and negative associations separately, because the standard GSEA algorithm doesn't cross at the zero point when associations change from positive to negative.

You shouldn't be using p-values to rank anything. A p-value gives information about the reliability of a result, but not its importance. See point #5 of this article about p-values:

https://amstat.tandfonline.com/doi/full/10.1080/00031305.2016.1154108#\_i31

---
Tags: r, deseq2, gsea
---


# Question
Title: Determine if an organism is microbic or unicellular
I am doing some metagenomic analysis of samples which can only contain microbes due to the experimental setup. For this, it would be useful to be able to distinguish whether a specific organism is a microbe or not, as it allows:

* to detect potentially false findings,
* to avoid an unnecessary inflation of the metagenomic database.

For animals and plants, this is straightforward, as they are never microbes. However, for fungi (and some other clades) it gets complicated. At the very least, I failed to find a simple list of microbial or non-microbial clades.

Ideally, I would like to have a function that maps a taxonomic ID to a boolean stating whether this organism is microbial. I acknowledge that this is utopic since *microbe* is not that well defined, there are facultatively multicellular organisms, and for some cases it may simply not be known. I am therefore happy with any reasonable approximation such as unicellularity. In particular, given my application, a small chance of falsely identifying an organism as microbial is not an issue.

I am currently working with the NCBI database accessed via the ETE 3 module for Python, but I am open to use other databases and tools.

# Answer
I would start by saying that your statement that no animals/plants are microbes is false. There are indeed microbial algae (plants) and microbial animals.

So it is indeed difficult. (I'll briefly note that unless your protocol specifically gets rid of environmental DNA, you will probably have some contaminating non-microbial DNA no matter what.)

However, we can use a different version of your proposition to make a simple (over-simple, most would argue) approximation: all prokaryotic organisms, e.g. bacteria and archaea, are microbial, and all eukaryotic organisms are not. This is not hard to do if you are already using NCBI taxonomy, as you can use e.g. this mapping. It has the advantage that all bacteria and archaea are *definitely* microbial.

Of course it's quite wrong, as you already know. In principle you could set up a complex decision tree sorting this out, though I imagine that's the thing you wish to avoid.

In terms of other genome databases, you might consider looking into GOLD, which has extensive metadata in a slightly less arcane setup than NCBI, which might be helpful. But there is a high rate of missingness there, and also there is still not necessarily a great "microbe"/"not microbe" mapping.

There are also existing databases of microorganisms, which may or may not satisfy your requirements; I have little familiarity with most of those so caveat emptor.

I believe (hope) that NCBI taxid is a good way to do lookups in all such databases, but it may not be true.
> 1 votes

---
Tags: public-databases, phylogenetics, taxonomy
---


# Question
Title: WGCNA co-expression network analysis with less than 20 samples
I have two cancer subtypes data. Subtype A is 14 samples and Subtype B is 23 samples. I'm interested in identifying the functions of some LncRNAs in the Subtype A group. For this I'm using all protein-coding genes and interested LncRNAs expression data (normalised counts).

From the WGCNA tutorial, I see that minimum number of samples needs to be 20. So, I guess using only Subtype A 14 samples for the analysis is not possible.

Can I use both Subtype A 14 samples and Subtype B 23 samples together for the analysis? In case if I use both, once I get the modules how do I know that the genes in modules are related to Subtype A or Subtype B? The pathways I get from the interested module is related to Subtype a or Subtype B? How this can be identified?

# Answer
It's hard to say without knowing how different the subtypes are. If you do a common WGCNA, you may find modules related to the differences between A and B as well as to drivers of expression variation in A and B (common and separate). Since A has fewer samples, the variation will have to be stronger to be reflected by WGCNA modules.

You can run the WGCNA, find modules that vary the most in A samples (these may or may not also vary in B samples), then see whether LncRNAs belong to these modules and have a high kME (module membership).
> 3 votes

---
Tags: r, rna-seq, bioconductor, gene-expression, wgcna
---


# Question
Title: Why is there a chloride ion in this 3D model?
Here we have a 3D model of a sequence that match with this sequence VVNIQKEIDRLNEVAKNLNESLIDLQELGKYEQYIK in Blast

the 3D structure is here

Why is there a chloride ion in this 3D model?

# Answer
> 3 votes
It is a crystallisation artefact. Namely, the protein are placed in a condition where they fall out of solution without aggregating a gloop (as happens in most well in xstal trial) and to pack as a crystal they need to be placed nice and orderly. It is coordinated by water so is not relevant. A lot of structures have these —DMSO and ions are the most commonly.

Two footnotes.

* 1ZV7 is a crystal structure, *i.e.* solved by X-ray diffraction. A model is generally used to described an *in silico* model, such as those from Swissmodel, phyre, I-Tasser etc. These are generated *ab initio* (from scratch) or threaded (mapped on top of) against crystal/EM/NMR structures. For Covid I-Tasser has a nice picture of its models.
* The best source of info about structures is the PDB. RCSB PDB is more commonly used, while PDBe is more tricker to navigate (or used to be). Both have the same data, but are better for the extras than NCBI.

---
Tags: protein-structure, covid-19, protein-protein-interaction, interactions, sars-cov-2
---


# Question
Title: How to tell if library generation for RNA-seq experiments are stranded vs not stranded
I would like to use RNA-seq data from the NIH GEO database. How can I tell if library generation for the RNA-seq experiments are stranded vs not stranded?

# Answer
> 3 votes
If you already have a reference transcriptome assembly, then you may run the Salmon tool for quantification using your RNA-seq samples. It has an automated library type inference function (https://salmon.readthedocs.io/en/latest/salmon.html#what-s-this-libtype). It might help you to determine whether the libs are stranded or not.

# Answer
> 2 votes
Unfortunately, the stranded nature of cDNA sequencing is not typically included in the metadata attached to a project, which means that downloading a small subset of the data and checking for strand consistency (as suggested by the other answer) is often the quickest way to work that out (when that information isn't included in the methods section of the associated paper). This knowledge is not particularly useful for casually browsing projects to find a stranded protocol that fits a particular study specification.

Sometimes you get lucky, but not often. GEO RNASeq datasets *usually* have links to the SRA accessions, and those accessions *usually* include annotation about the library type. For example, consider this dataset on GEO:

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE151535

It links to SRP265460 on SRA. Following the link for the first run for this project, I get to here:

https://www.ncbi.nlm.nih.gov/sra/SRX8436504\[accn\]

This includes the following metadata about the run:

> *Instrument:* NextSeq 500
> *Strategy:* RNA-Seq
> *Source:* TRANSCRIPTOMIC
> *Selection:* cDNA
> *Layout:* SINGLE
> *Construction protocol:* Replicates of oropharyngeal explants at different stages of development (37/38, 39, 41) were immersed in RNALater. Samples were homogenized with a Dounce homogenizer in RLT buffer containing ß-mercaptoethanol (Qiagen). The homogenate was then processed using an RNAEasy kit (Qiagen, Hilden, Germany). Replicate samples for each stage were analyzed for RNA quantity and quality on a Nanodrop 8000 spectrophotometer. RNA libraries were prepared for sequencing using standard Illumina protocols

But, as you can see, there's nothing about whether its stranded or not: Illumina has both stranded and unstranded "standard" protocols.

---
Tags: rna-seq, geo
---


# Question
Title: Compare nucleotide mutation with amino acid
To begin with I have two nucleotide sequence to compare I want to see how they might have evolved in terms of sequence similarity and differences. One of the way i can think of is align the sequences. But i m not getting how do i associate it with mutation.

Now question is how do i know about mutation. I have the respective amino acid data from which i know there been mutation which is making changes how the binding happens between two sequences.

So now i have to associate nucleotide mutation with amino acid sequence data.

Any suggestion or help would be really appreciated. As always if there is a `R` way to performing the analysis then it would be really good.

# Answer
> 1 votes
I would output the aligned nucleotide sequences and input them into an alignment editor such as `Sequtron` (Mac), `MEGAX`(particularly PC) or `Jalview` (Linux).

The alignment editors will allow you to flick between the triplet codon and the amino acid sequence. Therefore spotting which triplet codon mutations associate with each given amino acid mutation is easy.There is no `R` equivalent that I am aware of.

More generally I would do a non-synonymous / synonymous analysis which `fangorn` will perform. This assesses the number of mutations in which cause an amino acid change against the number of mutations that cause a 'silent' mutation (no amino acid change). A synonymous phylogeny will represent all mutations which don't change the tree. A non-synonymous phylogeny is the same as an amino tree.

There are character mapping approaches which will map an amino acid mutation onto a phylogeny as an alternative.

---
Tags: r, sequence-alignment, sequence-analysis
---


# Question
Title: What unit (TPM, FPKM/RPKM, or other) to use when working across samples
I have raw gene read counts and would like to perform an analysis across multiple samples. I've found conflicting info online on how this should be done. One commonality however is that FPKM/RPKM should not be used. On the other hand, there seems (to me at least) to be conflict about whether TPM is suitable. For example, this StatQuest video seems to suggest that TPM is suitable for analysis across samples. One person on this Biostar post (among other resources I've looked at) suggested that none of the units (FPKM, RPKM, TPM) are suitable for cross sample comparisons. I am quite confused. Can someone clear up the confusion?

# Answer
> 3 votes
Neither of RPKM, FPKM or TPM is a good choice. All these methods are similar, they only perform the operations in a slightly different order. Eventually they correct for gene length and sequencing depth. Correction for sequencing depth is necessary but not sufficient in many cases. Correction for length is typically not done for a differential analysis as it reduces the counts of many genes and therefore throws away statistical power. I will not describe this in detail since the StatQuest video series you link includes videos about DESeq2 and edgeR normalization procedures which extensively cover the normalization procedures and why these approaches are superior to RPKM/FPKM/TPM. In short: You have to correct for differences in library composition to compare between samples. Please watch these videos and read the manuals of the established RNA-seq analysis packages such as but not exclusively limma, edgeR and DESeq2.

---
Tags: differential-expression, gene-expression, normalization
---


# Question
Title: Is it conflicting: significant difference according to Wald test but not in PERMANOVA and ANOSIM to compare species abundances of metagenomic samples?
If there is no significant difference between two groups of metagenomic samples according to "Multivariate differential abundance tests" for example PERMANOVA and ANOSIM, does it makes sense to use "Univariate differential abundance testing" to see, for example, what species have significantly different abundances between groups?

The issue is that I found species that have significantly different abundances among samples of two groups using Wald test in DESeq2, but there is no significant difference between the two groups according to PERMANOVA and ANOSIM.

# Answer
> 0 votes
Yes, it does make sense to perform differential abundance testing after one doesn't find compositional differences with PERMANOVA / ANOSIM.

It seems to me you are drawing an (implicit) analogy to ANOVA and post-hoc tests like the Tukey Honest Significant Difference test. Intuitively, the analogy makes sense, however, these situations are quite different. PERMANOVA / ANOSIM and DESeq2 are operating on different aspects of the data, with different transformations and under different assumptions.

PERMANOVA and ANOSIM do not operate on the raw counts directly, rather, they test differences between groups by calculating some distance or dissimilarity between the samples. One can choose between several distance / dissimilarity measures, and the test result may be impacted by the chosen measure.

DESeq2 (and other packages), on the other hand, test for species abundances under a parametric model, operating on the raw counts directly.

Thus (sadly) one is not a natural, optional follow-up of the other, they are complementary and will tell different stories about the data.

---
Tags: differential-expression, deseq2, metagenome
---


# Question
Title: Pre-filtering genes for Principal Component Analysis
I have a raw counts data-set of 20,502 genes and 137 samples. I want to find out Principal Components which best explain variation between samples in different stages of tumor.

I am new to Machine learning and would like some help in selecting a criteria to reduce the number of genes before doing a PCA.

I have tried doing a DESeq2 analysis on the data-set but I am confused as to what to do next.

# Answer
> 3 votes
You can filter out genes that do not show variation across samples, they would not differentiate samples anyway. For the specifics on how to do so, please see the `Subpopulation Analysis` section of a nice single cell RNA-seq workflow from the Kharchenko Lab. I believe there won't be fundamental differences when applying this approach to bulk RNA-seq.

# Answer
> 0 votes
I think "the best" way to do this is going to depend on what your question is, unfortunately.

For example, let's say you're trying to find which genes vary by tumour stage. Variance of a given gene might be high across all samples, in which case it has no relation to the stage and might be useless for you. Alternatively, you may have a lowly-to-moderately expressed gene that does not vary much between samples of the same stage, but differs between stages, which is useful for you. However, because it's not highly-expressed, its variance may be small compared to the variance of highly-expressed genes and skipped over.

But if you're trying to do sequencing batch correction to identify batch effects, looking at raw sequencing counts and doing some type of quantile normalization might be all you need.

So I don't think there's a one-size-fits-all answer to your question, unless you're able to refine it a little bit

---
Tags: rna-seq, pca
---


# Question
Title: Calculating Synonymous and non-synonymous
```
library(seqinr)
library(ape)
library(phangorn)

```

alignment file data file sylvia.aln \<- read.alignment("abhi\_seq/all\_seq\_no\_gal.fas", format="fasta") sylvia.aln

```
kaks(sylvia.aln, verbose=TRUE)

Warning message:
In kaks(sylvia.aln) : sequence lengths are not a multiple of 3

```

It ends with this warning and there are no result.

When i ran one of the seqinr example data what I see there are no gaps in the sequence files. If i understand the results is due to the gap I suppose.

Any help or suggestion would be really appreciated.

# Answer
Woaa that is complicated data.

The reason it fails is easy, correcting it is harder. It fails because the open reading frame is lost during in silico translation, e.g. several gaps have 4 indels and the final nucleotides after the final indel are always 4. Thus it will invariably encounter a stop codon OR it will simply refuse to deal with an indel if it is incorporated into a triplet codon. There is no triplet codon model for an indel so the whole calculations stops.

Correcting data that complex ... the easiest way to do it is align you entire data set using amino acids, e.g. muscle. Have an identical nucleotide data set. MAKE SURE the names and positions are identical between your amino acid and nucleotide data, as a programmer this isn't hard to do however, everyone has their own scripts to do this. If you don't do this you next step will fall over.

Use a amino acid to nucleotide alignment translator. I use `transalign` (I think its EMBOSS). This will give you an inframe alignment. I noted many/all of you start codons are ATG which is a good sign. Check the nucleotide output and make sure there is no overhangs of the triplet codon at the terminal 3' end and check in `Jalview` or whatever your alignment editor of choice is that the translation is okay.

This will alignment work in the synonymous/non-synonymous analysis and all subsequent/complimentary analyses.

If you have an intron in protein your you'll need to remove it (hopefully not).

Finally you need to be careful of '3rd codon saturation' in your data because that heavily affects your trees and its certainly a risk given the heterogeneity of the indel structure, but first I would get an inframe alignment. There are checks for saturation (I'm sure even R can do this).
> 1 votes

---
Tags: r, phylogenetics
---


# Question
Title: Is outlier detection the right choice for determining if a value lies significantly away from the mean in a population in gene frequency?
I have the following dataframe that is comprised of enzyme substrates as the variable and then the unique number of times enzymes capiable of degrading this substrates appeared in my genomes and then,total count of times all of those enzymes appeard in my samples. For example, a 7 unique enzymes were found to be capable of degrading chitin, and they occured a total of 298 times across my genomes.

I wish to see if any particular substrate has significantly larger occurence in both unique and total, than the rest of the substrates. Would outlier detection be appropriate for this or would a simple T-test do?

# Answer
> 0 votes
I agree you have the right idea because that data doesn't look at all linear, so regression analysis looks really complicated. Your approach a straight comparison keeps it simple and is within the biological question you are asking,

However, I would use a non-parametric such as Wilcoxon instead of a T-test. I don't fully know the null distribution you are comparing it against. Wilcoxan doesn't require an assumption of normality. T-tests invariably require a transformation to normality and that gets complicated, but a log or sqrt transformation might be okay. If Wilcoxan produces a no significant difference you report it and move onto a T-test.

It is difficult to use outlier analysis with 9(?) observations, thats used if you have a way-out result and can afford to exclude it (\>3 standard deviations). Outlier analysis is to remove a single result that is messing up your analysis and doesn't affect the bulk of the data. You alot of observations to be considering outlier analysis.

---
Tags: statistics, outlier, chi-square
---


# Question
Title: Which module to select for Pathway analysis based on Module trait correlation and pvalue?
I have a total of 35 tumor samples classified into 4 subtypes. Subtype A, B, C, and D. I have RNAseq data. I'm interested in identifying modules related to each subtype with co-expression network analysis. I have read the WGCNA tutorial and create a `trait.csv` file like below:

Showing here some rows of the `trait.csv` file.

```
SamInfo <- read.csv("trait.csv")

Samples = rownames(datExpr);
traitRows = match(Samples, SamInfo$Samples);
datTraits = SamInfo[traitRows, -1];
rownames(datTraits) = SamInfo[traitRows, 1];

nGenes = ncol(datExpr);
#15648L

nSamples = nrow(datExpr);
#35L

#calculate eigengenes (1st principal component) of modules
MEs0 = moduleEigengenes(datExpr, dynamicColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, datTraits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);

#Will display correlations and their p-values
textMatrix = paste(signif(moduleTraitCor, 2), "\n(", signif(moduleTraitPvalue, 1), ")",
                   sep = "");
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3));

#Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor, xLabels = names(datTraits),
               yLabels = names(MEs), ySymbols = names(MEs), colorLabels = FALSE,
               colors = blueWhiteRed(50), textMatrix = textMatrix, setStdMargins = FALSE,
               cex.text = 0.5, zlim = c(-1,1), main = paste("Module-trait relationships"))

```

The Correlation between Subtypes and Modules Heatmap looks like below:

For functions identification (Pathway Analysis), from the above output for each subtype, should I select only the modules that were significantly positively correlated with subtypes, or Can I also use the significantly negatively correlated modules?

Is it good to use the modules that were not significant for pathway analysis?

In case if I can use positive and negative correlated modules, Can I combine those modules and use those genes for Pathway analysis?

# Answer
> 1 votes
If you can run pathway analysis programmatically on all modules, yes, do it on all modules. Otherwise, select the modules with strongest associations with each subtype, both positive and negative. I would keep the modules separate because they may have very different biological functions.

---
Tags: r, rna-seq, pathway, wgcna, coexpression
---


# Question
Title: Error: FASTA file not found using building kallisto index
I am trying to build an index on kallisto on mac using R studio and have saved my FASTA file to my desktop, however, when I run `kallisto index -i DI_mRNA.index DI_mRNA.fa.gz`, there is an error saying `Error: FASTA file not found DI_mRNA.fa.gz`

# Answer
> 2 votes
As mentioned in a comment, your working directory lacks the fasta file.

---
Tags: kallisto
---


# Question
Title: Clarification about how mapping to a reference genome works in RNA seq
When mapping reads to a reference genome, how is it possible to tell which part of the genome the read is referencing if the genome possibly has the same repeated sequences in different areas? What I mean by this is let's say for the sake of argument that you have a read "AGCTAGCT". If the reference genome has two places where the sequence "AGCTAGCT" is present, how do we tell which one it's coming from? I'm assuming that the length of the read is chosen long enough such that this doesn't happen but perhaps I'm wrong?

# Answer
As the other answers/comments have covered, some reads will inevitably align to multiple positions in the genome. This can come from a variety of sources, like aligning to repetitive regions of the genome, short read fragments (that will likely match multiple regions in the genome, even if it, in theory, came from a single location), and low base call qualities (the aligner doesn't have lots of confidence aligning it it to one location because it could map to another, with some likely mismatched calls in the read).

This isn't inherently a problem, depending on the analysis you want to do after aligning. In fact, having a single read that is compatible with multiple transcripts is one of the ideas behind the reference index used in kallisto, and the "transcript compatibility counts".

If you want to minimize this effect, two easy ways are to create longer reads on your sequencer, or to use paired-end reads (in theory, both ends should align relatively close to each other, so if one of the mates maps to multiple places, you can borrow from the confidence of the alignment of its mate).
> 1 votes

# Answer
You can't always know exactly where every read originated, even if you use long reads to sequence full-length transcripts (odds are there's at least one transcript that's duplicated in the genome). This is also not exactly a problem, such reads can either be discarded or reported as multimapping.
> 1 votes

---
Tags: rna-seq, reference-genome, reads
---


# Question
Title: merge rows (to columns) in R dataframe based on IDs
I have a dataframe with two columns: `$1 = transcript ID`, `$2 = expression counts (tpm)`. This comes from a merging of several kallisto counts (if it matters).

I would like to build, in the most basic R way (no tidyverse), a dataframe with as column $1 the unique transcript IDs, and the expression countsas following columns (as many as needed).

Sample input:

```
$ head test1
G01000001_1 0
G01000001_1 0.00855717
G01000001_1 0.0113467
G01000001_1 0.0121086
G01000001_1 0.0121108
G01000001_1 0.0148402
G01000001_1 0.0183057
G01000001_1 0.0268394
G01000001_1 0.0298547

```

Sample output:

```
$ head -3 test2

G01000001_1 0  0.00855717 0.0113467 0.0121086 0.0121108 0.0148402 0.0183057 0.0268394 0.0298547 0.0444686 0.046795 0.0547494 0.0640871 0.0674897 0.0719834 0.082262 0.0867695 0.094905 0.10269 0.11488 0.121549 0.131143 0.141432 0.145758 0.167886 0.184162 0.192757 0.200921 0.249595 0.283296 0.430173 0.473344 0.518981 0.674573 0.701431 0.845816 2.77479 6.28394
G01000002_1 0  0.0277237 0.0294273 0.029438 0.0379258 0.0606153 0.0614727 0.069608 0.142541 0.197644 0.218152 0.22193 0.227401 0.230179 0.259295 0.279463 0.289982 0.292763 0.303268 0.305259 0.308236 0.319988 0.532091 0.632095 0.707766 0.746436 0.821887 0.822496 0.997334 1.06616 1.21752 1.48436 10.0862 2.01763 3.94584 5.02089 7.55767 8.3501
G01000003_1 0  0.00265959 0.00322424 0.00920435 0.454749 0.692183 0.709159 0.948286 1.11827 1.39854 1.39985 1.5781 1.86343 10.2554 10.7185 11.0315 14.342 15.2383 2.13068 2.41792 2.52833 2.77868 3.56981 3.84991 3.97521 4.10445 4.16344 4.26902 4.97202 47.8293 47.9114 5.21789 5.21926 5.35712 5.49942 5.7085 5.87654 50.7484 51.3721 56.7279 7.35062 7.54485 7.64363e-05 8.1367 8.34828

```

This could be achieved in awk, but I would really like an R solution. For the logic, here is the awk way:

```
sort test1 | awk '{if(a!=$1) {a=$1; printf "\n%s%s",$0,FS} else {a=$1;$1="";printf $0 }} END {printf "\n" }' > test2

```

# Answer
For a data.frame named `foo`:

```
sink('some_output_file.txt')
uvals = unique(foo[,1])
for(uIDX in c(1:length(uvals))) {
    gname = uvals[uIDX]
    IDX = which(foo[,1] == gname)
    cat(sprintf("%s\t%s\n", gname, paste0(foo[IDX,2], collapse="\t")))
}
sink()

```
> 1 votes

# Answer
For that purpose, you may also use the **reshape2** R package which allows you to switch wide to long format for data.frames. However, for that, you need an additional column with sample identifiers, e.g.:

```
G01000001_1 Sample1 0
G01000001_1 Sample2 0.00855717
G01000001_1 Sample3 0.0113467
G01000002_1 Sample1 0.0121086
G01000002_1 Sample2 0.0121108
G01000002_1 Sample3 0.0148402
G01000003_1 Sample1 0.0183057
G01000003_1 Sample2 0.0268394
G01000003_1 Sample3 0.0298547

```

Then, load and convert the data in R:

```
library(reshape2)
a<-read.table("test")
dcast(a, V1~V2, value.var = "V3")

```
> 1 votes

---
Tags: r, text-processing, merge
---


# Question
Title: ERROR: (subscript) logical subscript too long
I am trying to remove a row of a dataframe in R however because the dataframe is so large when I use:

```
rawCounts <- rawCounts[rawCounts != geneID, ]

```

I get the error:

```
ERROR in rawCounts[rawCounts != geneID, ]: (subscript) logical subscript too long

```

Is there another way of removing a specific row in a large dataframe that is not so memory intensive?

Thanks in advance.

I have quite a bit of memory, memory.limit() returns:

```
> memory.limit()
[1] 16318

```

# Answer
I am answering my own question here.

I realised I should have been using

```
rawCounts <- rawCounts[rownames(rawCounts) != geneID, ]

```

Problem solved.
> 2 votes

---
Tags: r
---